text:
# -*- coding: utf-8 -*-
# # Runge-Kutta methods

using Plots
gr()

# We will implement in Julia different numerical methods to solve
#
# $$
# y'(t) = 1 - y(t)
# $$
#
# $ t \in [0,5] \qquad \mbox{ and } \qquad y(0) = 0 $

# ## Explicit Euler

"""

   euler(f::Function, t::Float64, y::Float64, h::Float64)

explicit euler method function that returns

```math
y^{n+1} = y^n + h \\cdot f(t^n, y^n)
```

"""
function euler(f::Function, t::Float64, y::Float64, h::Float64)
    t + h, y + h * f(t,y)
end
#----------------------------------------------------------------------------

?euler

# ## Runge-Kutta 2nd order

"""

   rk2(f::Function, t::Float64, y::Float64,  dt::Float64)

Runge-Kutta second order method function

```jldoctest
julia> f(t,y) = 1 - y
julia> rk2(f, 0.0, 0.0, 1.0)
(1.0, 0.5)
```

"""
function rk2(f::Function, t::Float64, y::Float64,  h::Float64)
    ỹ = y + h/2 * f(t,y)
    t + h, y + h * f(t+h/2,ỹ)
end
#----------------------------------------------------------------------------

?rk2

# ## Runge-Kutta 4th order

"""

   rk4(f::Function, t::Float64, y::Float64,  dt::Float64)

Runge-Kutta fourth order method function

[Runge–Kutta methods on Wikipedia](https://en.wikipedia.org/wiki/Runge–Kutta_methods)

"""
function rk4(f::Function, t::Float64, y::Float64,  dt::Float64)

    y₁ = dt * f(t,y)
    y₂ = dt * f(t+dt/2,y+y₁/2)
    y₃ = dt * f(t+dt/2,y+y₂/2)
    y₄ = dt * f(t+dt,y+y₃)

    t+dt, y+(y₁+2*y₂+2*y₃+y₄)/6

end
#----------------------------------------------------------------------------

?rk4
#----------------------------------------------------------------------------

# ## Solver function

"""

    solver(f::Function, Method::Function, t₀::Float64,
                y₀::Float64, dt::Float64, nsteps::Int64)

Solve numerically the equation ``y' = f(t, y)``

with `y(t₀)= y₀` and `nsteps` steps `h`

## Arguments
- `f::Function`: the function `f` of equation ``y' = f(t,y)``.
- `Method::Function`: numerical method from (tⁿ,yⁿ) returns ``(t^{n+1},y^{n+1})``

"""
function solver(f::Function,
                Method::Function,
                t₀::Float64,
                y₀::Float64, h::Float64, nsteps::Int64)

    t = zeros(Float64,nsteps)
    y = similar(t)

    t[1] = t₀
    y[1] = y₀

    for i in 2:nsteps
       t[i], y[i] = Method(f,t[i-1],y[i-1], h)
    end

    t, y

end
#----------------------------------------------------------------------------

?solver
#----------------------------------------------------------------------------

nsteps, tfinal   = 7, 5.0
t₀, x₀ = 0., 0.

dt = tfinal / (nsteps-1)
f(t, x) = 1 - x
#----------------------------------------------------------------------------

plot( solver(f, euler, t₀, x₀, dt, nsteps); marker = :o, label=:euler)
plot!(solver(f, rk2,   t₀, x₀, dt, nsteps); marker = :d, label=:rk2)
plot!(solver(f, rk4,   t₀, x₀, dt, nsteps); marker = :p, label=:rk4)
t = 0:0.1:5
plot!(t, 1 .- exp.(-t); line = 3, label = :exact)
#----------------------------------------------------------------------------

# # Callable object
#
# We want to build a numerical method to solve the problem and then use it
# as a function and still call the solver function.

#
#
# Runge-Kutta scheme can be built using Butcher tableau :
#
# \begin{array}{c|cccc}
#  c₁  &        &        &        &        \\
#  c₂  & a_{21} &        &        &        \\
#  c₃  & a_{31} & a_{32} &        &        \\
#  c₄  & a_{41} & a_{42} & a_{43} &        \\
#  \hline
#      & b_1    & b_2    & b_3    & b_4    \\
# \end{array}
#
# $$
# \forall i = 1, \dotsc, q, \begin{cases}t_{n,i} &= t_n + c_i h_n, \\
# y_{n,i} &= y_n + h_n \sum_{k = 1}^{i-1} a_{ik} p_{n,k}\\
# p_{n,i} &= f(t_{n,i}, y_{n,i}) \end{cases}
# $$
#
# $$
# y_{n+1} = y_n + h_n \sum_{k = 1}^q b_k p_{n,k}.
# $$
#

mutable struct RungeKutta

    q :: Int64
    a :: Array{Float64, 2}
    b :: Array{Float64, 1}
    c :: Array{Float64, 1}

    tn :: Vector{Float64}
    yn :: Vector{Float64}
    pn :: Vector{Float64}

    function RungeKutta( a::Array{Float64,2}, b::Vector{Float64}, c::Vector{Float64})

        q = length(c)
        @assert ( length(c) == size(a,1))
        @assert ( length(b) == size(a,2))
        tn = zeros(Float64, q)
        yn = zeros(Float64, q)
        pn = zeros(Float64, q)
        new( q, a, b, c, tn, yn, pn)

    end

end

function (rk::RungeKutta)(f::Function, t::Float64, y::Float64,  h::Float64)

    for i = 1:rk.q
        rk.tn[i] = t + rk.c[i] * h
        rk.yn[i] = y + h * sum([rk.a[i,k]*rk.pn[k] for k = 1:i-1])
        rk.pn[i] = f(rk.tn[i],rk.yn[i])
    end

    t + h, y + h * sum([rk.b[k]*rk.pn[k] for k in 1:rk.q ])

end
#----------------------------------------------------------------------------

function solver(f::Function,
                Method::RungeKutta,
                t₀::Float64,
                y₀::Float64, h::Float64, nsteps::Int64)

    t = zeros(Float64,nsteps)
    y = similar(t)

    t[1] = t₀
    y[1] = y₀

    for i in 2:nsteps
       t[i], y[i] = Method(f,t[i-1],y[i-1], h)
    end

    t, y

end
#----------------------------------------------------------------------------

a = [ 0   0   0 0;
      1/2 0   0 0;
      0   1/2 0 0;
      0   0   1 0]

b = [1/6 ,1/3, 1/3, 1/6]
c = [0   ,1/2, 1/2, 1  ]

rk4_new = RungeKutta(a, b, c)
#----------------------------------------------------------------------------

t = 0:0.1:5
plot(t, 1 .- exp.(-t),label = :exact)
plot!(solver(f, rk4_new, t₀, x₀, dt, nsteps), marker = :o,  label=:rk4_new)
#----------------------------------------------------------------------------

a = [  0   0  0  0;
     1/3   0  0  0;
    -1/3   1  0  0;
       1  -1  1  0]

b = [1/8, 3/8, 3/8, 1/8]
c = [0, 1/3,  2/3, 1]
rk4_38 = RungeKutta(a, b, c)
plot!(solver(f, rk4_38, t₀, x₀, dt, nsteps), marker = :r,  label=:rk4_38)
#----------------------------------------------------------------------------

# [Reference](https://www.juliabloggers.com/runge-kutta-methods/):Ordinary Differential Equation Solvers: Runge-Kutta Methods by Christina Lee

# # Creating expressions via interpolation
#
# Create a solver function with the method choosen at initialization.

macro add(x, y)
    return :($x + $y)
end

@add 2 3      ## or @add(2, 3)
#----------------------------------------------------------------------------

macro abs(x)
    return :( $x > 0 ? $x : -$x)
end
@abs(-2), @abs(2)
#----------------------------------------------------------------------------

macro make_method( meth)
    return quote
        function (f::Function, t₀::Float64,
                  x₀::Float64, dt::Float64, nsteps::Int64)

            t = zeros(Float64,nsteps)
            x = zeros(Float64,nsteps)

            t[1] = t₀
            x[1] = x₀

            for i in 2:nsteps
               t[i], x[i] = $meth(f,t[i-1],x[i-1], dt)
            end

            return t, x
        end
    end
end

rk4_solver = @make_method rk4

plot(rk4_solver(f, t₀, x₀, dt, nsteps))
#----------------------------------------------------------------------------

# [DifferentialEquations.jl](http://docs.juliadiffeq.org/latest/)
#
# [ODE Solvers](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#OrdinaryDiffEq.jl-1)

using OrdinaryDiffEq
using Plots

f(y,p,t) = 1.0 - y
y₀ = 0.0
t  = (0.0,5.0)
prob = ODEProblem(f,y₀,t)
euler  = solve(prob,Euler(), dt=1.0)
rk4  = solve(prob, RK4(), dt=1.0)
plot(euler,label="Euler")
plot!(rk4,label="RK4")
plot!(1:0.1:5, t->1. - exp(-t),lw=3,ls=:dash,label="True Solution!")
#----------------------------------------------------------------------------

using DifferentialEquations

function lorenz(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
#----------------------------------------------------------------------------

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob)
#----------------------------------------------------------------------------

plot(sol,vars=(1,2,3))

# ----------------------------------------------------------------------------

max_stars_repo_name: pnavaro/math-julia

max_issues_repo_name: pnavaro/math-julia

max_forks_repo_name: pnavaro/math-julia


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# for generating simulation objects based on a config file

mutable struct GenConfig
	inputPath::String
	outputPath::String
	mode::String # "all" or "calls"
	numCallsFiles::Int

	# output file names
	ambsFilename::String
	arcsFilename::String
	callsFilename::String
	hospitalsFilename::String
	mapFilename::String
	nodesFilename::String
	prioritiesFilename::String
	stationsFilename::String
	travelFilename::String

	# counts
	numAmbs::Int
	numCalls::Int
	numHospitals::Int
	numStations::Int

	# graph
	xNodes::Int # number of nodes in x direction
	yNodes::Int # number of nodes in y direction

	# map
	map::Map
	mapTrim::Float # fraction of map border to trim, to make sure objects generated on map are inside borders

	# misc
	startTime::Float
	maxCallArrivalTime::Float
	minLastCallArrivalTime::Float
	targetResponseDuration::Float
	offRoadSpeed::Float
	stationCapacity::Int
	travelModeSpeeds::Vector{Float}

	# call density raster
	callDensityRasterFilename::String
	cropRaster::Bool
	callRasterCellSeed::Int # seed for rng, will generate raster cell index
	callRasterCellLocSeed::Int # seed for rng, will generate location within raster cell

	# call related distributions and random number generators
	interarrivalTimeDistrRng::DistrRng
	priorityDistrRng::DistrRng
	dispatchDelayDistrRng::DistrRng
	onSceneDurationDistrRng::DistrRng
	transportDistrRng::DistrRng
	handoverDurationDistrRng::DistrRng

	# misc RNGs
	ambStationRng::AbstractRNG
	callLocRng::AbstractRNG
	hospitalLocRng::AbstractRNG
	stationLocRng::AbstractRNG

	# undefined
	travelTimeFactorDistrRng::DistrRng
	callRasterSampler::RasterSampler

	GenConfig() = new("", "", "", 1,
		"", "", "", "", "", "", "", "", "",
		nullIndex, nullIndex, nullIndex, nullIndex,
		nullIndex, nullIndex,
		Map(), 1e-6,
		nullTime, nullTime, nullTime, nullTime, nullTime, nullIndex, [],
		"", false, nullIndex, nullIndex)
end

function readGenConfig(genConfigFilename::String)::GenConfig
	# read gen config xml file
	genConfigFilename = genConfigFilename |> interpolateString |> abspath
	global genConfigFileDir = dirname(genConfigFilename)
	rootElt = xmlFileRoot(genConfigFilename)
	@assert(xName(rootElt) == "genConfig", string("xml root has incorrect name: ", xName(rootElt)))

	genConfig = GenConfig()

	genConfig.inputPath = containsElt(rootElt, "inputPath") ? joinPathIfNotAbs(genConfigFileDir, eltContentInterpVal(rootElt, "inputPath")) : genConfigFileDir
	genConfig.outputPath = joinPathIfNotAbs(genConfigFileDir, eltContentInterpVal(rootElt, "outputPath")) # output path can be absolute, or relative to genConfigFileDir
	genConfig.mode = eltContent(rootElt, "mode")
	@assert(in(genConfig.mode, ("all", "calls")))

	if genConfig.mode == "all"
		readGenConfigAll!(genConfig, rootElt)
	elseif genConfig.mode == "calls"
		readGenConfigCalls!(genConfig, rootElt)
	end

	return genConfig
end

# read xml rootElt to set fields of genConfig for generation mode 'all'
# mutates: genConfig
function readGenConfigAll!(genConfig::GenConfig, rootElt::XMLElement)
	readGenConfigMap!(genConfig, rootElt)
	readGenConfigCalls!(genConfig, rootElt)

	# output filenames
	simFilesElt = findElt(rootElt, "simFiles")
	simFilePath(filename::String) = joinPathIfNotAbs(genConfig.outputPath, eltContentInterpVal(simFilesElt, filename)) # filename can be absolute, or relative to output path
	genConfig.ambsFilename = simFilePath("ambulances")
	genConfig.arcsFilename = simFilePath("arcs")
	genConfig.hospitalsFilename = simFilePath("hospitals")
	genConfig.mapFilename = simFilePath("map")
	genConfig.nodesFilename = simFilePath("nodes")
	genConfig.prioritiesFilename = simFilePath("priorities")
	genConfig.stationsFilename = simFilePath("stations")
	genConfig.travelFilename = simFilePath("travel")
	@assert(genConfig.callsFilename != "") # should already be set in readGenConfigCalls!()

	# read sim parameters
	simElt = findElt(rootElt, "sim")

	# number of ambulances, hospitals, stations
	genConfig.numAmbs = eltContentVal(simElt, "numAmbs")
	genConfig.numHospitals = eltContentVal(simElt, "numHospitals")
	genConfig.numStations = eltContentVal(simElt, "numStations")

	# number of nodes in x and y direction for grid shaped graph
	graphElt = findElt(simElt, "graph")
	genConfig.xNodes = eltContentVal(graphElt, "xNodes")
	genConfig.yNodes = eltContentVal(graphElt, "yNodes")

	# misc values
	targetResponseDurationString = containsElt(simElt, "targetResponseDuration") ? "targetResponseDuration" : "targetResponseTime" # compat for "targetResponseTime"
	genConfig.targetResponseDuration = eltContentVal(simElt, targetResponseDurationString)
	genConfig.offRoadSpeed = eltContentVal(simElt, "offRoadSpeed") # km / day
	genConfig.stationCapacity = eltContentVal(simElt, "stationCapacity")
	travelModeSpeedsElt = findElt(simElt, "travelModeSpeeds")
	genConfig.travelModeSpeeds = (travelModeSpeedsElt == nothing ? [1.5 * genConfig.offRoadSpeed] : eltContentVal(travelModeSpeedsElt))

	# some defaults - should move to config file sometime
	genConfig.ambStationRng = MersenneTwister(0)
	genConfig.hospitalLocRng = MersenneTwister(4)
	genConfig.stationLocRng = MersenneTwister(5)
	genConfig.travelTimeFactorDistrRng = DistrRng(Uniform(1.0, 1.1); seed = 99)
end

# read xml rootElt to set fields of genConfig for generation mode 'calls'
# mutates: genConfig
function readGenConfigCalls!(genConfig::GenConfig, rootElt::XMLElement)
	genConfig.numCallsFiles = containsElt(rootElt, "numCallsFiles") ? eltContentVal(rootElt, "numCallsFiles") : 1
	@assert(genConfig.numCallsFiles >= 1)

	# output filename
	callsFilename = "calls.csv"
	if containsElt(rootElt, "simFiles")
		simFilesElt = findElt(rootElt, "simFiles")
		if containsElt(simFilesElt, "calls")
			callsFilename = eltContentInterpVal(simFilesElt, "calls")
		end
	end
	genConfig.callsFilename = joinPathIfNotAbs(genConfig.outputPath, callsFilename) # filename can be absolute, or relative to output path

	# read sim call parameters
	simElt = findElt(rootElt, "sim")

	# call distributions and random number generators
	callDistrsElt = findElt(simElt, "callDistributions")
	function callDistrsEltContent(distrName::String)
		distrElt = findElt(callDistrsElt, distrName)
		distr = eltContentVal(distrElt)
		seedAttr = attribute(distrElt, "seed")
		seed = (seedAttr == nothing ? nullIndex : eval(Meta.parse(seedAttr)))
		return DistrRng(distr; seed = seed)
	end
	genConfig.interarrivalTimeDistrRng = callDistrsEltContent("interarrivalTime")
	genConfig.priorityDistrRng = callDistrsEltContent("priority")
	genConfig.dispatchDelayDistrRng = callDistrsEltContent("dispatchDelay")
	genConfig.onSceneDurationDistrRng = callDistrsEltContent("onSceneDuration")
	if containsElt(callDistrsElt, "transport") && containsElt(callDistrsElt, "handoverDuration")
		genConfig.transportDistrRng = callDistrsEltContent("transport")
		genConfig.handoverDurationDistrRng = callDistrsEltContent("handoverDuration")
	else # compat
		@assert(containsElt(callDistrsElt, "transfer") && containsElt(callDistrsElt, "transferDuration"))
		@warn("`transfer` and `transferDuration` elements for call distributions in gen config are deprecated, use `transport` and `handoverDuration` instead.")
		genConfig.transportDistrRng = callDistrsEltContent("transfer")
		genConfig.handoverDurationDistrRng = callDistrsEltContent("transferDuration")
	end

	# call gen duration / count
	genConfig.startTime = eltContentVal(simElt, "startTime")
	@assert(genConfig.startTime >= 0)
	genConfig.numCalls = containsElt(simElt, "numCalls") ? eltContentVal(simElt, "numCalls") : nullIndex # can alternatively specify maxCallArrivalTime or minLastCallArrivalTime
	genConfig.maxCallArrivalTime = containsElt(simElt, "maxCallArrivalTime") ? eltContentVal(simElt, "maxCallArrivalTime") : nullTime
	genConfig.minLastCallArrivalTime = containsElt(simElt, "minLastCallArrivalTime") ? eltContentVal(simElt, "minLastCallArrivalTime") : nullTime
	@assert((genConfig.numCalls != nullIndex) + (genConfig.maxCallArrivalTime != nullTime) + (genConfig.minLastCallArrivalTime != nullTime) == 1, "Need exactly one of these values: numCalls, maxCallArrivalTime, minLastCallArrivalTime.")
	@assert(genConfig.startTime <= genConfig.maxCallArrivalTime || genConfig.maxCallArrivalTime == nullTime)
	@assert(genConfig.startTime <= genConfig.minLastCallArrivalTime || genConfig.minLastCallArrivalTime == nullTime)

	# call density raster
	callDensityRasterElt = findElt(simElt, "callDensityRaster")
	genConfig.callDensityRasterFilename = joinPathIfNotAbs(genConfig.inputPath, eltContentInterpVal(callDensityRasterElt, "filename"))
	genConfig.cropRaster = containsElt(callDensityRasterElt, "cropRaster") ? eltContentVal(callDensityRasterElt, "cropRaster") : false
	# seeds
	function callRasterSeedVal(seedName::String)
		seedAttr = attribute(callDensityRasterElt, seedName)
		return seedAttr == nothing ? nullIndex : eval(Meta.parse(seedAttr))
	end
	genConfig.callRasterCellSeed = callRasterSeedVal("cellSeed")
	genConfig.callRasterCellLocSeed = callRasterSeedVal("cellLocSeed")
	if containsElt(simElt, "map") readGenConfigMap!(genConfig, rootElt) end
	if isfile(genConfig.callDensityRasterFilename)
		raster = readRasterFile(genConfig.callDensityRasterFilename)
		if genConfig.cropRaster
			# crop raster to be within genConfig.map
			@assert(genConfig.map.xRange != nullDist) # fails if map has not been set
			mapTrimmed = trimmedMap(genConfig.map, genConfig.mapTrim)
			cropRaster!(raster, mapTrimmed)
		end
		genConfig.callRasterSampler = RasterSampler(raster, genConfig.callRasterCellSeed, genConfig.callRasterCellLocSeed)
	end

	# some defaults - should move to config file sometime
	genConfig.callLocRng = MersenneTwister(1)
end

# mutates: genConfig.map
function readGenConfigMap!(genConfig::GenConfig, rootElt::XMLElement)
	simElt = findElt(rootElt, "sim")
	@assert(containsElt(simElt, "map"))
	mapElt = findElt(simElt, "map")
	map = Map()
	map.xMin = eltContentVal(mapElt, "xMin")
	map.xMax = eltContentVal(mapElt, "xMax")
	map.xScale = eltContentVal(mapElt, "xScale")
	map.xRange = map.xMax - map.xMin
	map.yMin = eltContentVal(mapElt, "yMin")
	map.yMax = eltContentVal(mapElt, "yMax")
	map.yScale = eltContentVal(mapElt, "yScale")
	map.yRange = map.yMax - map.yMin
	@assert(map.xRange > 0 && map.yRange > 0)
	genConfig.map = map
end

function runGenConfig(genConfig::GenConfig; overwriteOutputPath::Bool = false, doPrint::Bool = true)
	outputPath = genConfig.outputPath # shorthand
	if isdir(outputPath) && !isempty(readdir(outputPath))
		if overwriteOutputPath
			doPrint && println("Deleting directory: ", outputPath)
			rm(outputPath; recursive=true)
		else
			error("Directory is not empty: ", outputPath)
		end
	end
	isdir(outputPath) || mkdir(outputPath)

	doPrint && println("Generation mode: ", genConfig.mode)
	if genConfig.mode == "all"
		# make all
		ambulances = makeAmbs(genConfig)
		calls = makeCalls(genConfig)
		hospitals = makeHospitals(genConfig)
		stations = makeStations(genConfig)
		travel = makeTravel(genConfig)

		# grid shaped graph
		graph = LightGraphs.grid([genConfig.xNodes, genConfig.yNodes])
		nodes = makeNodes(genConfig, graph)
		(arcs, travelTimes) = makeArcs(genConfig, graph, nodes)

		# save all
		doPrint && println("Saving output to: ", outputPath)
		writeAmbsFile(genConfig.ambsFilename, ambulances)
		writeArcsFile(genConfig.arcsFilename, arcs, travelTimes, "undirected")
		writeCallsFile(genConfig.callsFilename, genConfig.startTime, calls)
		writeHospitalsFile(genConfig.hospitalsFilename, hospitals)
		writeMapFile(genConfig.mapFilename, genConfig.map)
		writeNodesFile(genConfig.nodesFilename, nodes)
		writePrioritiesFile(genConfig.prioritiesFilename, repeat([genConfig.targetResponseDuration],numPriorities))
		writeStationsFile(genConfig.stationsFilename, stations)
		writeTravelFile(genConfig.travelFilename, travel)
	elseif genConfig.mode == "calls"
		runGenConfigCalls(genConfig, doPrint = doPrint)
	else
		error("Unrecognised generation mode")
	end

	return genConfig
end
function runGenConfig(genConfigFilename::String; overwriteOutputPath::Bool = false, doPrint::Bool = true)
	genConfig = readGenConfig(genConfigFilename)
	return runGenConfig(genConfig, overwriteOutputPath = overwriteOutputPath, doPrint = doPrint)
end

# generate calls, and generate locations using call density raster file (file name stored in genConfig)
function runGenConfigCalls(genConfig::GenConfig; doPrint::Bool = true, writeFile::Bool = true)
	numFiles = genConfig.numCallsFiles # shorthand
	@assert(numFiles >= 1)
	doPrint && println("Generating $numFiles call set(s).")
	callSets = [makeCalls(genConfig) for i = 1:numFiles]

	if writeFile
		if numFiles == 1
			calls = callSets[1]
			filename = genConfig.callsFilename
			doPrint && println("Saving calls file to: ", filename)
			writeCallsFile(filename, genConfig.startTime, calls)
		else
			@assert(numFiles > 1)
			doPrint && println("Saving calls files to: ", genConfig.outputPath)
			for j = 1:numFiles
				filename = joinpath(genConfig.outputPath, "calls_$j.csv")
				writeCallsFile(filename, genConfig.startTime, callSets[j])
			end
		end
	end

	return callSets
end

function makeAmbs(genConfig::GenConfig)
	ambulances = Vector{Ambulance}(undef, genConfig.numAmbs)
	for i = 1:genConfig.numAmbs
		ambulances[i] = Ambulance()
		ambulances[i].index = i
		ambulances[i].stationIndex = rand(genConfig.ambStationRng, 1:genConfig.numStations)
		ambulances[i].class = als
	end
	return ambulances
end

function makeArcs(genConfig::GenConfig, graph::LightGraphs.Graph, nodes::Vector{Node})
	numTravelModes = length(genConfig.travelModeSpeeds)
	@assert(numTravelModes >= 1)

	arcs = Vector{Arc}(undef, graph.ne)
	travelTimes = Array{Float,2}(undef,numTravelModes,length(arcs))

	i = 1
	for edge in LightGraphs.edges(graph)
		arcs[i] = Arc()
		arcs[i].index = i
		arcs[i].fromNodeIndex = edge.src
		arcs[i].toNodeIndex = edge.dst
		i = i + 1
	end

	for i = 1:numTravelModes, arc in arcs
		arc.distance = normDist(genConfig.map, nodes[arc.fromNodeIndex].location, nodes[arc.toNodeIndex].location)
		travelTimes[i,arc.index] = arc.distance / genConfig.travelModeSpeeds[i] * rand(genConfig.travelTimeFactorDistrRng)
	end

	return arcs, travelTimes
end

# make calls that are spatially randomly uniform in map, or distributed according to raster
function makeCalls(genConfig::GenConfig; rasterSampler::Union{RasterSampler,Nothing} = nothing)
	if (rasterSampler == nothing && isdefined(genConfig, :callRasterSampler)) rasterSampler = genConfig.callRasterSampler end
	calls = Call[]
	currentTime = genConfig.startTime
	while length(calls) < genConfig.numCalls || genConfig.numCalls == nullIndex
		currentTime += rand(genConfig.interarrivalTimeDistrRng) # apply time step
		if (currentTime > genConfig.maxCallArrivalTime && genConfig.maxCallArrivalTime != nullTime) break end

		call = Call()
		call.index = length(calls) + 1
		call.priority = Priority(rand(genConfig.priorityDistrRng))
		call.arrivalTime = currentTime
		call.dispatchDelay = rand(genConfig.dispatchDelayDistrRng)
		call.onSceneDuration = rand(genConfig.onSceneDurationDistrRng)
		call.transport = (rand(genConfig.transportDistrRng) == 1)
		call.hospitalIndex = nullIndex
		call.handoverDuration = rand(genConfig.handoverDurationDistrRng)
		if !call.transport call.handoverDuration = 0.0 end
		if rasterSampler == nothing
			call.location = randLocation(genConfig.map; trim = genConfig.mapTrim, rng = genConfig.callLocRng)
		else
			call.location = rasterRandLocations(rasterSampler, 1)[1]
		end

		push!(calls, call)

		if (currentTime >= genConfig.minLastCallArrivalTime && genConfig.minLastCallArrivalTime != nullTime) break end
	end

	return calls
end

function makeHospitals(genConfig::GenConfig)
	hospitals = Vector{Hospital}(undef, genConfig.numHospitals)
	for i = 1:genConfig.numHospitals
		hospitals[i] = Hospital()
		hospitals[i].index = i
		hospitals[i].location = randLocation(genConfig.map; trim = genConfig.mapTrim, rng = genConfig.hospitalLocRng)
	end
	return hospitals
end

function makeNodes(genConfig::GenConfig, graph::LightGraphs.Graph)
	nodes = Vector{Node}(undef, LightGraphs.nv(graph)) # should have length xNodes*yNodes
	map = genConfig.map # shorthand
	k = 1
	for j = 1:genConfig.yNodes, i = 1:genConfig.xNodes
		nodes[k] = Node()
		nodes[k].index = k
		nodes[k].location.x = map.xRange * ((i-0.5)/genConfig.xNodes) + map.xMin
		nodes[k].location.y = map.yRange * ((j-0.5)/genConfig.yNodes) + map.yMin
		k += 1
	end

	return nodes
end

function makeStations(genConfig::GenConfig)
	stations = Vector{Station}(undef, genConfig.numStations)
	for i = 1:genConfig.numStations
		stations[i] = Station()
		stations[i].index = i
		stations[i].location = randLocation(genConfig.map; trim = genConfig.mapTrim, rng = genConfig.stationLocRng)
		stations[i].capacity = genConfig.stationCapacity
	end
	return stations
end

function makeTravel(genConfig::GenConfig)
	travel = Travel() # note: will only partially fill travel fields, just enough for writeTravelFile()

	# create single travel mode, add to travel
	travelMode = TravelMode()
	travelMode.index = 1
	travelMode.offRoadSpeed = genConfig.offRoadSpeed
	travel.modes = [travelMode]
	travel.numModes = length(travel.modes)

	# use single travel mode in single travel mode set
	travel.numSets = 1
	travel.modeLookup = repeat([1],1,numPriorities) # modeLookup[setIndex, priority] gives travelModeIndex

	travel.setsStartTimes = [genConfig.startTime]
	travel.setsTimeOrder = [1]

	return travel
end

ext: jl


text:
using Documenter, Poly

makedocs(sitename="Poly.jl Documentation")

deploydocs(
    repo = "github.com/JuliaLabs/Poly.jl",
)


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

using ROOTFramework

using Compat.Test

@testset "TTree I/O" begin
    const scratchdir = mktempdir()
    try
        fname = joinpath(scratchdir, "out.root")

        tfile = open(TFile, fname, "recreate")
        ttree = create_ttree!(tfile, "data", "Data")

        idx = Ref{Int32}(0)
        v = Ref(0.0)

        create_branch!(ttree, "idx", idx)
        create_branch!(ttree, "v", v)

        for i in 1:1000000
            idx.x = i
            v.x = rand()
            push!(ttree)
        end

        close(tfile)
    finally
        rm(scratchdir; recursive = true)
    end
end


text: include("GastonPlot.jl")


text:
hasbackpropagation{F, LD}(::Type{MultilayerPerceptron{F,LD}}) = true

#implementation of the backpropagation algorithm for multilayer perceptrons.
function backpropagate!{F, LD}(mlp::MultilayerPerceptron{F, LD},
                               input::AbstractArray,
                               output::AbstractArray{F},
                               output_deltas::AbstractArray{F},
                               backpropstorage::BackpropStorage{MultilayerPerceptron{F, LD}},
                               input_deltas::VoidableDeltas{F} = nothing)

  l = length(LD) - 1
  #presume that the forward propagation has been done, and that the relevant
  #buffers have been executed.

  this_layer_deltas = output_deltas
  this_layer_output = output

  for l_idx = l:-1:2  #track back through the array.
    this_layer_input_ref = layer(backpropstorage, Val{l_idx})
    next_layer_delta_ref = backprop_layer(backpropstorage, Val{l_idx})

    #call backpropagation on the sublayer, which should be fully connected.
    backpropagate!(mlp.layers[l_idx], this_layer_input_ref, this_layer_output, this_layer_deltas, next_layer_delta_ref)

    #move the references over.
    this_layer_deltas = next_layer_delta_ref
    this_layer_output = this_layer_input_ref
  end

  #solve the last layer, and if the backpropagation isn't supposed to terminate here, pass on the next array.
  backpropagate!(mlp.layers[1], input, this_layer_output, this_layer_deltas, input_deltas)
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
function main()
    println("Executing all benchmarks")
    benchmark_downconvert();
    benchmark_carrier_replica();
    benchmark_code_replica();
    benchmark_correlate();
end

function plotall()
    println("Plotting all the .csv benchmarks under /data ...")
    plot_carrier_replica()
    plot_downconvert()
    plot_code_replica()
    plot_correlate()
end


text:
include("CoreExtended.jl")
include("TimeManager.jl")
include("LoggerManager.jl")

include("RessourceManager.jl")
include("FileManager.jl")
include("Environment.jl")
include("JLScriptManager.jl")
include("WindowManager.jl")

include("MatrixMath.jl")

""" TODO """
module JLGEngine

const ENGINE = JLGEngine

using CoreExtended

dir="JLGEngine/"

#=
dir=dirname(Base.@__FILE__)*"/JLGEngine/"

# search
searchdir(path,key) = filter(x->contains(x,key), readdir(path))
print("search:"); @time fileNames=searchdir(dir,".jl")

# sort
order=createSortedDict(Integer,Array{String,1})
print("sort:"); @time for fileName in fileNames
    path=dir*fileName
    #println("found(\"$fileName\")")
	open(path) do file
		content=readstring(file)
        id=1
         for x in eachmatch(r"using\s+(\w+)\.(\w+)", content)
             #f=x.captures[1]
             #order[id]=f
             #println(x.captures[2])
             id+=1
         end
         if haskey(order,id) push!(order[id],fileName) else order[id]=[fileName] end
	end
    println("-------------------")
    #include(path);
end

println("D: $order")

# include
ex=:()
print("prepare:"); @time for (_,a) in order
    for fileName in a
        println("include(\"$fileName\")")
        ex=:($(ex.args...); include($(dir*fileName)))
        #include(dir*fileName)
    end
end

print("include"); @time eval(ex)
println("---------------");
=#

""" TODO """
abstract type IComponent end

include(dir*"GraphicsManager.jl")
include(dir*"LibGL/LibGL.jl")

include(dir*"Management.jl")
include(dir*"StorageManager.jl")
include(dir*"MeshManager.jl")
include(dir*"ModelManager.jl")
include(dir*"ShaderManager.jl")
include(dir*"TransformManager.jl")
include(dir*"EntityManager.jl")
include(dir*"CameraManager.jl")
include(dir*"TextureManager.jl")
include(dir*"GameObjectManager.jl")
include(dir*"RenderManager.jl")

using .GraphicsManager

# set api references
""" TODO """
function init()
    GraphicsManager.setGraphicsDriver(LibGL)
    ShaderManager.init()
    StorageManager.init()
		TextureManager.init()
end

end #JLGEngine

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using PkgBenchmark

res = benchmarkpkg("TrajectoryOptimization")
export_markdown("benchmark/benchmark_results.md", res)


text:
# This is file provides the ability for the user to override the default prior.

# This code is `include()`'ed inside of `/scripts/venus.jl`, so it will **literarly replace** the
# function lnprior inside `src/model.jl` (all of it!). Be careful to make sure it accepts the same arguments!

using DiskJockey.model
import DiskJockey.model.lnprior_base
import DiskJockey.constants.ModelException

# for example, a replacement for the standard model might look like
function lnprior(pars::ParametersStandard, dpc_mu::Float64, dpc_sig::Float64, grid::Grid)
    lnp = lnprior_base(pars, dpc_mu, dpc_sig)

    # Modify this code here to enforce your prior!

    r_out = grid.Rs[end]/AU # [AU]
    # A somewhat arbitrary cutoff regarding the gridsize to prevent the disk from being too large
    # to fit on the model grid.
    if (3 * pars.r_c) > r_out
        throw(ModelException("User prior violated."))
    end

    return lnp
end

# while an example for the cavity model would look like this
function lnprior(pars::ParametersCavity, dpc_mu::Float64, dpc_sig::Float64, grid::Grid)
    lnp = lnprior_base(pars, dpc_mu, dpc_sig)

    # Modify this code here to enforce your prior!

    r_in = grid.Rs[1]/AU # [AU]
    r_out = grid.Rs[end]/AU # [AU]
    # A somewhat arbitrary cutoff regarding the gridsize to prevent the disk from being too large
    # to fit on the model grid.

    # Also check to make sure that r_cav is less than r_c but larger than r_in.
    if (3 * pars.r_c) > r_out || pars.r_cav < r_in || pars.r_cav > pars.r_c
        throw(ModelException("User prior violated."))
    end

    return lnp

end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using Statistics
using StatsBase: sample
using Distributions
using Random
using CUDA
using Revise
using EvoTrees
using EvoTrees: sigmoid, logit
using Plots

# prepare a dataset
Random.seed!(123)
features = rand(10_000) .* 5
X = reshape(features, (size(features)[1], 1))
Y = sin.(features) .* 0.5 .+ 0.5
Y = logit(Y) + randn(size(Y))
Y = sigmoid(Y)
𝑖 = collect(1:size(X,1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

################################
# linear
################################
params1 = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=200, nbins = 64,
    λ = 0.5, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = fit_evotree_gpu(params1, X_train, Y_train);
@time pred_train_linear = predict_gpu(model, X_train)

x_perm = sortperm(X_train[:,1])
plot(X_train, Y_train, msize = 1, mcolor = "gray", mswidth=0, background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color = "navy", linewidth = 1.5, label = "Linear")
# savefig("figures/regression_sinus_gpu.png")

params1 = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:mse,
    nrounds=200, nbins = 64,
    λ = 0.5, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = fit_evotree_gpu(params1, X_train, Y_train, print_every_n = 25);
@time model = fit_evotree_gpu(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n = 25);
@time pred_train_linear = predict_gpu(model, X_train)

################################
# Logistic
################################
params1 = EvoTreeRegressor(T=Float32,
    loss=:logistic, metric=:none,
    nrounds=200, nbins = 64,
    λ = 0.5, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = fit_evotree_gpu(params1, X_train, Y_train);
@time pred_train_linear = predict_gpu(model, X_train)

x_perm = sortperm(X_train[:,1])
plot(X_train, Y_train, msize = 1, mcolor = "gray", mswidth=0, background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color = "navy", linewidth = 1.5, label = "Logistic")
# savefig("figures/regression_sinus_gpu.png")

params1 = EvoTreeRegressor(T=Float32,
    loss=:logistic, metric=:logloss,
    nrounds=200, nbins = 64,
    λ = 0.5, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0)

@time model = fit_evotree_gpu(params1, X_train, Y_train, print_every_n = 25);
@time model = fit_evotree_gpu(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n = 25);
@time pred_train_linear = predict_gpu(model, X_train)

################################
# Gaussian
################################
params1 = EvoTreeGaussian(T=Float64,
    loss=:gaussian, metric=:gaussian,
    nrounds=200, nbins=64,
    λ = 1.0, γ=0.1, η=0.1,
    max_depth = 5, min_weight = 100.0,
    rowsample=0.5, colsample=1.0, rng=123)

@time model = fit_evotree_gpu(params1, X_train, Y_train, print_every_n = 25);
@time pred_train_gauss = predict_gpu(model, X_train)

pred_gauss = [Distributions.Normal(pred_train_gauss[i,1], pred_train_gauss[i,2]) for i in 1:size(pred_train_gauss,1)]
pred_q80 = quantile.(pred_gauss, 0.8)
pred_q20 = quantile.(pred_gauss, 0.2)

mean(Y_train .< pred_q80)
mean(Y_train .< pred_q20)

x_perm = sortperm(X_train[:,1])
plot(X_train[:, 1], Y_train, ms = 1, mcolor = "gray", mswidth=0, background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train_gauss[x_perm, 1], color = "navy", linewidth = 1.5, label = "mu")
plot!(X_train[:,1][x_perm], pred_train_gauss[x_perm, 2], color = "darkred", linewidth = 1.5, label = "sigma")
plot!(X_train[:,1][x_perm], pred_q20[x_perm, 1], color = "green", linewidth = 1.5, label = "q20")
plot!(X_train[:,1][x_perm], pred_q80[x_perm, 1], color = "green", linewidth = 1.5, label = "q80")
# savefig("figures/gaussian_sinus_gpu.png")

params1 = EvoTreeGaussian(T=Float32,
    loss=:gaussian, metric=:gaussian,
    nrounds=100, nbins=64,
    λ = 1.0, γ=0.1, η=0.1,
    max_depth = 5, min_weight = 20.0,
    rowsample=0.5, colsample=1.0, rng=123)

@time model = fit_evotree_gpu(params1, X_train, Y_train, print_every_n = 25);
@time model = fit_evotree_gpu(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n = 25);
@time pred_train_linear = predict_gpu(model, X_train)