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.
##########################################################################

mutable struct Location
	x::Float # longitude, or other
	y::Float # latitude, or other

	Location() = new(nullX, nullY)
	Location(x,y) = new(x,y)
end

mutable struct Point
	index::Int
	location::Location
	value::Any

	# node nearest to location
	nearestNodeIndex::Int
	nearestNodeDist::Float

	Point() = new(nullIndex, Location(), nothing,
		nullIndex, nullDist)
end

mutable struct Node
	index::Int
	location::Location
	offRoadAccess::Bool # if node can be used to get on-road and off-road

	fields::Dict{String,Any} # additional data, not used by simulation

	Node() = new(nullIndex, Location(), true,
		Dict{String,Any}())
end

mutable struct Arc
	index::Int
	fromNodeIndex::Int
	toNodeIndex::Int
	distance::Float # NaN if not set

	fields::Dict{String,Any} # additional data, not used by simulation

	Arc() = new(nullIndex, nullIndex, nullIndex, NaN,
		Dict{String,Any}())
end

# graph data for network (actually a digraph)
mutable struct Graph
	# parameters:
	isReduced::Bool
	nodes::Vector{Node}
	arcs::Vector{Arc}

	arcDists::Vector{Float} # arcDists[i] = arcs[i].distance

	light::LightGraphs.DiGraph
	fadjList::Vector{Vector{Int}} # shorthand for light.fadjlist; fadjList[i] gives nodes that are connected to node i by an outgoing arc from node i
	badjList::Vector{Vector{Int}} # shorthand for light.badjList; badjList[i] gives nodes that are connected to node i by an outgoing arc to node i

	# for full graph:
	nodePairArcIndex::SparseMatrixCSC{Int,Int} # for node indices i,j, nodePairArcIndex[i,j] should return arc index for node pair
	# cannot always use nodePairArcIndex for reduced graph, there may be more than one arc from node i to j, use spNodePairArcIndex instead

	Graph(isReduced::Bool) = new(isReduced, [], [],
		[],
		LightGraphs.DiGraph(), [], [],
		spzeros(Int, 0, 0))
end

# travel data for network
mutable struct NetTravel
	# parameters:
	isReduced::Bool
	modeIndex::Int
	arcTimes::Vector{Float} # arcTimes[i] gives travel time along arc i

	# for use in reduced graph:
	arcDists::Vector{Float} # arcDists[i] = distance of arc i; reference to Graph.arcDists, useful for checking that loaded rNetTravel is correct
	fadjList::Vector{Vector{Int}} # = rGraph.fadjList; useful for checking that loaded rNetTravel is correct
	spTimes::Array{FloatSpTime,2} # spTimes[i,j] = shortest path time between node i and j
	spDists::Array{FloatSpDist,2} # spDists[i,j] = shortest path distance between node i and j
	spFadjIndex::Array{IntFadj,2} # for shortest path from rNode i to j, spFadjIndex[i,j] gives the index (in fadjList[i], see Graph) of the successor rNode of i for this path
	spNodePairArcIndex::SparseMatrixCSC{Int,Int} # spNodePairArcIndex[i,j] = index of arc incident to nodes i,j, if it provides shortest travel time
	spFadjArcList::Vector{Vector{Int}} # spFadjArcList[i][j] gives index of rArc from rNode[i] to jth outgoing rNode of rNode[i] (i.e. rGraph.fadjList[i][j])

	# for use in full graph:
	fNodeToRNodeTime::Vector{Dict{Int,Float}} # fNodeToRNodeTime[i][j] gives time from fNode[i] to rNode[j], as long as rNode[j] is in fNodeToRNodes[i] (see type Network)
	fNodeFromRNodeTime::Vector{Dict{Int,Float}} # fNodeFromRNodeTime[i][j] gives time to fNode[i] from rNode[j], as long as rNode[j] is in fNodeFromRNodes[i] (see type Network)
	rArcFNodesTimes::Vector{Vector{Float}} # rArcFNodesTimes[i][k] gives travel time from rArc[i].fromNodeIndex to rArcFNodes[i][k]

	# useful object to/from node data, for full graph (see commonFNodes in Network):
	commonFNodeToFNodeTime::Array{Float,2} # commonFNodeToFNodeTime[i,j] gives shortest path time from commonFNodes[i] to fNode j
	fNodeToCommonFNodeTime::Array{Float,2} # fNodeToCommonFNodeTime[i,j] gives shortest path time from fNode i to commonFNodes[j]
	commonFNodeToFNodeDist::Array{Float,2} # commonFNodeToFNodeDist[i,j] gives shortest path distance from commonFNodes[i] to fNode j
	fNodeToCommonFNodeDist::Array{Float,2} # fNodeToCommonFNodeDist[i,j] gives shortest path distance from fNode i to commonFNodes[j]
	commonFNodeToFNodeRNodes::Array{Tuple{Int,Int},2} # commonFNodeToFNodeRNodes[i,j] gives start and end rNode for shortest path from commonFNodes[i] to fNode j
	fNodeToCommonFNodeRNodes::Array{Tuple{Int,Int},2} # fNodeToCommonFNodeRNodes[i,j] gives start and end rNode for shortest path from fNode i to commonFNodes[j]
	fNodeNearestHospitalIndex::Vector{Int} # fNodeNearestHospitalIndex[i] gives index of nearest hospital from fNode[i]

	NetTravel(isReduced::Bool) = new(isReduced, nullIndex, [],
		[], [], Array{FloatSpTime,2}(undef,0,0), Array{FloatSpDist,2}(undef,0,0), Array{IntFadj,2}(undef,0,0), spzeros(Int, 0, 0), [],
		[], [], [],
		Array{Float,2}(undef,0,0), Array{Float,2}(undef,0,0), Array{Float,2}(undef,0,0), Array{Float,2}(undef,0,0), Array{Tuple{Int,Int},2}(undef,0,0), Array{Tuple{Int,Int},2}(undef,0,0), [])
end

mutable struct Network
	fGraph::Graph # full graph
	rGraph::Graph # reduced graph

	# travel data for each graph:
	fNetTravels::Vector{NetTravel} # fNetTravels[i] gives fGraph travel data for travel mode i
	rNetTravels::Vector{NetTravel} # rNetTravels[i] gives rGraph travel data for travel mode i

	# for converting between full and reduced graph:
	rNodeFNode::Vector{Int} # rNodeFNode[i] returns index of rNode i in fGraph
	fNodeRNode::Vector{Int} # fNodeRNode[i] returns index of fNode i in rGraph, nullIndex if not in rGraph
	rArcFNodes::Vector{Vector{Int}} # rArcFNodes[i] gives array of fNode indices that belong to rArc i (ordered: fromNodeIndex -> toNodeIndex)
	fNodeRArcs::Vector{Vector{Int}} # fNodeRArcs[i] gives indices of rArcs that fNode[i] is on
	rArcFNodeIndex::Vector{Dict{Int,Int}} # rArcFNodeIndex[i][j] gives index that fNode[j] appears in rArc[i]; should be same as findall(rArcFNodes[i] .== j), so rArcFNodes[i][rArcFNodeIndex[i][j]] = j
	fNodeToRNodes::Vector{Vector{Int}} # fNodeToRNodes[i] gives indices of rNodes that fNode[i] can travel to, and is "near" ("near" meaning that the fNode is on an rArc incident to rNode, or fNode is rNode in rGraph; it is similar to adjacency)
	fNodeFromRNodes::Vector{Vector{Int}} # fNodeFromRNodes[i] gives indices of rNodes that can travel to fNode[i], and is "near"
	fNodeToRNodeNextFNode::Vector{Dict{Int,Int}} # fNodeToRNodeNextFNode[i][j] gives index of fNode after fNode[i] on path to rNode[j], where j is in fNodeToRNodes[i]. Needed to find path from an fNode to an rNode
	rArcFArcs::Vector{Vector{Int}} # rArcFArcs[i][j] gives the index of the jth fArc on rArcs[i]

	# distance
	fNodeToRNodeDist::Vector{Dict{Int,Float}} # fNodeToRNodeDist[i][j] gives distance from fNode[i] to rNode[j], as long as rNode[j] is in fNodeToRNodes[i]
	fNodeFromRNodeDist::Vector{Dict{Int,Float}} # fNodeFromRNodeDist[i][j] gives distance to fNode[i] from rNode[j], as long as rNode[j] is in fNodeFromRNodes[i]

	# fNodes that are common start/end points for travel (e.g, fNodes nearest to stations, hospitals):
	commonFNodes::Vector{Int} # list of common fNodes; = findall(isFNodeCommon)
	isFNodeCommon::Vector{Bool} # isFNodeCommon[i] = true if fNode i is common, false otherwise
	fNodeCommonFNodeIndex::Vector{Int} # fNodeCommonFNodeIndex[i] gives index of fNode i in commonFNodes, if isFNodeCommon[i] = true

	Network() = new(Graph(false), Graph(true),
		[], [],
		[], [], [], [], [], [], [], [], [],
		[], [],
		[], [], [])
end

# travel information, for on-road and off-road
mutable struct TravelMode
	index::Int

	offRoadSpeed::Float

	fNetTravel::NetTravel # reference to single fNetTravel in type Network
	rNetTravel::NetTravel # reference to single rNetTravel in type Network

	TravelMode() = new(nullIndex,
		nullTime,
		NetTravel(false), NetTravel(true))
end

# for storing travel modes, travel sets (a set of travel modes apply to each time period),
# and conditions for when to use each travel set/mode
mutable struct Travel
	numModes::Int # number of travel modes
	numSets::Int # number of sets of travel modes; each set may contain multiple travel modes e.g. for different combinations of travel priorities and ambulance classes. Different sets can overlap, by containing the same travel modes.

	modes::Vector{TravelMode}
	modeLookup::Array{Int,2} # modeLookup[i,j] = index of travel mode to use for travel set i, travel priority j. Change this variable according to modelling needs

	setsStartTimes::Vector{Float} # setsStartTimes[i] gives time at which travel set setsTimeOrder[i] should be started (setsStartTimes[i+1] or Inf gives end time)
	setsTimeOrder::Vector{Int} # setsTimeOrder[i] gives travel set to start using at time setsStartTimes[i]
	recentSetsStartTimesIndex::Int # index of most recently used value in setsStartTimes (and setsTimeOrder), should only ever increase in value

	Travel() = new(nullIndex, nullIndex,
		[], Array{Int,2}(undef,0,0),
		[], [], nullIndex)
end

# for storing ambulance routes
# routes include a path (on fGraph), and may start/end off the graph
mutable struct Route
	priority::Priority # travel priority
	travelModeIndex::Int

	# start and end locations and times
	startLoc::Location
	startTime::Float
	endLoc::Location
	endTime::Float

	# start and end nodes and times
	startFNode::Int # index of first fNode in route
	startFNodeTime::Float # time at which startFNode is reached
	startFNodeDist::Float # distance from startLoc to startFNode; can be different from result from normDist() if starting part way along an arc
	startRNode::Int # index of first rNode in route
	startRNodeTime::Float # time at which startRNode is reached
	endRNode::Int # index of last rNode in route
	endRNodeTime::Float # time at which endRNode is reached
	endFNode::Int # index of last fNode in route
	endFNodeTime::Float # time at which endFNode is reached

	firstRArc::Int # index of first rArc in route; = nullIndex if startFNode == endFNode

	## fields that vary throughout the route are below here
	## most of the route code assumes that these fields only change to move "forward" through the route

	status::RouteStatus
	recentUpdateTime::Float # time at which route was most recently updated

	# recent rArc visited ("recent" means it is from the most recent route update/query)
	recentRArc::Int # index of rArc recently visited
	recentRArcStartTime::Float # time that travel on rArc started / would have started if had travelled from first node of recentRArc (can start part way along the arc)
	recentRArcEndTime::Float # time that travel on rArc should end / would end if were travelling to last node of recentRArc (can end part way along the arc)

	# recent fNode visited on recentRArc
	recentRArcRecentFNode::Int # index of fNode recently visited on recentRArc, rArcFNodes[recentRArc][recentRArcRecentFNode] gives index of fNode in fGraph
	recentFNode::Int # = rArcFNodes[recentRArc][recentRArcRecentFNode]
	recentFNodeTime::Float # time that recentFNode was reached

	# next fNode to visit on recentRArc
	recentRArcNextFNode::Int # recentRArcRecentFNode + 1
	nextFNode::Int # = rArcFNodes[recentRArc][recentRArcNextFNode]
	nextFNodeTime::Float # time that nextFNode will be reached
	nextFNodeDist::Float # distance to reach nextFNode; this is only set after calling getRouteNextNodeDist!()

	Route() = new(nullPriority, nullIndex,
		Location(), nullTime, Location(), nullTime,
		nullIndex, nullTime, nullDist, nullIndex, nullTime, nullIndex, nullTime, nullIndex, nullTime,
		nullIndex,
		routeNullStatus, nullTime,
		nullIndex, nullTime, nullTime,
		nullIndex, nullIndex, nullTime,
		nullIndex, nullIndex, nullTime, nullDist)
end

mutable struct Event
	index::Int # index of event in list of events that have occurred (not for sim.eventList)
	parentIndex::Int # index of parent event
	form::EventForm # "type" is taken
	time::Float
	ambIndex::Int
	callIndex::Int
	stationIndex::Int

	Event() = new(nullIndex, nullIndex, nullEvent, nullTime, nullIndex, nullIndex, nullIndex)
end

mutable struct Ambulance
	index::Int
	status::AmbStatus
	stationIndex::Int
	callIndex::Int
	route::Route
	event::Event # next/current event, useful for deleting next event from eventList
	class::AmbClass # class/type of ambulance
	attributes::Dict{String,Any}
	# schedule::Schedule # to be added

	dispatchTime::Float # time at which amb was dispatched
	mobilisationTime::Float # time at which amb will mobilise if there is mobilisation delay; = nullTime if not mobilising

	# for animation:
	currentLoc::Location
	movedLoc::Bool
	destLoc::Location # destination

	# count statistics:
	numCallsTreated::Int # total number of calls that ambulance provided treatment on scene
	numCallsTransported::Int # total number of calls transported to hospital
	numDispatches::Int
	numDispatchesFromStation::Int # total number of dispatches while at station
	numDispatchesWhileMobilising::Int # total number of dispatches while mobilising
	numDispatchesOnRoad::Int # total number of dispatches while on road
	numDispatchesOnFree::Int # total number of dispatches after ambulance becomes free
	numRedispatches::Int # number of times that ambulance is redispatched from one call to another
	numMoveUps::Int
	numMoveUpsFromStation::Int
	numMoveUpsOnRoad::Int
	numMoveUpsOnFree::Int
	numMoveUpsReturnToPrevStation::Int # number of times that ambulance is told to return to the station from which it started move up (see moveUpFromStationIndex)

	# status statistics
	statusDurations::Dict{Union{AmbStatus,AmbStatusSet},Float} # duration spent in each status (or status set); only updated after setting status
	statusDistances::Dict{Union{AmbStatus,AmbStatusSet},Float} # distance travelled in each status (or status set); only updated after setting status
	statusTransitionCounts::Array{Int,2} # statusTransitionCounts[Int(status1), Int(status2)] gives number of times that amb transitioned from status1 to status2 (in AmbStatus enum)

	# for calculating statistics:
	statusSetTime::Float # time at which status was last set, even if set to same status value
	prevStatus::AmbStatus
	moveUpFromStationIndex::Int # index from which ambulance started move up; = nullIndex if not moving up from station

	# dispatch statistics:
	dispatchStartLocCounts::Dict{Location,Int} # dispatchStartLocCounts[loc] gives number of times that ambulance was dispatched from location 'loc'

	# # hospital and station statistics:
	# atHospitalDurations::Vector{Float} # atHospitalDurations[i] gives the duration that the ambulance spent at hospitals[i]
	# atHospitalCounts::Vector{Int} # atHospitalCounts[i] gives the number of times that the ambulance was at hospitals[i]
	# idleAtStationDurations::Vector{Int} # idleAtStationDurations[i] is the duration that the ambulance spent idle at stations[i]

	# # redispatch statistics
	# redispatchCounts::Array{Int,2} # redispatchCounts[Int(p1),Int(p2)] gives number of redispatches from call of priority p1 to call of priority p2

	Ambulance() = new(nullIndex, ambNullStatus, nullIndex, nullIndex, Route(), Event(), nullAmbClass, Dict(),
		nullTime, nullTime,
		Location(), false, Location(),
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		Dict(), Dict(), Array{Int,2}(undef,0,0),
		nullTime, ambNullStatus, nullIndex,
		Dict())
end

mutable struct Call
	index::Int
	priority::Priority
	transport::Bool # true if requires transport to hospital
	hospitalIndex::Int # hospital that call should be transported to. If hospitalIndex == nullIndex, will transport to nearest hospital. See also chosenHospitalIndex
	location::Location # where call occurs
	attributes::Dict{String,Any}

	# time/duration params:
	arrivalTime::Float # time at which call arrives
	dispatchDelay::Float # delay between call arrival and being ready to dispatch an ambulance
	onSceneDuration::Float # time spent at call location
	handoverDuration::Float # for hospital handover

	# node nearest to call location:
	nearestNodeIndex::Int
	nearestNodeDist::Float

	# call state
	status::CallStatus
	ambIndex::Int # index of responding ambulance

	# for animation:
	currentLoc::Location
	movedLoc::Bool

	# statistics:
	dispatchTime::Float # time at which final ambulance was dispatched (= arrivalTime + dispatchDelay, unless call queued/bumped)
	ambArrivalTime::Float # time at which ambulance arrives on-site
	hospitalArrivalTime::Float # time at which ambulance arrives at hospital
	numBumps::Int # total number of times that call gets bumped due to redispatch
	wasQueued::Bool # whether call was queued or not
	ambDispatchLoc::Location # location of responding ambulance at moment of dispatch
	ambStatusBeforeDispatch::AmbStatus # status of ambulance just before dispatch
	chosenHospitalIndex::Int # hospital (if any) that call was transported to

	# duration statistics:
	queuedDuration::Float # duration between being ready to dispatch (but no ambulances are free) and dispatching an ambulance
	bumpedDuration::Float # total duration spent waiting for ambulances that were later redirected away from call due to bumping
	waitingForAmbDuration::Float # total duration spent waiting for ambulances that have been dispatched, even if call bumped and ambulance did not reach call
	responseDuration::Float # duration between call arrival and ambulance arrival at call location
	ambGoingToCallDuration::Float # duration between dispatch of ambulance that responded and ambulance reaching call
	transportDuration::Float # duration for transporting call to hospital
	serviceDuration::Float # duration between call arrival and ambulance finishing with call

	# for calculating statistics:
	statusSetTime::Float # time at which status was last set, even if set to same status value

	Call() = new(nullIndex, nullPriority, true, nullIndex, Location(), Dict(),
		nullTime, nullTime, nullTime, nullTime,
		nullIndex, nullDist,
		callNullStatus, nullIndex,
		Location(), false,
		nullTime, nullTime, nullTime, 0, false, Location(), ambNullStatus, nullIndex,
		0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
		nullTime)
end

mutable struct Hospital
	index::Int
	location::Location
	attributes::Dict{String,Any}

	nearestNodeIndex::Int
	nearestNodeDist::Float

	# for statistics:
	numCalls::Int # total number of calls taken to hospital
	# numCallsTotalHandoverDuration::OffsetArray{Float,1,Vector{Float}} # numCallsTotalHandoverDuration[k] gives total duration that hospital had k calls for handover

	Hospital() = new(nullIndex, Location(), Dict(),
		nullIndex, nullDist,
		0)
end

mutable struct Station
	index::Int
	location::Location
	capacity::Int # maximum number of ambulances that station can hold
	attributes::Dict{String,Any}

	nearestNodeIndex::Int
	nearestNodeDist::Float

	# statistics:
	numIdleAmbsTotalDuration::OffsetArray{Float,1,Vector{Float}} # numIdleAmbsTotalDuration[k] gives total duration that station had k idle ambulances
	currentNumIdleAmbs::Int # current number of idle ambulances at station
	currentNumIdleAmbsSetTime::Float # time at currentNumIdleAmbs was last set (even if set to same value)

	# move up statistics:
	moveUpStartLocCounts::Dict{Location,Int} # locations (and counts) from which ambulances have started move up towards this station

	Station() = new(nullIndex, Location(), 0, Dict(),
		nullIndex, nullDist,
		OffsetVector(Float[],0), 0, nullTime,
		Dict())
end

mutable struct MobilisationDelay
	use::Bool # true if using mobilisation delay

	distrRng::DistrRng
	expectedDuration::Float

	MobilisationDelay() = new(false, DistrRng(Normal(0,0), seed = 0), 0.0)
	MobilisationDelay(use, distrRng, expectedDuration) = new(use, distrRng, expectedDuration)
end

# conditions that determine ambulance redispatch behaviour
mutable struct Redispatch
	allow::Bool # true if redispatch is allowed, false otherwise
	conditions::Array{Bool,2} # conditions[Int(p1),Int(p2)] returns true if ambulance may be redispatched from its current call (priority p1) to another call (priority p2). Change this var to match modelling needs.

	Redispatch() = new(false, fill(false,numPriorities,numPriorities))
	function Redispatch(::Type{Val{:default}})
		redispatch = Redispatch()
		redispatch.allow = true
		@assert(!any(redispatch.conditions)) # all values should be false so far
		redispatch.conditions[Int(lowPriority), Int(highPriority)] = true
		redispatch.conditions[Int(medPriority), Int(highPriority)] = true
		return redispatch
	end
end

mutable struct Map
	xMin::Float
	xMax::Float
	yMin::Float
	yMax::Float
	xRange::Float # xMax - xMin
	yRange::Float # yMax - yMin
	xScale::Float # convert delta(x) to distance
	yScale::Float # convert delta(y) to distance
	# xRangeDist::Float # xRange * xScale
	# yRangeDist::Float # yRange * yScale

	Map() = new(nullX, nullX, nullY, nullY, nullDist, nullDist, nullDist, nullDist)
	Map(xMin, xMax, yMin, yMax, xScale, yScale) = new(xMin, xMax, yMin, yMax, xMax - xMin, yMax - yMin, xScale, yScale)
end

# grid search rectangle, part of type Grid
# for keeping track of search progress,
# while looking for nearest node in grid
mutable struct GridSearchRect
	# distance of search rectangle borders from a given location
	xDist::Vector{Float}
	yDist::Vector{Float}

	# min and max range of grid indices to search next
	ixSearch::Vector{Int}
	iySearch::Vector{Int}

	# min and max grid indices of search rectangle, for each direction
	# contains indices of next search also
	ixSearched::Vector{Int}
	iySearched::Vector{Int}

	GridSearchRect() = new([nullDist, nullDist], [nullDist, nullDist],
		[nullIndex, nullIndex], [nullIndex, nullIndex],
		[nullIndex, nullIndex], [nullIndex, nullIndex])
	GridSearchRect(ix, iy) = new([nullDist, nullDist], [nullDist, nullDist],
		[ix, ix], [iy, iy],
		[ix, ix], [iy, iy])
end

# grid rectangle, part of type Grid
mutable struct GridRect
	nodeIndices::Vector{Int}

	GridRect() = new([])
end

# grid breaks map region into rectangles,
# each rectangle stores graph node indices (from full graph)
# this is for quickly finding the nearest node to a location
mutable struct Grid
	nx::Int # number of divisions in x direction
	ny::Int # number of divisions in y direction
	xRange::Float # size of divisions in x direction, = map.xRange / nx
	yRange::Float # size of divisions in y direction, = map.yRange / ny
	xRangeDist::Float # distance of divisions in x direction, = xRange * map.xScale
	yRangeDist::Float # distance of divisions in y direction, = yRange * map.yScale
	rects::Array{GridRect,2} # rectangles
	searchRect::GridSearchRect

	Grid() = new(nullIndex, nullIndex, nullDist, nullDist, nullDist, nullDist, Array{GridRect,2}(undef, 0, 0), GridSearchRect())
	function Grid(map::Map, nx, ny)
		grid = new(nx, ny, nullDist, nullDist, nullDist, nullDist, Array{GridRect,2}(undef, nx, ny))
		for i = 1:nx
			for j = 1:ny
				grid.rects[i,j] = GridRect()
			end
		end
		grid.xRange = map.xRange / nx
		grid.yRange = map.yRange / ny
		grid.xRangeDist = grid.xRange * map.xScale
		grid.yRangeDist = grid.yRange * map.yScale
		grid.searchRect = GridSearchRect()
		return grid
	end
end

# simple raster type
# stores data for rectangular grid of cells,
# cell (i,j) has centre x[i], y[j], and value z[i,j]
# requires at least 2 cells in each direction
mutable struct Raster
	# parameters:
	x::Vector{Float}
	y::Vector{Float}
	z::Array{Float,2} # z[i,j] corresponds with x[i], y[j]

	nx::Int # length(x)
	ny::Int # length(y)
	dx::Float # x step size
	dy::Float # y step size

	Raster() = new([], [], Array{Float,2}(undef,0,0),
		0, 0, 0.0, 0.0)
	function Raster(x, y, z)
		nx = length(x)
		ny = length(y)
		@assert(nx > 1)
		@assert(ny > 1)
		@assert((nx, ny) == size(z))
		dx = (maximum(x) - minimum(x))/(nx-1)
		dy = (maximum(y) - minimum(y))/(ny-1)
		@assert(dx > 0)
		@assert(dy > 0)
		return new(x, y, z, nx, ny, dx, dy)
	end
end

# for generating random locations within a raster,
# where raster z values are proportional to probability of location being
# within corresponding cell (i.e. z values are for categorical distribution),
# and locations are generated uniformly within the cell
mutable struct RasterSampler
	raster::Raster
	cellDistrRng::DistrRng # to generate index of raster cell
	cellLocRng::AbstractRNG # rng used when generating location within raster cell

	function RasterSampler(raster::Raster, cellRng::AbstractRNG, cellLocRng::AbstractRNG)
		@assert(all(raster.z .>= 0)) # otherwise probability values will be negative
		p = raster.z[:] / sum(raster.z) # pdf for z
		cellSampler = sampler(Categorical(p))
		cellDistrRng = DistrRng(cellSampler, cellRng)
		return new(raster, cellDistrRng, cellLocRng)
	end
	function RasterSampler(raster::Raster, cellSeed::Int, cellLocSeed::Int)
		cellRng = (cellSeed >= 0 ? MersenneTwister(cellSeed) : MersenneTwister(rand(UInt32)))
		cellLocRng = (cellLocSeed >= 0 ? MersenneTwister(cellLocSeed) : MersenneTwister(rand(UInt32)))
		return RasterSampler(raster, cellRng, cellLocRng)
	end
end

# to model call demand for a single priority
# can use multiple demand modes to model change in demand with time, using Demand type
mutable struct DemandMode
	# parameters:
	index::Int # mode index
	rasterIndex::Int # index of demand raster in Demand.rasters
	priority::Priority # demand priority
	arrivalRate::Float # demand per day

	raster::Raster # reference to Demand.rasters[rasterIndex]
	rasterMultiplier::Float # = arrivalRate / sum(raster.z), to scale raster.z values to match arrivalRate

	DemandMode() = new(nullIndex, nullIndex, nullPriority, nullTime,
		Raster(), 0.0)
end

# To model demand (calls).
# Stores demand modes, demand sets (a set of demand modes apply to each time period),
# and conditions for when to use each demand set/mode.
mutable struct Demand
	initialised::Bool
	numRasters::Int # number of demand rasters
	numModes::Int # number of demand modes
	numSets::Int # number of sets of demand modes; each set may contain multiple demand modes e.g. for different combinations of call priorities and ambulance classes. Different sets can overlap, by containing the same demand modes.

	rasters::Vector{Raster} # demand rasters, to model spatial demand. Raster cell z values are not actual demand, but are proportional to demand
	rasterFilenames::Vector{String} # rasterFilenames[i] gives filename of rasters[i]

	modes::Vector{DemandMode}
	modeLookup::Array{Int,2} # modeLookup[i,j] = index of demand mode to use for demand set i, demand priority j. Change this variable according to modelling needs

	setsStartTimes::Vector{Float} # setsStartTimes[i] gives time at which demand set setsTimeOrder[i] should be started (setsStartTimes[i+1] or Inf gives end time)
	setsTimeOrder::Vector{Int} # setsTimeOrder[i] gives demand set to start using at time setsStartTimes[i]
	recentSetsStartTimesIndex::Int # index of most recently used value in setsStartTimes (and setsTimeOrder), should only ever increase in value

	Demand() = new(false, nullIndex, nullIndex, nullIndex,
		[], [],
		[], Array{Int,2}(undef,0,0),
		[], [], nullIndex)
end

# For a given set of points, coverage time, and travel mode,
# for each station store the points that can be reached by
# travelling from the station within the coverage time.
mutable struct PointsCoverageMode
	# parameters:
	index::Int
	points::Vector{Point} # reference to DemandCoverage.points
	coverTime::Float # any demand that can be reached within this time gets additional coverage of 1, otherwise 0. Usually equal to target cover time minus dispatch delay.
	travelMode::TravelMode # reference to a travel mode

	pointSets::Vector{Vector{Int}} # pointSets[i] has set of all point indices covered by the same unique set of stations
	stationSets::Vector{Vector{Int}} # stationSets[i] = unique set of stations for pointSets[i]
	stationsCoverPointSets::Vector{Vector{Int}} # stationsCoverPointSets[i] = indices of pointSets covered by station i

	PointsCoverageMode() = new(nullIndex, [], nullTime, TravelMode(),
		[], [], [])
end

mutable struct DemandCoverage
	# params
	coverTimes::Dict{Priority,Float} # coverTimes[p] gives the target cover time for demand of priority p
	rasterCellNumRows::Int # number of rows of points to create per demand raster cell
	rasterCellNumCols::Int # number of columns of points to create per demand raster cell

	initialised::Bool
	points::Vector{Point} # demand is aggregated to points, same points are used for all demand rasters
	nodesPoints::Vector{Vector{Int}} # nodesPoints[i] gives indices of points for which node i is the nearest node
	rastersPointDemands::Vector{Vector{Float}} # rastersPointDemands[i][j] is demand at points[j] for Demand.rasters[i]

	pointsCoverageModes::Vector{PointsCoverageMode}
	pointsCoverageModeLookup::Vector{Dict{Float,Int}} # pointsCoverageModeLookup[TravelMode.index][coverTime] gives index of PointsCoverageMode
	pointSetsDemands::Array{Vector{Float},2} # pointSetsDemands[PointsCoverageMode.index, DemandMode.rasterIndex] gives relative demand values for each point set in PointsCoverageMode.pointSets, for Demand.rasters[rasterIndex]. Note that this needs to be multiplied by DemandMode.rasterMultiplier to get absolute (instead of relative) demand values.

	DemandCoverage() = new(Dict(), 0, 0,
		false, [], [], [],
		[], [], Array{Vector{Float},2}(undef,0,0))

	function DemandCoverage(coverTimes::Dict{Priority,Float}, rasterCellNumRows::Int, rasterCellNumCols::Int)
		dc = demandCoverage = DemandCoverage()
		dc.coverTimes = coverTimes
		dc.rasterCellNumRows = rasterCellNumRows
		dc.rasterCellNumCols = rasterCellNumCols
		return demandCoverage
	end

	function DemandCoverage(demandCoverage::DemandCoverage)
		dc = demandCoverage # shorthand
		return DemandCoverage(dc.coverTimes, dc.rasterCellNumRows, dc.rasterCellNumCols)
	end
end

# Mean and half the width of the confidence interval of the mean.
# Half-width is also known as half-length, or margin of error.
struct MeanAndHalfWidth
	mean::Float
	halfWidth::Float
	MeanAndHalfWidth(m, hw) = new(convert(Float,m), convert(Float,hw))
end

# cover bound simulation replication
mutable struct CoverBoundSimRep
	numCalls::Int
	numFreeAmbsCount::Vector{Int} # numFreeAmbsCount[i] = number of times that calls arrived while there were i free ambs
	queuedCallDurations::Vector{Float} # for calls that arrived when there were no free ambs, record how long each call was queued
	bound::Float # upper bound on fraction of in time responses

	CoverBoundSimRep() = new(0, [], [], 0)
end

# cover bound simulation
mutable struct CoverBoundSim
	numReps::Int
	numAmbs::Int

	warmUpDuration::Float
	minLastCallArrivalTime::Float

	reps::Vector{CoverBoundSimRep}
	bound::MeanAndHalfWidth # upper bound on fraction of in time responses

	interarrivalTimeDistrRng::DistrRng
	ambBusyDurationRng::AbstractRNG

	CoverBoundSim() = new(0, 0,
		nullTime, nullTime,
		[], MeanAndHalfWidth(0,0))
end

# cover bound data for a single combination of conditions (e.g. single demand mode and single travel mode)
mutable struct CoverBoundMode
	# parameters
	index::Int
	responseTravelMode::TravelMode # travel mode for ambulance response
	# demandPriorities::Vector{Priority} # all demand priorities that use this cover bound mode

	# travel times
	stationsToNodesTimes::Array{Float,2} # stationsToNodesTimes[i,j] = travel time from station i to node j
	stationsToPointsTimes::Array{Float,2} # stationsToPointsTimes[i,j] = travel time from station i to demand point j
	pointsToHospitalTimes::Vector{Float} # pointsToHospitalTimes[i] = travel time from point i to nearest hospital

	# group demand into sets
	# use in conjunction with demandCoverage, e.g. demandCoverage.nodesPoints
	nodeSets::Vector{Vector{Int}} # nodeSets[i] gives node indices that have same order in which ambulances can reach from each station
	nodeSetsStationList::Vector{Vector{Int}} # nodeSetsStationList[i] has list of stations ordered (from closest in time to furthest) from nodeSets[i]

	# distributions related to service
	distrs::Dict{String,Distribution} # distrs[name] gives distribution; name from Set("onSceneDuration", "handoverDuration")
	approxDistrs::Dict{String,Distribution} # approxDistrs[name] gives lower bound approximation (discrete, non-parametric distribution) of distrs[name]
	# cdfs::Dict{String,Vector{Float}} # to add? for "d1" and "d2"
	transportProb::Float # probability that call requires transport to hospital; assume call is taken to nearest hospital

	CoverBoundMode() = new(nullIndex, TravelMode(),
		Array{Float,2}(undef,0,0), Array{Float,2}(undef,0,0), [],
		[], [],
		Dict(), Dict(), 1.0)
end

# cover bound - a bound on optimal move-up performance
mutable struct CoverBound
	modes::Vector{CoverBoundMode}
	modeLookup::Vector{Int} # modeLookup[i] gives index of mode (in modes) to use for travel mode i. Will assume for now that there is only one demand set (demand is static) and one travel set (i.e. travel times are static).

	# parameters
	ambBusyDurationsToSample::Vector{Float} # will find upper bound on probability of amb being busy for each duration in this list
	accountForQueuedDurations::Bool # for queued calls, calculate maximum coverage according to time remaining to reach call on time; set to false for original cover bound which calculates coverage of queued calls as if the call was queued for no time
	queuedDurationsToSample::Vector{Float} # will solve MCLP for different coverage times based on how long the calls were queued; use if accountForQueuedDurations == true

	dispatchDelay::Float # assumed constant for simpler modelling; add any mobilisation delay to this

	# amb busy duration
	ambBusyDurationProbUpperBounds::Array{Float,2} # ambBusyDurationProbUpperBounds[i,j] = upper bound on probability of dispatched amb being busy for duration <= ambBusyDurationsToSample[i] for j free ambulances
	ambBusyDurationLowerBoundDistrs::Vector{Sampleable} # ambBusyDurationLowerBoundDistrs[i] is sampleable distribution of amb busy duration for i free ambs

	# coverage
	numAmbsMaxCoverageFrac::Vector{Float} # numAmbsMaxCoverageFrac[i] gives maximum fraction of demand that can be covered with i free ambulances (from 1:numStations), from solving MCLP
	queuedDurationsMaxCoverageFrac::Vector{Float} # queuedDurationsMaxCoverageFrac[i] = maximum coverage for call that was queued for a duration of queuedDurationsToSample[i]; use if accountForQueuedDurations == true

	sim::CoverBoundSim

	CoverBound() = new([], [],
		[], false, [],
		0.0,
		Array{Float,2}(undef,0,0), [],
		[], [],
		CoverBoundSim())
end

# move up data types
abstract type MoveUpDataType end
mutable struct EmptyMoveUpData <: MoveUpDataType end

# compliance table data
mutable struct CompTableData <: MoveUpDataType
	# parameters:
	compTable::CompTable # compTable[i,j] = number of ambulances to place at station j, with i total free ambs

	compTableStationSlots::Vector{Vector{Int}} # sum(compTableStationSlots[i] .== j) == compTable[i,j]

	# arrays for recycling:
	ambMovable::Vector{Bool} # ambMovable[i] = true if ambulance i can be moved up, false otherwise

	CompTableData() = new(CompTable(undef,0,0),
		[],
		[])
end

# ddsm - dynamic double standard model
mutable struct DdsmData <: MoveUpDataType
	# parameters:
	coverFractionTargetT1::Float # fraction of demand to cover within coverTimes[1]
	travelTimeCost::Float
	slackWeight::Float
	coverTimeDemandPriorities::Vector{Priority} # coverTimeDemandPriorities[i] is demand priority for coverTimes[i]

	coverTimes::Vector{Float} # target coverage response durations; coverTimes[1] <= coverTimes[2]
	options::Dict{Symbol,Any}

	DdsmData() = new(0.0, 0.0, 0.0, [],
		[], Dict())
end

# dmexclp - dynamic maximum expected coverage location problem
mutable struct DmexclpData <: MoveUpDataType
	# parameters:
	busyFraction::Float # fraction for which each ambulance is busy, approximate
	demandWeights::Dict{Priority,Float} # weight of each demand priority on the objective function
	# some other relevant parameters are stored in sim: demand, demandCoverage, responseTravelPriorities

	marginalBenefit::Vector{Float} # marginalBenefit[i] = benefit of adding an ith ambulance to cover single demand, calculated from busyFraction

	# arrays for recycling:
	stationNumFreeAmbs::Vector{Int} # stationNumFreeAmbs[i] = number of free ambulances assigned to station i
	stationMarginalCoverages::Vector{Float} # stationMarginalCoverages[i] gives extra coverage provided from placing newly freed ambulance at station i
	# pointSetsCoverCounts::Vector{Vector{Int}} # pointSetsCoverCounts[i][j] = number of free ambulances covering node set j, for demand.pointsCoverageModes i

	DmexclpData() = new(0.0, Dict(),
		[],
		[], [])
end

# multi-compliance table data
mutable struct MultiCompTableData <: MoveUpDataType
	# parameters:
	multiCompTable::MultiCompTable # multiCompTable[i][j][k] is jth option for number of ambs to assign to station k with i free ambs

	MultiCompTableData() = new(MultiCompTable())
end

mutable struct PriorityListData <: MoveUpDataType
	# parameters:
	priorityList::Vector{Int} # priorityList[i] gives station index that the ith free ambulance should be moved to

	# arrays for recycling:
	stationNumFreeAmbs::Vector{Int} # number of free ambulances assigned to each station

	PriorityListData() = new([],
		[])
end

mutable struct ZhangIpData <: MoveUpDataType
	# parameters:
	marginalBenefits::Vector{Vector{Float}} # marginalBenefits[i][j] = benefit of adding a jth ambulance to station i
	stationCapacities::Vector{Int} # stationCapacities[i] is the number of ambulances that station i can hold on completion of move up, should be <= stations[i].capacity
	travelTimeCost::Float # travel time cost multiplier
	onRoadMoveUpDiscountFactor::Float # discount of travel cost of move up of ambulance on-road and with "regret" travel time <= regretTravelTimeThreshold
	regretTravelTimeThreshold::Float
	expectedHospitalHandoverDuration::Float

	stationSlots::Vector{Int}
	benefitSlots::Vector{Float}
	marginalBenefitsDecreasing::Bool # true if benefit of adding ambulance j to a station is < benefit of adding ambulance j-1, for all stations
	stationSlotsOrderPairs::Array{Int,2} # stationSlotsOrderPairs[i,1:2] gives two stationSlots indices, first should be filled (with ambulance) before second

	ZhangIpData() = new([], [], 1.0, 1.0, nullTime, nullTime,
		[], [], false, Array{Int,2}(undef,0,0))
end

mutable struct Temp0Data <: MoveUpDataType
	# parameters:
	busyFraction::Float # ambulance busy fraction - should remove this, make marginalBenefit a parameter
	travelTimeCost::Float # travel time cost multiplier
	maxIdleAmbTravelTime::Float # max travel time for idle ambulances. 0.021 (days) is about 30 minutes
	maxNumNearestStations::Int # number of nearest stations to consider for each ambulance (may include ambulance's station)

	marginalBenefit::Vector{Float} # marginalBenefit[i] = benefit of adding an ith ambulance to a station

	Temp0Data() = new(0.0, 0.0, Inf, 0,
		[])
end

mutable struct Temp1Data <: MoveUpDataType
	# parameters:
	benefit::Array{Vector{Float},2} # benefit[i1,i2][k] = benefit from having k ambulances collectively at stations i1 and i2
	busyFraction::Float # ambulance busy fraction - should remove this, make marginalBenefit a parameter
	stationPairs::Vector{Vector{Int}} # stationPairs[i] gives ith station pair, which is a pair of station indices
	travelTimeCost::Float # travel time cost multiplier
	maxIdleAmbTravelTime::Float # max travel time for idle ambulances
	maxNumNearestStations::Int # number of nearest stations to consider for each ambulance (may include ambulance's station)

	marginalBenefit::Array{Vector{Float},2} # marginal benefit values from ambulance to station allocations, calculated from 'benefit'

	Temp1Data() = new(Array{Vector{Float},2}(undef,0,0), 0.0, Vector{Vector{Int}}(), 0.0, Inf, 0,
		Array{Vector{Float},2}(undef,0,0))
end

mutable struct Temp2Data <: MoveUpDataType
	# parameters:
	benefit::Array{Array{Float,2},2} # benefit[i1,i2][k1,k2] = benefit from having k1 and k2 ambulances at stations i1 and i2, respectively
	busyFraction::Float # ambulance busy fraction - should remove this, make marginalBenefit a parameter
	stationPairs::Vector{Vector{Int}} # stationPairs[i] gives ith station pair, which is a pair of station indices
	travelTimeCost::Float # travel time cost multiplier
	maxIdleAmbTravelTime::Float # max travel time for idle ambulances
	maxNumNearestStations::Int # number of nearest stations to consider for each ambulance (may include ambulance's station)

	marginalBenefit::Array{Array{Float,2},2} # marginal benefit values from ambulance to station allocations, calculated from 'benefit'

	Temp2Data() = new(Array{Array{Float,2},2}(undef,0,0), 0.0, Vector{Vector{Int}}(), 0.0, Inf, 0,
		Array{Array{Float,2},2}(undef,0,0))
end

mutable struct MoveUpData
	# parameters:
	useMoveUp::Bool
	moveUpModule::MoveUpModule # indicates move up module to be used

	# move up data types:
	compTableData::CompTableData
	ddsmData::DdsmData
	dmexclpData::DmexclpData
	multiCompTableData::MultiCompTableData
	priorityListData::PriorityListData
	zhangIpData::ZhangIpData
	temp0Data::Temp0Data
	temp1Data::Temp1Data
	temp2Data::Temp2Data

	MoveUpData() = new(false, nullMoveUpModule,
		CompTableData(), DdsmData(), DmexclpData(), MultiCompTableData(), PriorityListData(), ZhangIpData(), Temp0Data(), Temp1Data(), Temp2Data())
end

mutable struct File
	name::String
	path::String # includes name
	iostream::IOStream
	checksum::UInt

	# File() = new()
	File() = new("", "", IOStream(""), 0)
end

mutable struct EventsFile
	io::IOStream # write simulation trace of events to this file
	eventFilter::Dict{EventForm,Bool} # eventFilter[eventForm] = true if saving events of eventForm to file, false otherwise

	EventsFile() = new(IOStream(""), Dict([e => true for e in instances(EventForm)]))
end

mutable struct Resimulation
	# parameters:
	use::Bool # true if resimulating (will follow event trace), false otherwise
	timeTolerance::Float

	events::Vector{Event} # note that events[i].index may not necessarily equal i (events[i].index >= i)
	eventsChildren::Vector{Vector{Event}} # eventsChildren[i] gives events that are children of sim event i; length(eventsChildren) may be less than number of events if events file does not include all events
	prevEventIndex::Int # index of previous event in events field

	eventFilter::Dict{EventForm,Bool} # eventFilter[eventForm] = true if resimulating with events of eventForm
	doDispatch::Bool # true if can resimulate dispatch decisions, i.e. eventFilter[ambDispatched] == true
	doMoveUp::Bool # true if can resimulate move-up decisions, i.e. eventFilter[ambMoveUpToStation] == true

	Resimulation() = new(false, 0.0,
		[], [], nullIndex,
		Dict(), false, false)
end

# statistics for a single ambulance, or multiple ambulances
mutable struct AmbulanceStats
	ambIndex::Int # for single ambulance

	numCallsTreated::Int # total number of calls that ambulance provided treatment on scene
	numCallsTransported::Int # total number of calls transported to hospital

	numDispatches::Int # all dispatches
	numDispatchesFromStation::Int # total number of dispatches while at station
	numDispatchesWhileMobilising::Int # total number of dispatches while mobilising
	numDispatchesOnRoad::Int # total number of dispatches while on road
	numDispatchesOnFree::Int # total number of dispatches after ambulance becomes free
	numRedispatches::Int # number of times that ambulance is redispatched from one call to another

	numMoveUps::Int
	numMoveUpsFromStation::Int
	numMoveUpsOnRoad::Int
	numMoveUpsOnFree::Int
	numMoveUpsReturnToPrevStation::Int

	statusDurations::Dict{Union{AmbStatus,AmbStatusSet},Float} # duration spent in each status (or status set)
	statusDistances::Dict{Union{AmbStatus,AmbStatusSet},Float} # distance travelled in each status (or status set)
	statusTransitionCounts::Array{Int,2} # statusTransitionCounts[Int(status1), Int(status2)] gives number of times that amb transitioned from status1 to status2

	AmbulanceStats() = new(nullIndex,
		0, 0,
		0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0,
		Dict(), Dict(), Array{Int,2}(undef,0,0))
end

# statistics for a single call, or multiple calls
mutable struct CallStats
	callIndex::Int # for single call
	numCalls::Int # for multiple calls

	numQueued::Int
	numBumped::Int # number of calls bumped once or more
	numBumps::Int # number of bumps from all calls
	numTransports::Int
	numResponsesInTime::Int # number of in time ambulance responses

	# durations
	totalDispatchDelay::Float # delay between call arrival and being ready to dispatch an ambulance
	totalOnSceneDuration::Float # duration spent at call location
	totalHandoverDuration::Float # for hospital handover
	totalQueuedDuration::Float # duration between being ready to dispatch (but no ambulances are free) and dispatching an ambulance
	totalBumpedDuration::Float # total duration spent waiting for ambulances that were later redirected away from call due to bumping
	totalWaitingForAmbDuration::Float # total duration spent waiting for ambulances that have been dispatched, even if call bumped and ambulance did not reach call
	totalResponseDuration::Float # duration between call arrival and ambulance arrival at call location
	totalAmbGoingToCallDuration::Float # duration of ambulance that responded to reach call
	totalTransportDuration::Float # duration for transporting call to hospital
	totalServiceDuration::Float # duration between call arrival and ambulance finishing with call

	# totalQueueLengthDuration::OffsetArray{Float,1,Vector{Float}} # totalQueueLengthDuration[k] = total duration spent with k queued calls
	responseDurationHist::Histogram

	CallStats() = new(nullIndex, 0,
		0, 0, 0, 0, 0,
		0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
		NullHist())
end

# statistics for a single hospital, or multiple hospitals
mutable struct HospitalStats
	hospitalIndex::Int # for single hospital

	numCalls::Int # total number of calls taken to hospital

	HospitalStats() = new(nullIndex,
		0)
end

# statistics for single station, or multiple stations
mutable struct StationStats
	stationIndex::Int # for single station

	numIdleAmbsTotalDuration::OffsetArray{Float,1,Vector{Float}} # numIdleAmbsTotalDuration[k] gives total duration that station had k idle ambulances

	# numMoveUpsOnFree::Int # number of move-ups to this station of ambulances that have just become free
	# numMoveUpsOnRoad::Int # number of move-ups to this station that were from ambulances that were driving on the road
	# numMoveUpsFromHospitals::Vector{Int} # numMoveUpsFromHospitals[i] is the number of move-ups to this station of ambulances from hospitals[i]
	# numMoveUpsFromStations::Vector{Int} # numMoveUpsFromStations[i] is the number of move-ups of ambulances from stations[i] to this station
	# numMoveUpsToStations::Vector{Int} # numMoveUpsToStations[i] is the number of move-ups of ambulances from this station to stations[i]
	# also count how many move-ups were actually completed, and how many were cut short?

	StationStats() = new(nullIndex,
		OffsetVector(Float[],0))
		# 0, 0, [], [], [])
end

# simulation statistics for a period
mutable struct SimPeriodStats
	# statistics are for [startTime, endTime)
	startTime::Float # start time for stats
	endTime::Float # end time for stats
	duration::Float # = endTime - startTime

	ambulances::Vector{AmbulanceStats}
	hospitals::Vector{HospitalStats}
	stations::Vector{StationStats}

	# aggregate values
	ambulance::AmbulanceStats # for all ambulances
	call::CallStats # statistics of calls that arrived within [startTime, endTime)
	callPriorities::Dict{Priority,CallStats} # callPriorities[p] is statistics of priority p calls that arrived within [startTime, endTime)
	hospital::HospitalStats # for all hospitals
	station::StationStats # for all stations

	# numAmbsFreeTotalDuration::OffsetArray{Float,1,Vector{Float}} # numAmbsFreeTotalDuration[k] gives total duration that k ambulances were free (working but not busy)

	SimPeriodStats() = new(nullTime, nullTime, nullTime,
		[], [], [],
		AmbulanceStats(), CallStats(), Dict(), HospitalStats(), StationStats())
end

mutable struct SimStats
	captures::Vector{SimPeriodStats} # from sim start to some time
	periods::Vector{SimPeriodStats} # arbitrary start and end times; calculated from captures, periods[i] = captures[i] - captures[i-1]

	# params
	doCapture::Bool # true if capturing statistics during simulation
	warmUpDuration::Float # first capture will be at sim.startTime + warmUpDuration (unless warmUpDuration = 0)
	periodDurationsIter::Stateful # iterator that returns duration of next stats period
	nextCaptureTime::Float # time at which statistics will next be captured; updated according to periodDurationsIter

	# sim timestamps
	simStartTime::Float
	simEndTime::Float
	warmUpEndTime::Float # = simStartTime + warmUpDuration
	lastCallArrivalTime::Float # useful for ignoring cool-down period

	# misc
	recordDispatchStartLocCounts::Bool # if true, will save starting location of ambulance dispatches in ambulance.dispatchStartLocCounts
	recordMoveUpStartLocCounts::Bool # if true, will save starting location of ambulance move up in station.moveUpStartLocCounts

	# response duration histogram
	recordResponseDurationHist::Bool # if true, will save call response durations to histogram stored in callStats.responseDurationHist
	responseDurationHistBinWidth::Float # width (duration) of bins for callStats.responseDurationHist
	responseDurationHistNumBins::Int # number of bins for callStats.responseDurationHist

	SimStats() = new([], [],
		true, 0.0, Stateful([]), Inf,
		nullTime, nullTime, nullTime, nullTime,
		false, false,
		false, nullTime, 0)
end

mutable struct Simulation
	startTime::Float
	time::Float # time of most recent event, or time of recent animation frame if animating
	endTime::Float # calculated after simulating

	# world:
	net::Network
	travel::Travel
	map::Map
	grid::Grid

	ambulances::Vector{Ambulance}
	calls::Vector{Call}
	hospitals::Vector{Hospital}
	stations::Vector{Station}

	# shorthand:
	numAmbs::Int # length(ambulances)
	numCalls::Int # length(calls)
	numHospitals::Int # length(hospitals)
	numStations::Int # length(stations)

	eventList::Vector{Event} # events to occur now or in future
	eventIndex::Int # index of event in events that have occurred
	queuedCallList::Vector{Call} # keep track of queued calls. Calls can be queued after call arrivalTime + dispatchDelay

	resim::Resimulation

	mobilisationDelay::MobilisationDelay

	# decision logic
	addCallToQueue!::Function
	findAmbToDispatch!::Function
	redispatch::Redispatch
	moveUpData::MoveUpData

	# demand
	demand::Demand
	demandCoverage::DemandCoverage

	responseTravelPriorities::Dict{Priority,Priority} # responseTravelPriorities[p] gives the travel priority for responding to call of priority p
	targetResponseDurations::Vector{Float} # targetResponseDurations[Int(priority)] gives maximum desired response duration for call of given priority

	# for animation:
	currentCalls::Set{Call} # all calls between arrival and service finish at current time
	previousCalls::Set{Call} # calls in currentCalls for previous frame

	stats::SimStats

	reps::Vector{Simulation} # simulation replications, should have calls and stats defined
	numReps::Int # length(reps)
	isRunnable::Bool # whether or not simulation object can be run / simulated

	# files/folders:
	inputPath::String
	outputPath::String
	inputFiles::Dict{String,File} # given input file name (e.g. "ambulances"), returns file information
	outputFiles::Dict{String,File}

	eventsFile::EventsFile

	writeOutput::Bool # true if writing output files during simulation (e.g. events file), false otherwise (e.g. animating)
	initialised::Bool # true if simulation has been initialised and so can be run, false otherwise
	used::Bool # true if simulation has started running (and so restarting would require copying from backup)
	complete::Bool # true if simulation has ended (no events remaining)
	animating::Bool # true if being used for animation, false otherwise

	backup::Simulation # copy of simulation, for restarts (does not include a backup of all fields in order to save on memory, see backup! function for missing fields)

	configRootElt::XMLElement

	Simulation() = new(nullTime, nullTime, nullTime,
		Network(), Travel(), Map(), Grid(),
		[], [], [], [],
		0, 0, 0, 0,
		[], 0, [],
		Resimulation(),
		MobilisationDelay(),
		nullFunction, nullFunction, Redispatch(Val{:default}), MoveUpData(),
		Demand(), DemandCoverage(),
		Dict(), [],
		Set(), Set(),
		SimStats(),
		[], 0, true,
		"", "", Dict(), Dict(),
		EventsFile(),
		false, false, false, false, false)
end


text:
using GnuplotRecipes
using Documenter

DocMeta.setdocmeta!(GnuplotRecipes, :DocTestSetup, :(using GnuplotRecipes, Gnuplot); recursive=true)

makedocs(;
    modules=[GnuplotRecipes],
    authors="Michal Sojka <michal.sojka@cvut.cz> and contributors",
    repo="https://github.com/wentasah/GnuplotRecipes.jl/blob/{commit}{path}#{line}",
    sitename="GnuplotRecipes.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://wentasah.github.io/GnuplotRecipes.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/wentasah/GnuplotRecipes.jl",
    devbranch="main",
)

ext: jl


text:
"""
    eaCV( channel :: Choi, method :: Symbol = :primal )

Numerically evaluates the entanglement-assisted communication value for the
quantum channel ``\\mathcal{N}`` given in the [`Choi`](@ref) operator representation.
This quantity is defined as:

```math
\\begin{aligned}
    \\text{eaCV}(\\mathcal{N}) & = \\max_{\\Omega^{AB}} \\text{Tr}[\\Omega^{AB}J_{\\mathcal{N}}] \\\\
    \\hspace{1cm}
    & \\text{s.t.} \\quad  \\text{Tr}_A[\\Omega^{AB}] = \\mathbb{I}^B; \\\\
    \\hspace{1cm}
    & \\qquad \\Omega^{AB} \\in \\text{POS}(A\\;:\\;B).
\\end{aligned}
```

Note that for a channel ``\\mathcal{N}`` with output Hilbert space dimension ``d_B``,
the entanglement-assisted CV is bound as

```math
\\text{eaCV}(\\mathcal{N}) \\leq d_B \\cdot \\text{cv}(\\mathcal{N}).
```

The `method` parameter accepts the values:
* `:primal` - Solves the primal problem (see [`eaCVPrimal`](@ref)).
* `:dual` - Solves the dual problem (see [`eaCVDual`](@ref)).

# Returns

A Tuple containing `(max_cv, optimizer)` where `max_cv` is the communication value
and `optimizer` is the matrix which achieves the optimum for the particular problem.
"""
eaCV(channel :: Choi, method::Symbol=:primal) = eaCV(channel, Val(method))
eaCV(channel :: Choi, ::Val{:primal}) = eaCVPrimal(channel.JN, channel.in_dim, channel.out_dim)
eaCV(channel :: Choi, ::Val{:dual}) = eaCVDual(channel.JN, channel.in_dim, channel.out_dim)

"""
    eaCVPrimal(
        ρ :: Matrix{<:Number},
        dimA :: Int,
        dimB :: Int
    ) :: Tuple{Float64, Matrix{ComplexF64}}

This function solves the SDP
```math
\\max \\{ \\langle \\rho, X \\rangle :  \\text{Tr}_{A}(X) = I_{B} , X \\succeq 0 \\}
```
and returns the optimal value and the optimizer, X.
This is the primal problem of the SDP corresponding to the entanglement-assisted commmunication value.
It is related to the channel min-entropy ``H_{\\min}(A|B)_{\\mathcal{J}(\\mathcal{N})}`` by
``\\text{eaCV}(\\mathcal{N}) = 2^{-H_{\\min}(A|B)_{\\mathcal{J}(\\mathcal{N})}}``.
(See [Section 6.1 of this reference](https://arxiv.org/abs/1504.00233) for further details about
the min-entropy.)
Note: we label the primal as the maximization problem.
"""
function eaCVPrimal(ρ :: AbstractArray, dimA :: Int, dimB :: Int) :: Tuple{Float64,  Matrix{ComplexF64}}
    X = HermitianSemidefinite(dimA*dimB)
    objective = real(tr(ρ' * X))
    constraint = partialtrace(X, 1, [dimA,dimB]) == Matrix{Float64}(I,dimB,dimB)
    problem = maximize(objective,constraint)
    qsolve!(problem)
    return problem.optval, X.value
end

"""
    eaCVDual(
        ρ :: Matrix{<:Number},
        dimA :: Int,
        dimB :: Int
    ) :: Tuple{Float64, Matrix{ComplexF64}}

This function solves the SDP
```math
\\min \\{ \\text{Tr}(Y) :  I_{A} \\otimes Y \\succeq \\rho, Y \\in \\text{Herm}(B) \\}
```
and returns the optimal value and the optimizer, Y. This is the dual problem for the SDP
corresponding to the entanglement-assisted commmunication value.
It is related to the channel min-entropy ``H_{\\min}(A|B)_{\\mathcal{J}(\\mathcal{N})}`` by
``\\text{eaCV}(\\mathcal{N}) = 2^{-H_{\\min}(A|B)_{\\mathcal{J}(\\mathcal{N})}}``.
(See [Section 6.1 of this reference](https://arxiv.org/abs/1504.00233) for further details about
the min-entropy.)
Note: we label the primal as the maximization problem.
"""
function eaCVDual(ρ :: AbstractArray, dimA :: Int, dimB :: Int) :: Tuple{Float64,  Matrix{ComplexF64}}
    identMat = Matrix{Float64}(I, dimA, dimA)
    Y = HermitianSemidefinite(dimB)
    objective = real(tr(Y))
    constraint = [kron(identMat , Y) ⪰ ρ]
    problem = minimize(objective,constraint)
    qsolve!(problem)
    return problem.optval, Y.value
end

"""
    pptCV( channel :: Choi, method :: Symbol = :primal )

Numerically evaluates the positive partial transpose (PPT) relaxation of the
communication value for the given [`Choi`](@ref) operator representation of a
quantum channel.
This optimization problem is expressed as,

```math
\\begin{aligned}
    \\text{pptCV}(\\mathcal{N}) &=  \\max_{\\Omega^{AB}} \\text{Tr}[\\Omega^{AB}J_{\\mathcal{N}}] \\\\
    \\hspace{1cm}
    & \\text{s.t.} \\quad \\text{Tr}_A[\\Omega^{AB}] = \\mathbb{I}^B; \\\\
    \\hspace{1cm}
    & \\qquad \\Omega^{AB} \\in \\text{PPT}(A\\;:\\;B),
\\end{aligned}
```

where ``\\text{PPT}(A \\;:\\; B)`` indicates the cone of matrices with a positve
partial transpose.
The ``\\text{PPT}(A\\;:\\;B)`` cone contains the ``\\text{SEP}(A\\;:\\;B)`` cone
therfore, the ``\\text{pptCV}`` is an upper bound for the ``\\text{cv}``

```math
\\text{cv}(\\mathcal{N}) \\leq \\text{pptCV}(\\mathcal{N}),
```

where equality is found for the channels whose separable cone is the same as their
PPT cone.

The `method` parameter accepts the values:
* `:primal` - Solves the primal problem (see [`pptCVPrimal`](@ref)).
* `:dual` - Solves the dual problem (see [`pptCVDual`](@ref)).

# Returns

A Tuple containing `(max_cv, optimizers...)` where `max_cv` is the communication
value and optimizers is the optimal value for each of the optimization paraameters.
"""
pptCV(channel :: Choi, method::Symbol=:primal) = pptCV(channel, Val(method))
pptCV(channel :: Choi, ::Val{:primal}) = pptCVPrimal(channel.JN, channel.in_dim, channel.out_dim)
pptCV(channel :: Choi, ::Val{:dual}) = pptCVDual(channel.JN, channel.in_dim, channel.out_dim)

"""
    pptCVPrimal(
        ρ :: Matrix{<:Number},
        dimA :: Int,
        dimB :: Int
    ) :: Tuple{Float64,  Matrix{ComplexF64}}

This function solves the SDP
```math
\\max \\{ \\langle \\rho, X \\rangle :  \\text{Tr}_{A}(X) = I_{B} , \\Gamma(X) \\succeq 0, X \\succeq 0 \\}
```
where ``\\Gamma( \\cdot)`` is the partial transpose with respect to the second system, and
returns the optimal value and the optimizer, X.
This is the primal problem for the SDP relaxation of the channel value. The relaxation is to the PPT cone.
This has various interpretations. Note: we label the primal as the maximization problem.
"""
function pptCVPrimal(ρ :: AbstractArray, dimA :: Int, dimB :: Int) :: Tuple{Float64,  Matrix{ComplexF64}}
    X = HermitianSemidefinite(dimA*dimB)
    objective = real(tr(ρ' * X))
    constraints = [partialtrace(X, 1, [dimA,dimB]) == Matrix{Float64}(I,dimB,dimB),
                   partialtranspose(X,2,[dimA,dimB]) ⪰ 0]
    problem = maximize(objective,constraints)
    qsolve!(problem)
    return problem.optval, X.value
end

"""
    pptCVDual(
        ρ :: Matrix{<:Number},
        dimA :: Int,
        dimB :: Int
    ) :: Tuple{Float64,  Matrix{ComplexF64}, Matrix{ComplexF64}}

This function solves the SDP
```math
\\min \\{ \\text{Tr}(Y_{1}) : I_{A} \\otimes Y_{1} - \\Gamma(Y_{2}) \\succeq \\rho, Y_{2} \\succeq 0, Y_{1} \\in \\text{Herm}(B) \\}
```
where `` \\Gamma( \\cdot)`` is the partial transpose with respect to the second system,
and returns the optimal value and optimizer, ``(Y_1 , Y_2 )``.
This is the dual problem for the SDP relaxation of the channel value. The relaxation is to the PPT cone.
This has various interpretations. Note: we label the primal as the maximization problem.
"""
function pptCVDual(ρ :: AbstractArray, dimA :: Int, dimB :: Int, dual=true :: Bool) :: Tuple{Float64,  Matrix{ComplexF64}, Matrix{ComplexF64}}
    identMat = Matrix{Float64}(I, dimA, dimA)
    Y1 = ComplexVariable(dimB,dimB)
    Y2 = HermitianSemidefinite(dimA*dimB)
    objective = real(tr(Y1))
    constraints = [kron(identMat,Y1) - partialtranspose(Y2, 2 , [dimA,dimB]) ⪰ ρ,
                   Y1' - Y1 == zeros(dimB,dimB)] #Forces Hermiticity
    problem = minimize(objective,constraints)
    qsolve!(problem)
    return problem.optval, Y1.value, Y2.value
end

"""
    pptCVMultiplicativity(
        channel1 :: Choi,
        channel2 :: Choi;
        singular_method::Symbol = :primal,
        parallel_method::Symbol = :dual
    ) :: Vector

This function takes the [`Choi`](@ref) operators of two channels
`channel1` (``\\mathcal{N}_{A_{1} \\to B_{1}}``) and `channel2`
(``\\mathcal{M}_{A_{2} \\to B_{2}}``) and returns as an array:

```
[
    pptCV( channel1 ),
    pptCV( channel2 ),
    pptCV( channel1 ⊗ channel2 ),
    pprt( channel1 ⊗ channel2 ) - pptCV( channel1 ) * pptCV( channel2 )
]
```

By default, it uses [`pptCVPrimal`](@ref) for the single channel uses, as this
provides a lower bound, and [`pptCVDual`](@ref) for the parallel case, as
this is always an upper bound.
These defaults can be overridden with the keyword args `singular_method` and
`parallel_method` which each accept the symbol values `:primal` and `:dual`.
"""
function pptCVMultiplicativity(
    channel1 :: Choi,
    channel2 :: Choi;
    singular_method::Symbol = :primal,
    parallel_method::Symbol = :dual
) :: Vector
    # singular channel use
    cv_1, = pptCV(channel1, singular_method)
    cv_2, = pptCV(channel2, singular_method)

    # parallel channel use
    par_cv, = pptCV(parChoi(channel1, channel2), parallel_method)

    return [cv_1,cv_2,par_cv, par_cv - (cv_1 * cv_2)]
end

"""
    pptCVMultiplicativity(
            JN :: Matrix,
            Ndin :: Int,
            Ndout :: Int,
            JM :: Matrix,
            Mdin :: Int,
            Mdout :: Int;
            step1isdual = false :: Bool,
            step2isprimal = false :: Bool
    ) :: Vector

This function takes the Choi operators of two channels
``\\mathcal{N}_{A_{1} \\to B_{1}}`` and ``\\mathcal{M}_{A_{2} \\to B_{2}}``
along with their input and output dimensions and returns ``\\text{pptCV}(\\mathcal{N})``,
``\\text{pptCV}(\\mathcal{M})``, and ``\\text{pptCV}(\\mathcal{N}\\otimes \\mathcal{M})``.
By default, it uses [`pptCVPrimal`](@ref) for the single channel values, as this
provides a lower bound, and [`pptCVDual`](@ref) for the parallel case, as
this is always an upper bound. If the dimension is such that the dual can't be
used, there is an optional argument for using [`pptCVPrimal`](@ref). There is
also an optional argument to use [`pptCVDual`](@ref) for single channel values.
"""
function pptCVMultiplicativity(
    JN :: Matrix,
    Ndin :: Int,
    Ndout :: Int,
    JM :: Matrix,
    Mdin :: Int,
    Mdout :: Int;
    step1isdual = false :: Bool,
    step2isprimal = false :: Bool
) :: Vector

    if step1isdual
        cv_1, opt_1 = pptCVPrimal(JN, Ndin, Ndout)
        cv_2, opt_2 = pptCVPrimal(JM, Mdin, Mdout)
    else
        cv_1, opt_11, opt_21 = pptCVDual(JN, Ndin, Ndout)
        cv_2, opt_12, opt_22 = pptCVDual(JM, Mdin, Mdout)
    end

    perm_vec = [1,3,2,4]
    dimVec = [Ndin, Ndout, Mdin, Mdout]
    tot_din = Ndin*Mdin
    tot_dout = Ndout*Mdout
    par_choi = permuteSubsystems(kron(JN,JM),perm_vec,dimVec)
    if step2isprimal
        par_cv, opt1 = pptCVPrimal(par_choi,tot_din,tot_dout)
    else
        par_cv, opt1, opt2 = pptCVDual(par_choi,tot_din,tot_dout)
    end

    return [cv_1,cv_2,par_cv, par_cv - (cv_1 * cv_2)]
end

"""
    WHIDLP(
        d1 :: Int64,
        d2 :: Int64,
        λ :: Union{Int,Float64}
    ) :: Float64

This function implements the linear program (LP) to determine the communication
value of the Werner-Holevo channel tensored with the identity channel, when the
problem is relaxed to optimizing over the PPT cone. (See cite for derivation).
d1,d2 are the input-output dimensions of the Werner-Holevo and identity channel
respectively. λ is the parameter defining the Werener-Holevo channel.
"""
function WHIDLP(d1 :: Int64, d2 :: Int64, λ :: Union{Int,Float64} ) :: Float64
    #This is the vector of variables and we use alphabetical order x[1] = w, x[2] = x,...
    v = Variable(4)
    #Given the number of constraints we define the problem and add constraints
    objective = (v[1]+v[3]*d2+(2λ-1)*(v[2]+v[4]*d2))
    problem = maximize(objective)
    problem.constraints += [
        0 <= v[1]-v[2]+d2*v[3]-d2*v[4],
        0 <= v[1]-v[2],
        0 <= v[1]+v[2]+d2*v[3]+d2*v[4],
        0 <= v[1]+v[2],
        0 <= v[1] + d1*v[2] - v[3] - d1*v[4],
        0 <= v[1] - v[3],
        0 <= v[1] + d1*v[2] + v[3] + d1*v[4],
        0 <= v[1] + v[3],
        1 == d1*d2*v[1]+d2*v[2]+d1*v[3] + v[4]
    ]
    qsolve!(problem)

    return d1*d2*problem.optval
end

"""
    generalWHLPConstraints(
        n :: Int,
        d :: Int,
        λ_vec :: Union{Vector{Float64},Vector{Int64}}
    ) :: Tuple{Matrix{Float64},Matrix{Float64},Matrix{Float64},Matrix{Float64}}

This function returns the linear program constraints for calculating the PPT
communication value of the Werner-Holevo channels run in parallel for arbitrary
``n``. ``n`` is the number of Werner-Holevo channels, ``d`` is the dimension of
every Werner-Holevo channel (assumed to be the same), and ``\\lambda_{\\text{vec}}``
is such that ``\\lambda_{\\text{vec}}[i]`` is the λ parameter for the ``i^{th}``
Werner-Holevo channel. The returned matrices represent the linear maps enforcing
the constraints on the optimizer. ``A`` is the poistivity constraint, ``B`` is
the PPT constraint, ``g`` is the trace constraint, and ``a`` defines the objective
function.

!!! warning
    It takes ``O(n2^{2n})`` steps to generate. If one wants a large dimension,
    we suggest you save the resulting constraints.
"""
function generalWHLPConstraints(n :: Int, d :: Int, λ_vec :: Union{Vector{Float64},Vector{Int64}}) :: Tuple{Matrix{Float64},Matrix{Float64},Matrix{Float64},Matrix{Float64}}
    #Quick sanity checks
    if n > 11
        println("WARNING: You are trying to generate constraints at a size where the time it will take is non-trivial.")
    elseif length(λ_vec) != n
        throw(DomainError(λ_vec, "λ_vec must have the length of n."))
    elseif !all(λ_vec -> λ_vec >= 0 && λ_vec <= 1, λ_vec)
        throw(DomainError(λ_vec, "λ_vec must contain values in [0,1]"))
    end

    #Note that the objective function will need to be scaled by d after the calculation
    A = zeros(2^n,2^n)
    B = zeros(2^n,2^n)
    g = zeros(1,2^n)
    #a = zeros(2^n,1)
    ζ = ones(2^n,1)
    for s in [0:2^n - 1;]
        #Trace condition
        s_string = digits(Int8, s, base=2, pad=n) |> reverse
        w_s = sum(s_string) #This is the hamming weight
        g[s+1] = d^(n-w_s) #This is d to the power of the hamming weight
        for j in [0:2^n - 1;]
            j_string = digits(Int8, j, base=2, pad=n) |> reverse
            #This is the positivity condition
            #This is the hamming weight of bit and of s and j
            w_sj = sum(digits(Int8,(s&j),base=2,pad=n))
            A[j+1,s+1] = _WH_coeff_sign(w_sj)
            #This is for the objective function
            if j <= n -1
                ζ[s+1] = ζ[s+1]*_WH_lambda_coeff(j+1,s_string[j+1],λ_vec)
            end
            #This whole loop is the ppt constraint
            B_nonzero = true
            for i in [0:n-1;]
                #This is for the ppt constraint
                if s_string[i+1] == 1 && j_string[i+1] == 0
                    B_nonzero = false
                end
            end
            if B_nonzero
                B[j+1,s+1] = d^(w_s) #d^(-1. *(n - w_s)) original scaling
            else
                B[j+1,s+1] = 0
            end
        end
    end
    a = A*ζ
    return A, B, g, a
end
#These are helper functions
function _WH_coeff_sign(x)
    isodd(x) ? -1 : 1
end
function _WH_lambda_coeff(i,bit,λ_vec)
    return bit == 0 ? λ_vec[i] : (1-λ_vec[i])
end

"""
    wernerHolevoCVPPT(
        n :: Int64
        d :: Int64,
        A :: Matrix{Float64},
        B :: Matrix{Float64},
        g :: Matrix{Float64},
        a :: Matrix{Float64}
    ):: Tuple{Float64, Matrix{Float64}}

This function evaluates the linear program for the PPT relaxation of the communication
value of the Werner-Holevo channel. The LP is written
```math
    \\max \\{\\langle a, v \\rangle : Ax \\geq 0 , Bx \\geq 0 , \\langle g , v \\rangle = 1 \\}
```
This function takes as inputs: ``n``, the number of Werner-Holevo channels,
``d``, the dimension of every Werner-Holevo Channel, and the constraints
``A,B,g,a`` which are obtained from [`generalWHLPConstraints`](@ref) outputs.
It returns the cvPPT value and the optimizer.

!!! warning
    For ``n \\geq 10`` the solver may be slow.
"""
function wernerHolevoCVPPT(
        n :: Int64,
        d :: Int64,
        A :: Matrix{Float64},
        B :: Matrix{Float64},
        g :: Matrix{Float64},
        a :: Matrix{Float64}
    ) :: Tuple{Float64, Matrix{Float64}}

    v = Variable(2^n)
    objective = a' * v
    problem = maximize(objective)
    problem.constraints += [
        A * v >= 0 ; B * v >= 0 ; g * v == 1
    ]
    qsolve!(problem)
    cvPPT = d^n * problem.optval
    return cvPPT, v.value
end
"""
    twoSymCVPrimal(channel :: Choi) :: Tuple{Float64,  Matrix{ComplexF64}}
    twoSymCVPrimal(
        ρ :: Matrix{<:Number},
        dimA :: Int,
        dimB :: Int
    ) :: Tuple{Float64,  Matrix{ComplexF64}}

Given the [`Choi`](@ref) operator representation of a channel, or alternatively,
the Choi matrix `ρ` and its input output dimensions, this function solves the SDP
```math
\\max \\{ \\langle \\rho, X \\rangle :  \\text{Tr}_{A}(X) = I_{B} ,
                                        \\Gamma^{B_{1}}(X) \\succeq 0,
                                        X = (I_{A} \\otimes \\mathbb{F}_{B})X(I_{A} \\otimes \\mathbb{F}_{B}),
                                        X \\succeq 0 \\}
```
where ``\\Gamma^{B_{1}}( \\cdot)`` is the partial transpose with respect to the second system,
and returns the optimal value and the optimizer, X. The conditions on X demand it is two-symmetric, i.e.
an element of the lowest level of the [DPS hierarchy.](https://arxiv.org/abs/quant-ph/0308032)
Note: we label the primal as the maximization problem.

!!! warning "Runs Out of Memory Easily"
    This function  will run out of memory for the tensor product of even qutrit to qutrit channels.
"""
twoSymCVPrimal(channel :: Choi) :: Tuple{Float64,  Matrix{ComplexF64}} = twoSymCVPrimal(channel.JN, channel.in_dim, channel.out_dim)
function twoSymCVPrimal(ρ :: AbstractArray, dimA :: Int, dimB :: Int) :: Tuple{Float64,  Matrix{ComplexF64}}
    X = HermitianSemidefinite(dimA*dimB^2)
    idA = Matrix(1I, dimA, dimA); idB = Matrix(1I, dimB, dimB)
    F = swapOperator(dimB)
    objective = real(tr(kron(ρ,idB)' * X))
    constraints = [partialtrace(partialtrace(X, 1, [dimA,dimB,dimB]),2,[dimB,dimB]) == Matrix{Float64}(I,dimB,dimB),
                   X - (kron(idA,F)*X*kron(idA,F)') == 0,
                   partialtranspose(X,2,[dimA,dimB,dimB]) ⪰ 0]
    problem = maximize(objective,constraints)
    qsolve!(problem)
    return problem.optval, X.value
end
#I don't write the dual program because it's too big to be useful anyways

"""
    threeSymCVPrimal(channel :: Choi) :: Tuple{Float64,  Matrix{ComplexF64}}
    threeSymCVPrimal(
        ρ :: Matrix{<:Number},
        dimA :: Int,
        dimB :: Int
    ) :: Tuple{Float64,  Matrix{ComplexF64}}

Given the [`Choi`](@ref) operator representation of a channel, or alternatively,
the Choi matrix `ρ` and its input output dimensions, this function solves the SDP
```math
\\max \\{ \\langle \\rho, X \\rangle :  \\text{Tr}_{A}(X) = I_{B} ,
                                        \\Gamma^{B_{1}}(X) \\succeq 0,
                                        X = (I_{A} \\otimes W_{\\pi})X(I_{A} \\otimes W_{\\pi}^{*}) \\quad \\forall \\pi \\in \\mathcal{S_{3}},
                                        X \\succeq 0 \\}
```
where ``\\Gamma^{B_{1}}( \\cdot)`` is the partial transpose with respect to the second system,
``W_{\\pi}`` is the unitary that permutes the subspaces according to permutation ``\\pi``.
The function returns the optimal value and the optimizer, X. The conditions on X demand it is two-symmetric, i.e.
an element of the lowest level of the [DPS hierarchy.](https://arxiv.org/abs/quant-ph/0308032)
Note: we label the primal as the maximization problem.
"""
threeSymCVPrimal(channel :: Choi) :: Tuple{Float64,  Matrix{ComplexF64}} = threeSymCVPrimal(channel.JN, channel.in_dim, channel.out_dim)
function threeSymCVPrimal(ρ :: AbstractArray, dimA :: Int, dimB :: Int) :: Tuple{Float64,  Matrix{ComplexF64}}
    """
    This is a helper function for threeSymCVPrimal because permuteSubsystems would need
    a reshape(::Variable, ::NTuple), but Convex only supports reshape(::Variable, ::Int, ::Int).
    Obviously it is not ideal.
    """
    function swapCons(X,dimA,dimB,t)
        idA,idB = Matrix(1I, dimA, dimA), Matrix(1I,dimB,dimB)
        if t == 1 #Permutation [1,2,4,3]
            return kron(idA,kron(idB,swapOperator(dimB)))*X*kron(idA,kron(idB,swapOperator(dimB)))'
        elseif t == 2 #Permutation [1,3,2,4]
            return kron(idA,kron(swapOperator(dimB),idB))*X*kron(idA,kron(swapOperator(dimB),idB))'
        elseif t == 3 #Permutation [1,3,4,2]
            Xp = kron(idA,kron(swapOperator(dimB),idB))*X*kron(idA,kron(swapOperator(dimB),idB))' #Swaps 3 and 2
            return kron(idA,kron(idB,swapOperator(dimB)))*Xp*kron(idA,kron(idB,swapOperator(dimB)))' #Swaps 2 and 4
        elseif t == 4 #Permutation [1,4,2,3]
            Xp = kron(idA,kron(idB,swapOperator(dimB)))*X*kron(idA,kron(idB,swapOperator(dimB)))' #Swaps 3 and 4
            X = kron(idA,kron(swapOperator(dimB),idB))*Xp*kron(idA,kron(swapOperator(dimB),idB))' #Swaps 2 and 4
        else #Permutation [1,4,3,2]
            Xp = kron(idA,kron(swapOperator(dimB),idB))*X*kron(idA,kron(swapOperator(dimB),idB))' #Swaps 3 and 2
            Xpp = kron(idA,kron(idB,swapOperator(dimB)))*Xp*kron(idA,kron(idB,swapOperator(dimB)))' #Swaps 2 and 4
            return kron(idA,kron(swapOperator(dimB),idB))*Xpp*kron(idA,kron(swapOperator(dimB),idB))' #Swaps 3 and 4
        end
    end
    X = HermitianSemidefinite(dimA*dimB^3)
    idB12 = Matrix(1I, dimB^2, dimB^2)
    dims = [dimA,dimB,dimB,dimB]
    objective = real(tr(kron(ρ,idB12)' * X))
    constraints = [partialtrace(partialtrace(X, 1, [dimA,dimB,dimB^2]),2,[dimB,dimB^2]) == Matrix{Float64}(I,dimB,dimB),
                   X - swapCons(X,dimA,dimB,1) == 0,
                   X - swapCons(X,dimA,dimB,2) == 0,
                   X - swapCons(X,dimA,dimB,3) == 0,
                   X - swapCons(X,dimA,dimB,4) == 0,
                   X - swapCons(X,dimA,dimB,5) == 0,
                   partialtranspose(X,2,dims) ⪰ 0]
    problem = maximize(objective,constraints)
    qsolve!(problem)
    return problem.optval, X.value
end
#I don't write the dual program because it's too big to be useful anyways

max_forks_repo_name: ChitambarLab/cv-channel


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    GSLIB

Rotation convention of the GSLIB software.
"""
struct GSLIB <: RotationConvention end
axesseq(::Type{GSLIB}) = :ZXY
orientation(::Type{GSLIB}) = (:CW,:CCW,:CCW)
angleunits(::Type{GSLIB}) = :DEG
mainaxis(::Type{GSLIB}) = :Y
isextrinsic(::Type{GSLIB}) = false


text:
using NormalisingFlows, Base.Test, Optim, Nabla, QuadGK, Distributions

import PDMats
import NormalisingFlows: apply, logdetJ

const __transforms = [DiagAffine, Affine, Planar, Radial]

@testset "NormalisingFlows" begin
    include("normal.jl")
    include("invertibles.jl")
    include("flows.jl")
end


text:
using ASCIIPlots
using Base.Test

A = eye(10)
p = imagesc(A)
@test p.s == "\n\n\t@#                  \n\t  @#                \n\t    @#              \n\t      @#            \n\t        @#          \n\t          @#        \n\t            @#      \n\t              @#    \n\t                @#  \n\t                  @#\n\n"

A = eye(20)
imagesc(A)

A = eye(10)
imagesc(A)
A[10, 1] = 1
A[1, 10] = 1
imagesc(A)

A = randn(10, 10)
imagesc(A)


text:
using Pkg
packages = keys(Pkg.installed())
if !in("Vec", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/Vec.jl.git"))
end
if !in("Records", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/Records.jl.git"))
end
if !in("AutomotiveSimulator", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/AutomotiveSimulator.jl.git"))
end
if !in("AutomotiveVisualization", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/AutomotiveVisualization.jl.git"))
end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
module TranscodingStreams

export
    TranscodingStream,
    Noop,
    NoopStream

const ByteData = Union{Vector{UInt8},Base.CodeUnits{UInt8}}

include("memory.jl")
include("buffer.jl")
include("error.jl")
include("codec.jl")
include("state.jl")
include("stream.jl")
include("io.jl")
include("noop.jl")
include("transcode.jl")
include("testtools.jl")

end # module


text:
__precompile__(true)

"""
The module for [Xorshift Family](@ref).

Provide 11 RNG types:

- [`Xorshift64`](@ref)
- [`Xorshift64Star`](@ref)
- [`Xorshift128`](@ref)
- [`Xorshift128Star`](@ref)
- [`Xorshift128Plus`](@ref)
- [`Xorshift1024`](@ref)
- [`Xorshift1024Star`](@ref)
- [`Xorshift1024Plus`](@ref)
- [`Xoroshiro128`](@ref)
- [`Xoroshiro128Star`](@ref)
- [`Xoroshiro128Plus`](@ref)
"""
module Xorshifts

export Xorshift64, Xorshift64Star
include("xorshift64.jl")

export Xorshift128, Xorshift128Star, Xorshift128Plus
include("xorshift128.jl")

export Xorshift1024, Xorshift1024Star, Xorshift1024Plus
include("xorshift1024.jl")

export Xoroshiro128, Xoroshiro128Star, Xoroshiro128Plus
include("xoroshiro128.jl")

include("docs.jl")

end


text:
using Pkg
Pkg.add(Pkg.PackageSpec(; url="https://github.com/JuliaGaussianProcesses/JuliaGPsDocs.jl"))

## Build the docs
using AugmentedGPLikelihoods
using JuliaGPsDocs

JuliaGPsDocs.generate_examples(AugmentedGPLikelihoods)

using Documenter
using DocumenterCitations
using Literate

DocMeta.setdocmeta!(
    AugmentedGPLikelihoods, :DocTestSetup, :(using AugmentedGPLikelihoods); recursive=true
)

bib = CitationBibliography(joinpath(@__DIR__, "references.bib"))

makedocs(
    bib;
    modules=[AugmentedGPLikelihoods],
    authors="Théo Galy-Fajou <theo.galyfajou@gmail.com> and contributors",
    repo="https://github.com/JuliaGaussianProcesses/AugmentedGPLikelihoods.jl/blob/{commit}{path}#{line}",
    sitename="AugmentedGPLikelihoods.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaGaussianProcesses.github.io/AugmentedGPLikelihoods.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Likelihoods" => map(readdir(joinpath(@__DIR__, "src", "likelihoods"))) do x
            joinpath("likelihoods", x)
        end,
        "Examples" => JuliaGPsDocs.find_generated_examples(AugmentedGPLikelihoods),
        "Misc" => "misc.md",
        "References" => "references.md",
    ],
    strict=true,
)

deploydocs(;
    repo="github.com/JuliaGaussianProcesses/AugmentedGPLikelihoods.jl",
    devbranch="main",
    push_preview=true,
)


text:
using BinDeps
using Pkg

@BinDeps.setup

url = "https://github.com/ajaymendez/rocksdb/archive/c_checkpoint.zip"

librocksdb = library_dependency("librocksdb")

provides(Sources, URI(url), librocksdb, unpacked_dir="rocksdb-c_checkpoint")

builddir = BinDeps.builddir(librocksdb)
srcdir = joinpath(BinDeps.depsdir(librocksdb),"src", "rocksdb-c_checkpoint")
libdir = BinDeps.libdir(librocksdb)

if Sys.isapple()
    libfile = joinpath(libdir,librocksdb.name*".dylib")
    libname = "librocksdb.dylib"
elseif Sys.isunix()
    libfile = joinpath(libdir,librocksdb.name*".so")
    libname = "librocksdb.so"
end

provides(BuildProcess,
    (@build_steps begin
        GetSources(librocksdb)
        CreateDirectory(libdir)
        @build_steps begin
            ChangeDirectory(srcdir)
            FileRule(libfile, @build_steps begin
                     `make shared_lib`
                     `cp $(libname)  ../../usr/lib/`
            end)
        end
    end), librocksdb, os = :Unix)

@BinDeps.install Dict(:librocksdb => :librocksdb)


text:
# This interface the data stored in data/processed/rath2017__data fetched from
# Rath 2017 (https://pure.mpg.de/pubman/item/item_2508673_4)

# A single structure with all the data.
const _RATH_BUNDLE = Dict()
function _load_rath_bundle()

    files = [_cont_culture_proc_file(exp) for exp in EXPS]
    push!(files, _stand_medium_proc_file())
    for file in files
        !isfile(file) && return _RATH_BUNDLE
    end

    empty!(_RATH_BUNDLE)
    for exp in EXPS

        data = get!(_RATH_BUNDLE, exp, Dict())

        # std medium
        stdm_conv = CSV.read(_stand_medium_proc_file(), DataFrame);
        for (i, met) in enumerate(stdm_conv[!, :id])
            val = stdm_conv[i, :conc]
            err = 0.0
            unit = stdm_conv[i, :unit]
            # id,                   val,             err,      unit
            data["c$(met)"] = Dict("val" => val, "err" => err, "unit" => unit)
        end

        #cont cul data
        cul_data_convs = CSV.read(_cont_culture_proc_file(exp), DataFrame)
        for (i, id) in enumerate(cul_data_convs[!, :id])
            val = cul_data_convs[i, :val]
            err = cul_data_convs[i, :err]
            unit = cul_data_convs[i, :unit]
            # id,       val, err, unit
            data[id] = Dict("val" => val, "err" => err, "unit" => unit)
        end

        # ξ
        ξval = data["Xv"]["val"] / data["D"]["val"]
        ξerr = ξval * (data["Xv"]["err"]/ data["Xv"]["val"])
        ξunit = "gDW * hr/ L"
        data["ξ"] = Dict("val" => ξval, "err" => ξerr, "unit" => ξunit)

    end
    return _RATH_BUNDLE
end

# interface
# Will define the val, qval, sval, sval, err and unit function to access Rath data
_parse_id(prefix, id) = (string(prefix) == "q" && string(id) == "μ") ? "μ" : string(prefix, id) # handle biomass
function _define_interface()
    for base_fun in [:val, :err, :unit]
        key = string(base_fun)
        @eval begin
            $base_fun(id, exp) = _RATH_BUNDLE[string(exp)][string(id)][$key]
            function $base_fun(id, exp::AbstractString, deflt)
                try
                    $base_fun(id, exp);
                catch err
                    err isa KeyError && return deflt
                    rethrow(err)
                end
            end
            $base_fun(id, exps::Vector, deflt) = [$base_fun(id, exp, deflt) for exp in exps]
            $base_fun(id, exps::Vector) = [$base_fun(id, exp) for exp in exps]
            $base_fun(id) = $base_fun(id, $EXPS) # all the experiments
        end

        for prefix in ["q", "c", "s"]
            fun = Symbol(string(prefix, base_fun))
            @eval $fun(id) = $base_fun(_parse_id($prefix, id))
            @eval $fun(id, exp) = $base_fun(_parse_id($prefix, id), exp)
            @eval $fun(id, exps, deflt) = $base_fun(_parse_id($prefix, id), exps, deflt)
        end
    end
end


text:
using Plots
using FFTW

cos_taylor(x) = 1-( x^2/2 )+( x^4/24.0 )-( x^6/720.0 )+( x^8/40320.0 );

"""
Calculate cos(x) with a five-term Taylor series, using only the quarter-waves closest to the origin for increased accuracy.
"""
function fast_cos( x::Float64 )
        if x < -π/2
            x = -π/2 + abs(x+π/2);
            return -cos_taylor(x);
        elseif π/2 < x
            x = π/2 - abs(x-π/2);
            return -cos_taylor(x)
        else
            return cos_taylor(x)
        end
end

x = LinRange(-π, π, 1000);

cosx_taylor = cos_taylor.(x);
cosx_native = cos.(x);

delta = cosx_native-cosx_taylor

# Time domain comparison.
plot( 	x, cosx_native,
	 	label = "Julia native cosine",
        size = (1024, 720),)
plot!( 	x, cosx_taylor,
	 	label = "5-term Taylor series",)
plot!( 	x, delta*100,
	 	label = "Difference*100",)

savefig("figures/cos/time_domain")

taylor_cos2 = fast_cos.(x);
delta2 = cosx_native-taylor_cos2

# Taylor vs first quarter taylor.
plot( 	x, delta*100,
	 	label = "Difference",
        size = (1024, 720),)
plot!( 	x, delta2*100,
	 	label = "Difference (first quarter)",)

savefig("figures/cos/time_domain_comparison2")

"""
Wrap sample to a -π =< sample < π range.
"""
function wrap_phase( sample::Float64 )
    while π <= sample
        sample = sample - 2π;
    end
    while sample < -π
        sample = sample + 2π;
    end
    return sample;
end

# Frequency domain analysis
x = LinRange(0, 1024π, 2^14);
x = wrap_phase.(x);
cosx_taylor = cos_taylor.(x);
cosx_native = cos.(x);

"""
Calculate normalized power spectrum for real signals.
"""
function frequecny_domain_analysis( signal_t )
    power_f = abs.(fft(signal_t)).^2;
    normalized_power = power_f/maximum(power_f);
    power_dB = 10*log10.(normalized_power);
    return power_dB[1:Int(size(power_dB)[1]/2)]
end

omega = LinRange(0, 1, Int(size(x)[1]/2));

# Frequency domain.
plot( 	omega, frequecny_domain_analysis(cosx_native),
	 	label = "Julia native cosine",
        ylims = [-100, 0],
        size = (1024, 720),)
xlabel!("Frequency relative to Fs")
ylabel!("Normalized power [dBC]")

plot!( 	omega, frequecny_domain_analysis(cosx_taylor),
	 	label = "5-term Taylor series",)

savefig("figures/cos/frequency_domain")

# Frequency domain first quarter.
plot( 	omega, frequecny_domain_analysis(cosx_native),
	 	label = "Julia native cosine",
        ylims = [-100, 0],
        size = (1024, 720),)
xlabel!("Frequency relative to Fs")
ylabel!("Normalized power [dBC]")

plot!( 	omega, frequecny_domain_analysis(fast_cos.(x)),
	 	label = "5-term Taylor, first quarter",)

savefig("figures/cos/frequency_domain2")

max_stars_repo_name: ErikBuer/Integer-Array

max_issues_repo_name: ErikBuer/Integer-Array

max_forks_repo_name: ErikBuer/Integer-Array


text:
# Tabix Reader
# ============

function Base.read(filename::AbstractString, ::Type{Tabix})
    return open(input -> read(input, Tabix), filename)
end

function Base.read(input_::IO, ::Type{Tabix})
    input = BGZFStream(input_)

    # magic
    T = read(input, UInt8)
    B = read(input, UInt8)
    I = read(input, UInt8)
    x = read(input, UInt8)
    if T != UInt8('T') || B != UInt8('B') || I != UInt8('I') || x != 0x01
        error("invalid tabix magic bytes")
    end

    n_refs = read(input, Int32)
    format = read(input, Int32)
    col_seq = read(input, Int32)
    col_beg = read(input, Int32)
    col_end = read(input, Int32)
    meta = read(input, Int32)
    skip = read(input, Int32)

    l_nm = read(input, Int32)
    data = read(input, UInt8, l_nm)
    names = split(String(data), '\0', keep=false)

    indexes = read_indexes(input, n_refs)

    if !eof(input)
        n_no_coor = Nullable(read(input, UInt64))
    else
        n_no_coor = Nullable{UInt64}()
    end

    return Tabix(
        format,
        (col_seq, col_beg, col_end),
        meta,
        skip,
        names,
        indexes,
        n_no_coor)
end


text:
export jumble_iter

function jumble_iter(text::String)
    return SubwordIter(word_to_bag(text))
end


text:
using Flux, Flux.Tracker
using Distributions

include("../ADVI/ADVI.jl")

mutable struct State
  w
  m
  s
  State() = new()
end

function State(K::Integer)
  state = State()

  # NOTE: this works! Notice that dimensions do matter
  #       especially for dirichlet stuff. Needs to be a tuple for size.
  #       So, the minimum must be (1, K-1) as the last dim is the dirichlet
  #       dim. Maybe needs a FIXME?
  state.m = ADVI.ModelParam(K, "real")
  state.w = ADVI.ModelParam((1, K - 1), "simplex",
                            m=fill(0.0, 1, K-1),
                            s=fill(exp(-2), 1, K-1))
  state.s = ADVI.ModelParam(K, "positive",
                            m=fill(0.0, K),
                            s=fill(exp(-2), K))

  return state
end

function rsample(state::State)
  real = State()
  tran = State()

  for key in fieldnames(State)
    f = getfield(state, key)

    r = ADVI.rsample(f)
    t = ADVI.transform(f, r)

    setfield!(real, key, r)
    setfield!(tran, key, t)
  end

  return real, tran
end

function loglike(tran::State, y::Vector{Float64})
  K = length(tran.m)
  N = length(y)

  # NOTE: Dirichlet parameters need special treatment because it changes dims!
  w = reshape(tran.w, 1, K)
  m = reshape(tran.m, 1, K)
  s = reshape(tran.s, 1, K)
  y_rs = reshape(y, N, 1)
  return sum(ADVI.lpdf_gmm(y_rs, m, s, w, dims=2, dropdim=true))
end

function logprior(real::State, tran::State, mp::State)
  lp = 0.0

  lp += sum(ADVI.compute_lpdf(Normal(1.85, 5), tran.m))

  lp += sum(ADVI.compute_lpdf(Gamma(1, 0.1), tran.s))
  lp += sum(ADVI.logabsdetJ(mp.s, real.s, tran.s))

  K = length(tran.w)
  lp += sum(ADVI.compute_lpdf(Dirichlet(ones(K) / K), tran.w))
  lp += sum(ADVI.logabsdetJ(mp.w, real.w, tran.w))

  return lp
end

function logq(real::State, mp::State)
  lq = 0.0
  for key in fieldnames(State)
    vp = getfield(mp, key)
    r = getfield(real, key)
    lq += sum(ADVI.log_q(vp, r))
  end
  return lq
end

function compute_elbo(mp::State, y::Vector{Float64}, N::Integer,
                      metrics::Dict{Symbol, Vector{Float64}})
  real, tran = rsample(mp)

  ll = loglike(tran, y) * N / length(y)
  lp = logprior(real, tran, mp)
  lq = logq(real, mp)

  elbo = ll + lp - lq
  append!(metrics[:elbo], elbo.data)

  return elbo
end


text:
struct Foo
    bar
    baz::Int
    qux::Float64
end

foo = Foo("Hello",23,1.5)

@show typeof(foo)
@show fieldnames(foo)
@show foo.bar
@show foo.baz
@show foo.qux

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module FixedEffects

##############################################################################
##
## Dependencies
##
##############################################################################
using Base: @propagate_inbounds
using LinearAlgebra
using StatsBase
using GroupedArrays
using Requires
using Printf

##############################################################################
##
## Load files
##
##############################################################################
include("utils/lsmr.jl")
include("utils/progressbar.jl")

include("FixedEffect.jl")
include("AbstractFixedEffectSolver.jl")
include("FixedEffectSolvers/FixedEffectLinearMap.jl")
include("FixedEffectSolvers/FixedEffectSolverCPU.jl")

has_CUDA() = false
function __init__()
	    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
	    if CUDA.functional()
		    has_CUDA() = true
	    	include("FixedEffectSolvers/FixedEffectSolverGPU.jl")
	    end
	end
end

include("precompile.jl")
_precompile_()

##############################################################################
##
## Exported methods and types
##
##############################################################################

export FixedEffect,
AbstractFixedEffectSolver,
solve_residuals!,
solve_coefficients!

end  # module FixedEffects


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "XML2"
version = v"2.9.10"

# Collection of sources required to build XML2Builder
sources = [
    ArchiveSource("https://github.com/GNOME/libxml2/archive/v$(version).tar.gz",
                  "3bdffb4d728e4dc135b3210bf2a2cebb76548b820a5617c68abb7b83654066dd"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/libxml2-*
./autogen.sh --prefix=${prefix} --host=${target} \
    --without-python \
    --disable-static \
    --with-zlib=${prefix} \
    --with-iconv=${prefix}
make -j${nproc}
make install

# Remove heavy doc directories
rm -rf ${prefix}/share/{doc/libxml2-*,gtk-doc}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libxml2", :libxml2),
    ExecutableProduct("xmlcatalog", :xmlcatalog),
    ExecutableProduct("xmllint", :xmllint),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("Libiconv_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
#############################################################################
# matrixfrac.jl
# implements the atom for x^T*P^{-1}*x, where P is a positive semidefinite
# matrix.
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

struct MatrixFracAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr, AbstractExpr}
    size::Tuple{Int, Int}

    function MatrixFracAtom(x::AbstractExpr, P::AbstractExpr)
        if x.size[2] != 1
            error("first argument of matrix frac must be a vector")
        elseif P.size[1] != P.size[2]
            error("second argument of matrix frac must be square")
        elseif x.size[1] != P.size[1]
            error("sizes must agree for arguments of matrix frac")
        end
        children = (x, P)
        return new(:matrixfrac, hash(children), children, (1,1))
    end
end

function sign(m::MatrixFracAtom)
    return Positive()
end

function monotonicity(m::MatrixFracAtom)
    return (NoMonotonicity(), NoMonotonicity())
end

function curvature(m::MatrixFracAtom)
    return ConvexVexity()
end

function evaluate(m::MatrixFracAtom)
    x = evaluate(m.children[1])
    return x'*inv(evaluate(m.children[2]))*x
end

matrixfrac(x::AbstractExpr, P::AbstractExpr) = MatrixFracAtom(x, P)
matrixfrac(x::Value, P::AbstractExpr) = MatrixFracAtom(Constant(x), P)
matrixfrac(x::AbstractExpr, P::Value) = MatrixFracAtom(x, Constant(P))

function conic_form!(m::MatrixFracAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, m)
        x = m.children[1]
        P = m.children[2]
        t = Variable()
        # the matrix [t x'; x P] has Schur complement t - x'*P^{-1}*x
        # this matrix is PSD <=> t >= x'*P^{-1}*x
        p = minimize(t, [t x'; x P] ⪰ 0)
        cache_conic_form!(unique_conic_forms, m, p)
    end
    return get_conic_form(unique_conic_forms, m)
end


text:
"""
    init_zero_controls(t)

Return initial state of uniformly-zero climate controls.

See also: [`Controls`](@ref)
"""
function init_zero_controls(t::Array{Float64,1})
    c = Controls(
        zeros(size(t)),
        zeros(size(t)),
        zeros(size(t)),
        zeros(size(t))
    )
    return c
end


text:
# # Using custom structs
#
# This tutorial shows how to use custom structs on the GPU. Our example will be a one dimensional
# interpolation. Lets start with the CPU version:
using CUDA

struct Interpolate{A}
    xs::A
    ys::A
end

function (itp::Interpolate)(x)
    i = searchsortedfirst(itp.xs, x)
    i = clamp(i, firstindex(itp.ys), lastindex(itp.ys))
    @inbounds itp.ys[i]
end

xs_cpu = [1.0, 2.0, 3.0]
ys_cpu = [10.0,20.0,30.0]
itp_cpu = Interpolate(xs_cpu, ys_cpu)
pts_cpu = [1.1,2.3]
result_cpu = itp_cpu.(pts_cpu)

# Ok the CPU code works, let's move our data to the GPU:
itp = Interpolate(CuArray(xs_cpu), CuArray(ys_cpu))
pts = CuArray(pts_cpu);
# If we try to call our interpolate `itp.(pts)`, we get an error however:
# ```
# ...
# KernelError: passing and using non-bitstype argument
# ...
# ```
# Why does it throw an error? Our calculation involves
# a custom type `Interpolate{CuArray{Float64, 1}}`.
# At the end of the day all arguments of a CUDA kernel need to
# be bitstypes. However we have
isbitstype(typeof(itp))
# How to fix this? The answer is, that there is a conversion mechanism, which adapts objects into
# CUDA compatible bitstypes.
# It is based on the [Adapt.jl](https://github.com/JuliaGPU/Adapt.jl) package and basic types like `CuArray` already participate in this mechanism. For custom types,
# we just need to add a conversion rule like so:
import Adapt
function Adapt.adapt_structure(to, itp::Interpolate)
    xs = Adapt.adapt_structure(to, itp.xs)
    ys = Adapt.adapt_structure(to, itp.ys)
    Interpolate(xs, ys)
end
# Now our struct plays nicely with CUDA.jl:
result = itp.(pts)
# It works, we get the same result as on the CPU.
@assert CuArray(result_cpu) == result
# Alternatively instead of defining `Adapt.adapt_structure` explictly, we could have done
# ```julia
# Adapt.@adapt_structure Interpolate
# ```
# which expands to the same code that we wrote manually.


text:
#!/usr/bin/env julia

using Base.Threads, FLoops

function _picalc(numsteps)
    slice = 1.0 / numsteps

    @floop @simd for i = 1:numsteps
        x = (i - 0.5) * slice
        @reduce(sum += 4.0 / (1.0 + x ^ 2))
    end

    return sum * slice
end

function picalc(numsteps)

    println("Calculating PI using:")
    println("  ", numsteps, " slices")
    println("  ", nthreads(), " thread(s)")

    start = time()
    mypi = _picalc(numsteps)
    elapsed = time() - start

    println("Obtained value of PI: ", mypi)
    println("Time taken: ", elapsed, " seconds")

end

numsteps = if length(ARGS) > 0
    parse(Int, ARGS[1])
else
    1_000_000_000
end

# Warm up kernel
_picalc(10)

# Run the full example
picalc(numsteps)