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

# Cover bound - an upper bound on the fraction of responses in time for optimal move-up.
# From paper: "A bound on the performance of an optimal ambulance redeployment policy".

# mutates: coverBound, coverBoundSim
function calcCoverBound!(sim::Simulation;
		coverBound::Union{CoverBound,Nothing} = nothing,
		coverBoundSim::Union{CoverBoundSim,Nothing} = nothing, # need this if coverBound.sim is not set
		ambBusyDurationsToSample::Vector{Float} = Float[],
		queuedDurationsToSample::Vector{Float} = Float[],
		doPrint::Bool = false)

	@assert(ambBusyDurationsToSample != [] && all(ambBusyDurationsToSample .> 0))
	@assert(issorted(ambBusyDurationsToSample))

	if coverBound == nothing
		coverBound = initCoverBound(sim, doPrint = doPrint)
		@assert(coverBoundSim != nothing)
		coverBound.sim = coverBoundSim
	end
	cb = coverBound # shorthand
	cbs = coverBoundSim = coverBound.sim # shorthand

	if cb.ambBusyDurationsToSample != ambBusyDurationsToSample || length(cb.ambBusyDurationsToSample) != size(cb.ambBusyDurationProbUpperBounds, 1)
		cb.ambBusyDurationsToSample = ambBusyDurationsToSample
		cb.ambBusyDurationProbUpperBounds = calcAmbBusyDurationProbUpperBounds(sim, cb, doPrint = doPrint)
	else
		# user should ensure that ambBusyDurationProbUpperBounds are for ambBusyDurationsToSample
		doPrint && println("Using existing coverBound.ambBusyDurationProbUpperBounds.")
	end

	calcAmbBusyDurationLowerBoundDistrs!(cb) # create amb busy duration distributions

	cb.accountForQueuedDurations = !isempty(queuedDurationsToSample)
	if cb.accountForQueuedDurations
		if cb.queuedDurationsToSample != queuedDurationsToSample || length(cb.queuedDurationsToSample) != length(cb.queuedDurationsMaxCoverageFrac)
			cb.queuedDurationsToSample = queuedDurationsToSample
			cb.queuedDurationsMaxCoverageFrac = calcQueuedDurationsMaxCoverageFrac(sim, cb, doPrint = doPrint)
		else
			# user should ensure that queuedDurationsMaxCoverageFrac are for queuedDurationsToSample
			doPrint && println("Using existing coverBound.queuedDurationsMaxCoverageFrac.")
		end
	end

	simulateCoverBound!(cb)

	return coverBound
end

# initialise coverBound based on sim
# set kwarg pMedianRelax = true if using linear relaxation of p-median problem, false for binary
function initCoverBound(sim::Simulation; pMedianRelax::Bool = true, doPrint::Bool = false)

	coverBound = CoverBound()

	sim.demand.initialised || initDemand!(sim)
	sim.demandCoverage.initialised || initDemandCoverage!(sim)

	@assert(sim.demand.numSets == 1) # otherwise would have to calculate cover bound for each demand set
	@assert(sim.travel.numSets == 1) # otherwise would have to calculate cover bound for each travel set

	genConfig = readGenConfig(sim.inputFiles["callGenConfig"].path) # used for service time related distributions
	# will assume that all demand priorities have the same distributions

	coverBound.modeLookup = fill(nullIndex, sim.travel.numModes)

	numModes = 0
	for demandPriority in priorities
		responseTravelPriority = sim.responseTravelPriorities[demandPriority]
		travelMode = getTravelMode!(sim.travel, responseTravelPriority, sim.startTime) # assumes sim.travel.numSets == 1

		# check if there is already another coverBoundMode with the same response travel priority
		if coverBound.modeLookup[travelMode.index] != nullIndex
			continue # can skip this demandPriority, it will give the same coverBoundMode as another
		end

		coverBoundMode = initCoverBoundMode(sim, travelMode)
		numModes += 1
		coverBoundMode.index = numModes
		coverBound.modeLookup[travelMode.index] = numModes
		push!(coverBound.modes, coverBoundMode)

		setCoverBoundModeDistrs!(coverBoundMode, genConfig)
		calcApproxDistrs!(coverBoundMode) # pMax::Float = 0.999, dt::Float = 1/60/60/24
	end

	# will use deterministic dispatch delay, to simplify implementation
	dispatchDelayDistr = genConfig.dispatchDelayDistrRng.d
	if var(dispatchDelayDistr) != 0
		@warn("Have changed input dispatch delay distribution, as this implementation of the cover bound assumes deterministic dispatch delay.")
		dispatchDelayDistr = Normal(mean(dispatchDelayDistr), 0)
	end
	coverBound.dispatchDelay = mean(dispatchDelayDistr)
	@assert(coverBound.dispatchDelay >= 0)

	coverBound.numAmbsMaxCoverageFrac = calcNumAmbsMaxCoverageFrac(sim)

	return coverBound
end

function initCoverBoundMode(sim::Simulation, responseTravelMode::TravelMode)
	# shorthand
	currentTime = sim.startTime
	hospitalTravelPriority = lowPriority # default

	# shorthand:
	@unpack stations, numStations = sim
	numNodes = length(sim.net.fGraph.nodes)
	points = sim.demandCoverage.points
	numPoints = length(points)
	nodesPoints = sim.demandCoverage.nodesPoints

	# calculate travel time from every station to every node and point
	travelMode = responseTravelMode # shorthand
	stationsToNodesTimes = zeros(Float, numStations, numNodes)
	stationsToPointsTimes = zeros(Float, numStations, numPoints)
	for i = 1:numStations
		station = stations[i] # shorthand

		# get travel time from station to nearest node
		(node1, dist1) = (station.nearestNodeIndex, station.nearestNodeDist)
		time1 = offRoadTravelTime(travelMode, dist1) # time to reach nearest node from station

		for j = 1:numNodes
			pathTravelTime = shortestPathTravelTime(sim.net, travelMode.index, node1, j)
			stationsToNodesTimes[i,j] = time1 + pathTravelTime
			for k in nodesPoints[j] # indices of points for which node j is the nearest node
				point = points[k] # shorthand
				@assert(point.nearestNodeIndex == j)
				time2 = offRoadTravelTime(travelMode, point.nearestNodeDist) # time to reach point from nearest node
				stationsToPointsTimes[i,k] = time1 + pathTravelTime + time2
			end
		end
	end

	# group nodes based on the order in which each station can reach them (closest to furthest)
	nodeSetsDict = Dict{Vector{Int}, Vector{Int}}() # nodeSetsDict[perm] gives the node indices that can be reached by stations in order stations[perm]
	for j = 1:numNodes
		stationsToNodeTimes = view(stationsToNodesTimes, :, j)
		perm = sortperm(stationsToNodeTimes)
		v = get!(nodeSetsDict, perm, Int[])
		push!(v, j)
	end
	nodeSets = collect(values(nodeSetsDict))
	nodeSetsStationList = collect(keys(nodeSetsDict)) # nodeSetsStationList[i] = list of station indices in order (closest to furthest) that they cover nodeSets[i]

	# calculate travel time from every demand point to the closest hospital
	travelMode = getTravelMode!(sim.travel, hospitalTravelPriority, currentTime) # assumes travel.numSets == 1
	pointsToHospitalTimes = zeros(Float, numPoints)
	hospitalToNearestNodeTime = [offRoadTravelTime(travelMode, h.nearestNodeDist) for h in sim.hospitals]
	for j = 1:numNodes
		hospitalIndex = travelMode.fNetTravel.fNodeNearestHospitalIndex[j]
		node2 = sim.hospitals[hospitalIndex].nearestNodeIndex
		pathTravelTime = shortestPathTravelTime(sim.net, travelMode.index, j, node2)
		time2 = hospitalToNearestNodeTime[hospitalIndex]
		for k in nodesPoints[j] # indices of points for which node j is the nearest node
			point = points[k] # shorthand
			time1 = offRoadTravelTime(travelMode, point.nearestNodeDist) # time to reach point from nearest node
			@assert(pointsToHospitalTimes[k] == 0) # value should not yet be initialised
			pointsToHospitalTimes[k] = time1 + pathTravelTime + time2
		end
	end

	coverBoundMode = CoverBoundMode()
	coverBoundMode.stationsToNodesTimes = stationsToNodesTimes
	coverBoundMode.stationsToPointsTimes = stationsToPointsTimes
	coverBoundMode.pointsToHospitalTimes = pointsToHospitalTimes
	coverBoundMode.nodeSets = nodeSets
	coverBoundMode.nodeSetsStationList = nodeSetsStationList
	coverBoundMode.responseTravelMode = responseTravelMode

	return coverBoundMode
end

# mutates: coverBoundMode.distrs, coverBoundMode.transportProb
function setCoverBoundModeDistrs!(coverBoundMode::CoverBoundMode, genConfig::GenConfig)
	distrs = coverBoundMode.distrs # shorthand
	distrs["onSceneDuration"] = genConfig.onSceneDurationDistrRng.d
	distrs["transport"] = genConfig.transportDistrRng.d
	distrs["handoverDuration"] = genConfig.handoverDurationDistrRng.d
	coverBoundMode.transportProb = distrs["transport"].p # probability from bernoulli distribution
end

# mutates: coverBoundMode.distrs, coverBoundMode.transportProb
function setCoverBoundModeDistrs!(coverBoundMode::CoverBoundMode, sim::Simulation)
	genConfig = readGenConfig(sim.inputFiles["callGenConfig"].path)
	setCoverBoundModeDistrs!(coverBoundMode, genConfig)
end

# Create a discrete, non-parametric distribution that approximates the given distribution at the given x values.
# Set lowerBound = true (default) to create distribution that underestimates values, lowerBound = false to overestimate.
# Set removeTrailingZeros = true (default is false) to remove zeros at start and end of cdf of distribution.
function approxDistr(d::Distribution, x::Vector{Float}; lowerBound::Bool = true, removeTrailingZeros::Bool = false)
	x = issorted(x) ? x : sort(x)
	P = [cdf(d,v) for v in x] # probability of sampling value <= v for each v in x
	p = lowerBound ? (vcat(P[2:end],1) - P) : (P - vcat(0,P[1:end-1])) # convert cdf to pdf
	if removeTrailingZeros
		i = 1
		while p[i] == 0 i += 1 end
		j = length(p)
		while p[j] == 0 j-= 1 end
		x = x[i:j]
		p = p[i:j]
	end
	return DiscreteNonParametric(x,p)
end

# Convolute two discrete, non-parametric distributions.
# Requires that x values of the distributions have same step sizes.
function convolute(d1::DiscreteNonParametric, d2::DiscreteNonParametric)
	# shorthand
	x1, p1 = d1.support, d1.p
	x2, p2 = d2.support, d2.p
	n1, n2 = length(x1), length(x2)

	if n1 == 1 || n2 == 1
		x = x1 .+ x2
		p = p1 .* p2
		return DiscreteNonParametric(x,p)
	end

	# check that distribution of x values in the distributions match
	function getDistributionStepSize(d::DiscreteNonParametric)
		x = d.support
		@assert(issorted(x))
		n = length(x)
		dx = (x[end] - x[1]) / (n - 1)
		@assert(dx > 0)
		for i = 1:n @assert(isapprox(x[i], x[1] + dx*(i-1))) end
		return dx
	end
	dx1 = getDistributionStepSize(d1)
	dx2 = getDistributionStepSize(d2)
	@assert(isapprox(dx1, dx2))
	dx = (dx1 + dx2) / 2

	# calculate x for convolution
	xMin = x1[1] + x2[1]
	xMax = x1[end] + x2[end]
	nx = round(Int, (xMax - xMin) / dx) + 1
	x = collect(range(xMin; stop = xMax, length = nx))
	@assert(isapprox((x[end] - x[1]) / (nx - 1), dx))

	# calculate p for convolution
	p = zeros(Float, nx) # note: calculating p from p1 * p2' is inefficient
	for i = 1:n1, j = 1:n2
		p[i+j-1] += p1[i] * p2[j]
	end
	@assert(isapprox(sum(p), 1))

	return DiscreteNonParametric(x,p)
end

# Calculate distributions needed for lower bound on amb busy duration distribution.
# Mutates: coverBoundMode.approxDistrs
function calcApproxDistrs!(coverBoundMode::CoverBoundMode; pMax::Float = 0.999, dt::Float = 1/60/60/24)
	@assert(0 <= pMax < 1)
	@assert(dt > 0)

	function approxDistr2(d::Distribution, pMax::Float, dx::Float; kwargs...)
		@assert(0 <= pMax < 1)
		@assert(dx > 0)
		xMax = quantile(d, pMax)
		xMax = ceil(xMax / dx) * dx
		x = collect(0:dx:xMax)
		return approxDistr(d, x; kwargs...)
	end

	kwargs = (lowerBound = true, removeTrailingZeros = true)
	approxDistrs = coverBoundMode.approxDistrs # shorthand
	for (name, distr) in coverBoundMode.distrs
		if name == "transport" continue end
		approxDistrs[name] = approxDistr2(distr, pMax, dt; kwargs...)
	end

	approxDistrs["d1"] = deepcopy(approxDistrs["onSceneDuration"]) # approximation of on-scene duration
	approxDistrs["d2"] = convolute(approxDistrs["d1"], approxDistrs["handoverDuration"]) # lower bound approximation of d1 and handover duration
	# Could create a weighted sum of distributions d1 and d2 so that there is only one amb busy duration distribution to sample from,
	# but that is more work (for me) than just having the two distributions and accounting for the relative probability of using each.

	return nothing
end

function makeCdf(d::DiscreteNonParametric)
	@assert(issorted(d.support))
	y = cumsum(d.p)
	@assert(isapprox(y[end], 1))
	return y
end

# return indices (i,j) such that xs[i] <= x <= xs[j], for increasing xs
# if x is outside range of xs, then i = 0 or j = length(xs)+1
function binarySearch(xs::Vector{Float}, x::Float)::Tuple{Int,Int}
	# @assert(issorted(xs, lt=<=)) # values should be strictly increasing # slow
	n = length(xs)

	# edge cases
	if x < xs[1] return (0,1)
	elseif x == xs[end] return (n,n)
	elseif x > xs[end] return (n,n+1)
	end

	i = 1; j = n
	while i < j-1
		k = div(i+j, 2)
		xs[k] <= x ? i = k : j = k
	end

	return xs[i] == x ? (i,i) : (i,j)
end

# Given distribution and its cdf (vector of increasing probability values) and x,
# find values (floor and ceil) in domain (d.support) that are closest to x and return corresponding cdf values.
function binarySearch(d::DiscreteNonParametric, cdf::Vector{Float}, x::Float)::Tuple{Float,Float}
	n = length(cdf)
	@assert(length(d.p) == n)
	(i,j) = binarySearch(d.support, x)
	if i == 0 return (0.0, cdf[1])
	elseif j == n+1 return (cdf[end], 1.0)
	else return (cdf[i], cdf[i+1])
	end
end

# For the given targetAmbBusyDuration, return a dict stationsNodeSetsP, where stationsNodeSetsP[l][j] is an upper bound
# on the fraction of total demand that can be served by an ambulance that is busy for duration <= targetAmbBusyDuration,
# by assigning station j to serve node set that has stations ordered by list l (indices of nearest to furthest stations).
function calcNodeSetsAmbBusyDurationProbs(sim::Simulation, coverBound::CoverBound, targetAmbBusyDuration::Float)
	@assert(targetAmbBusyDuration >= 0)

	# shorthand:
	@unpack stations, numStations = sim
	points = sim.demandCoverage.points
	numPoints = length(points)
	nodesPoints = sim.demandCoverage.nodesPoints
	rastersPointDemands = sim.demandCoverage.rastersPointDemands

	@assert(sim.demand.numSets == 1) # otherwise would have to calculate cover bound for each demand set
	totalDemand = sum(mode -> mode.arrivalRate, sim.demand.modes) # assumes sim.demand.numSets == 1

	stationsNodeSetsP = Dict{Vector{Int}, Vector{Float}}() # stationsNodeSetsP[l][j] = fraction of total demand that can be served by an ambulance that is busy for duration <= targetAmbBusyDuration, by assigning station j to serve node set that has stationList == l
	for coverBoundMode in coverBound.modes, stationList in coverBoundMode.nodeSetsStationList
		get!(stationsNodeSetsP, stationList, zeros(Float, sim.numStations))
	end

	for coverBoundMode in coverBound.modes
		# shorthand
		@unpack stationsToPointsTimes, pointsToHospitalTimes, nodeSets, nodeSetsStationList, transportProb = coverBoundMode
		d1 = coverBoundMode.approxDistrs["d1"]
		d2 = coverBoundMode.approxDistrs["d2"]
		numNodeSets = length(nodeSets)

		# create cdf for quick look-up of probability values from approximate distributions
		cdf1 = makeCdf(d1)
		cdf2 = makeCdf(d2)

		# calculate probabilities for amb being busy for <= targetAmbBusyDuration, for each pair of station and point
		stationsPointsP = zeros(Float, numStations, numPoints) # stationsPointsP[i,j] = probability of amb busy duration (to serve point j) <= targetAmbBusyDuration from station i
		for i = 1:numStations, j = 1:numPoints
			# need to calculate probability that amb will be busy for <= targetAmbBusyDuration while serving point, based on two possible outcomes:
			# - no transport to hospital, look at probability of targetAmbBusyDuration >= durations: response travel + on-scene.
			# - transport to hospital, look at probability of targetAmbBusyDuration >= durations: response travel + on-scene + transport + at-hospital.
			# calculate the above probabilities and multiply by the probability of that outcome, sum for total probability of amb being busy for time <= targetAmbBusyDuration.
			t1 = stationsToPointsTimes[i,j]
			t2 = t1 + pointsToHospitalTimes[j]
			p1 = binarySearch(d1, cdf1, targetAmbBusyDuration - t1)[2] # probability of amb busy duration <= targetAmbBusyDuration, assuming no transport to hospital
			p2 = binarySearch(d2, cdf2, targetAmbBusyDuration - t2)[2] # probability of amb busy duration <= targetAmbBusyDuration, assuming transport to hospital
			stationsPointsP[i,j] = (1-transportProb) * p1 + (transportProb) * p2
		end

		# get all demand modes which give the right travel mode (based on response) for this coverBoundMode
		demandModes = DemandMode[]
		for (demandPriority, travelPriority) in sim.responseTravelPriorities
			travelMode = getTravelMode!(sim.travel, travelPriority, sim.startTime) # assumes travel.numSets == 1
			if travelMode == coverBoundMode.responseTravelMode
				demandMode = getDemandMode!(sim.demand, demandPriority, sim.startTime) # assumes demand.numSets == 1
				push!(demandModes, demandMode)
			end
		end

		pointDemandFracs = sum(mode -> rastersPointDemands[mode.rasterIndex] * mode.rasterMultiplier, demandModes) / totalDemand # pointDemandFracs[i] is fraction of demand at point i

		# group demand points
		for (i,nodeSet) in enumerate(nodeSets)
			stationList = nodeSetsStationList[i]
			nodeSetsP = stationsNodeSetsP[stationList]
			for j in nodeSet # j = node index
				for k in nodesPoints[j] # indices of points for which node j is the nearest node
					for s = 1:numStations
						nodeSetsP[s] += stationsPointsP[s,k] * pointDemandFracs[k]
					end
				end
			end
		end
	end

	return stationsNodeSetsP
end

# Get an upper bound on amb busy duration distribution for each
# combination of amb busy duration and number of free ambulances.
function calcAmbBusyDurationProbUpperBounds(sim, coverBound::CoverBound;
		ambBusyDurationsToSample::Vector{Float} = coverBound.ambBusyDurationsToSample,
		numAmbsList::Vector{Int} = collect(1:sim.numStations), pMedianRelax::Bool = true, doPrint::Bool = false)

	@assert(issorted(ambBusyDurationsToSample))
	@assert(issorted(numAmbsList))

	pMedianOptions = copy(pMedianDefaultOptions)
	pMedianOptions[:x_bin] = !pMedianRelax

	numTargets = length(ambBusyDurationsToSample) # shorthand
	n = length(numAmbsList) # shorthand
	ambBusyDurationProbUpperBounds = zeros(Float, numTargets, length(numAmbsList)) # ambBusyDurationProbUpperBounds[i,j] = upper bound on probability of dispatched amb being busy for duration <= ambBusyDurationsToSample[i] for j free ambulances
	for (i,t) in enumerate(ambBusyDurationsToSample)
		doPrint && print("\rTarget amb busy duration $i of $numTargets, 0 ambulances of $n ")

		stationsNodeSetsP = calcNodeSetsAmbBusyDurationProbs(sim, coverBound, t)
		stationsP = hcat(values(stationsNodeSetsP)...)

		for (j,a) in enumerate(numAmbsList)
			doPrint && print("\rTarget amb busy duration $i of $numTargets, $j ambulances of $n ")

			results = Dict()
			solvePMedian(a, -stationsP; options = pMedianOptions, results = results)
			ambBusyDurationProbUpperBounds[i,j] = -results[:cost]

			# # Testing: group node sets further, as we can assume that each station will have at least its a^th-to-last preferred station selected.
			# # Maybe only do this if pMedianRelax = false? Otherwise could lead to weaker upper bounds.
			# m = sim.numStations - a + 1
			# stationsReducedNodeSetsP = Dict{Vector{Int}, Vector{Float}}()
			# for (stationList, nodeSetsP) in stationsNodeSetsP
				# x = get!(stationsReducedNodeSetsP, stationList[1:m], zeros(Float, sim.numStations))
				# x .+= nodeSetsP
			# end
			# stationsReducedP = hcat(values(stationsReducedNodeSetsP)...)
			# results = Dict()
			# solvePMedian(a, -stationsReducedP; options = pMedianOptions, results = results)
			# ambBusyDurationProbUpperBounds[i,j] = -results[:cost]
		end
	end
	doPrint && println()

	# fix any values in ambBusyDurationProbUpperBounds that are incorrectly larger than they should be (not sure why this happens),
	# can have cases where decreasing the number of ambulances or target amb busy duration can (incorrectly) slightly increase the value from p-median problem.
	p = ambBusyDurationProbUpperBounds # shorthand
	# should have p[i,j] <= p[i+1,j] and p[i,j] <= p[i,j+1]
	m, n = size(p)
	for i = m:-1:1, j = n:-1:1
		if i < m p[i,j] = min(p[i,j], p[i+1,j]) end
		if j < n p[i,j] = min(p[i,j], p[i,j+1]) end
	end

	return ambBusyDurationProbUpperBounds
end

# create amb busy duration distribution
# requires fields of coverBound to already be populated: ambBusyDurationsToSample, ambBusyDurationProbUpperBounds
# mutates: coverBound.ambBusyDurationLowerBoundDistrs
function calcAmbBusyDurationLowerBoundDistrs!(coverBound::CoverBound)
	coverBound.ambBusyDurationLowerBoundDistrs = Sampleable[]
	ambBusyDurations = vcat(0, coverBound.ambBusyDurationsToSample)
	for j = 1:size(coverBound.ambBusyDurationProbUpperBounds,2)
		p1 = Float.(vcat(0, coverBound.ambBusyDurationProbUpperBounds[:,j], 1))
		p2 = p1[2:end] - p1[1:end-1]
		distr = DiscreteNonParametric(ambBusyDurations, p2)
		push!(coverBound.ambBusyDurationLowerBoundDistrs, sampler(distr))
	end
end

# solve MCLP to determine maximum coverage achievable with different numbers of free ambulances (1:sim.numStations)
function calcNumAmbsMaxCoverageFrac(sim::Simulation)
	@assert(sim.demand.numSets == 1)
	@assert(sim.travel.numSets == 1)

	# solve MCLP for each possible number of ambulances
	numAmbsMaxCoverageFrac = zeros(Float, sim.numStations) # numAmbsMaxCoverageFrac[i] gives maximum fraction of demand coverage (calculated from the Maximal Coverage Location Problem) achievable with i free ambulances
	# note that having more ambs than stations has no benefit for MCLP; once every station has an amb then max coverage has been achieved
	arrivalRate = sum([mode.arrivalRate for mode in sim.demand.modes]) # assumes demand.numSets == 1
	results = Dict() # for MCLP results
	for i = 1:sim.numStations
		solveMclp!(sim; numAmbs = i, results = results)
		numAmbsMaxCoverageFrac[i] = results[:objVal] / arrivalRate
	end

	return numAmbsMaxCoverageFrac
end

# for each call queued duration in queuedDurationsToSample,
# calculate the maximum coverage achievable within the remaining response time for a single ambulance
function calcQueuedDurationsMaxCoverageFrac(sim::Simulation, coverBound::CoverBound;
	queuedDurationsToSample::Vector{Float} = coverBound.queuedDurationsToSample,
	demandWeights::Dict{Priority,Float} = Dict([p => 1.0 for p in priorities]),
	doPrint::Bool = false)

	@assert(queuedDurationsToSample != [])
	@assert(issorted(queuedDurationsToSample))
	@assert(queuedDurationsToSample[1] >= 0)

	@assert(all(x -> x >= 0, values(demandWeights)))
	@assert(sum(values(demandWeights)) > 0)
	for p in keys(demandWeights)
		if demandWeights[p] == 0
			delete!(demandWeights, p)
		end
	end
	demandPriorities = collect(keys(demandWeights))

	@assert(sim.demand.numSets == 1)
	@assert(sim.travel.numSets == 1)

	numStations = sim.numStations
	numPoints = length(sim.demandCoverage.points)
	arrivalRate = sum(mode -> mode.arrivalRate, sim.demand.modes) # assumes demand.numSets == 1

	# get data for each demand priority
	# cannot necessarily combine these (e.g. combine by common travel modes), as different demand priorities may have different target response times
	demandPrioritiesData = Dict()
	currentTime = sim.startTime
	for demandPriority in demandPriorities
		responseTravelPriority = sim.responseTravelPriorities[demandPriority]
		demandMode = getDemandMode!(sim.demand, demandPriority, currentTime) # assumes demand.numSets == 1
		travelMode = getTravelMode!(sim.travel, responseTravelPriority, sim.startTime) # assumes sim.travel.numSets == 1
		demandPriorityData = Dict()
		demandPriorityData[:coverBoundMode] = coverBound.modes[coverBound.modeLookup[travelMode.index]]
		demandPriorityData[:rasterPointDemands] = sim.demandCoverage.rastersPointDemands[demandMode.rasterIndex] * demandMode.rasterMultiplier * demandWeights[demandPriority]
		demandPrioritiesData[demandPriority] = demandPriorityData
	end

	maxCoverageFracs = zeros(Float, length(queuedDurationsToSample)) # maxCoverageFracs[i] = mclp value for queuedDurationsToSample[i]
	# stationsNumAmbsList = []
	for (i,queuedDuration) in enumerate(queuedDurationsToSample)
		doPrint && print("\rQueued duration $i of $(length(queuedDurationsToSample))")

		# calculate demand covered by different station sets
		pointData = Dict{Vector{Bool}, Float}() # pointData[stationsCoverPoint] = pointDemand
		for demandPriority in demandPriorities
			coverTime = sim.targetResponseDurations[Int(demandPriority)] - (coverBound.dispatchDelay + queuedDuration)
			if coverTime > 0
				demandPriorityData = demandPrioritiesData[demandPriority]
				rasterPointDemands = demandPriorityData[:rasterPointDemands]
				stationsCoverPoints = demandPriorityData[:coverBoundMode].stationsToPointsTimes .<= coverTime # would need to change this if dispatch delay is not deterministic
				for j = 1:numPoints
					stationsCoverPoint = stationsCoverPoints[:,j] # stationSet[i] = true if station i covers point j
					get!(pointData, stationsCoverPoint, 0.0)
					pointData[stationsCoverPoint] += rasterPointDemands[j]
				end
			end
		end
		pointStations = findall.(collect(keys(pointData)))
		pointDemands = collect(values(pointData))
		# point j has demand pointDemands[j] and is covered by stations pointStations[j]

		if isempty(pointStations) || (length(pointStations) == 1 && !any(pointStations[1]))
			break # cannot reach any points in time
		end

		# solve mclp
		results = Dict()
		stationsNumAmbs = solveMclp(1, pointDemands, pointStations, results = results)
		maxCoverageFracs[i] = results[:objVal] / arrivalRate
		# push!(stationsNumAmbsList, stationsNumAmbs)
	end

	# # if call was queued for no time, coverage should match that of MCLP for one amb and full target response time (minus dispatch delay)
	# if queuedDurationsToSample[1] == 0 @assert(isapprox(maxCoverageFracs[1], coverBound.numAmbsMaxCoverageFrac[1])) end

	return maxCoverageFracs # max probability of covering call that has been queued for queuedDurationsToSample[i] is maxCoverageFracs[i]
end

function initCoverBoundSim(; numReps::Int = 1, numAmbs::Int = 0, warmUpDuration::Float = 0.0, minLastCallArrivalTime::Float = nullTime,
		interarrivalTimeDistrRng::Union{DistrRng,Nothing} = nothing,
		arrivalRate::Float = 0.0, interarrivalTimeSeed::Int = 0, # use these if interarrivalTimeDistrRng == nothing
		ambBusyDurationSeed::Int = 1)

	@assert(numReps >= 1)
	@assert(numAmbs >= 1)
	@assert(warmUpDuration >= 0)
	@assert(minLastCallArrivalTime > warmUpDuration)

	if interarrivalTimeDistrRng == nothing
		@assert(arrivalRate > 0)
		interarrivalTimeDistrRng = DistrRng(Exponential(1/arrivalRate), seed = interarrivalTimeSeed)
	else
		@assert(isa(interarrivalTimeDistrRng.d, Exponential))
		@assert(mean(interarrivalTimeDistrRng.d) != Inf)
	end

	ambBusyDurationRng = MersenneTwister(ambBusyDurationSeed)

	coverBoundSim = CoverBoundSim()
	@pack! coverBoundSim = numReps, numAmbs, warmUpDuration, minLastCallArrivalTime, interarrivalTimeDistrRng, ambBusyDurationRng

	return coverBoundSim
end

# mutates: coverBoundSim.interarrivalTimeDistrRng, coverBoundSim.ambBusyDurationRng
function simulateRepCoverBound!(coverBound::CoverBound)
	cbs = coverBound.sim # shorthand
	@assert(cbs.warmUpDuration > 0)
	@assert(cbs.interarrivalTimeDistrRng != nothing)
	@assert(cbs.numAmbs >= 1)
	@assert(0 <= cbs.warmUpDuration < cbs.minLastCallArrivalTime)

	numStations = length(coverBound.ambBusyDurationLowerBoundDistrs)
	@assert(numStations >= 1)

	# generate call arrival times
	# could do this within sim loop, but will do here for simplicity
	callArrivalTimes = Float[]
	t = 0.0
	while t < cbs.minLastCallArrivalTime
		t += rand(cbs.interarrivalTimeDistrRng)
		push!(callArrivalTimes, t)
	end
	numCalls = length(callArrivalTimes)
	@assert(numCalls >= 1)

	# run simple sim, count number of free ambulances at moment of each dispatch
	currentTime = 0.0
	eventList = Event[]
	addEvent!(eventList, form = callArrives, time = callArrivalTimes[1])
	numFreeAmbs = cbs.numAmbs
	numFreeAmbsCount = zeros(Int, cbs.numAmbs) # numFreeAmbsCount[i] = how many times calls arrive when there are i free ambulances
	numQueuedCalls = 0
	earliestQueuedCallIndex = nullIndex # of the currently queued calls, index of the call that arrived earliest
	queuedCallDurations = Float[] # duration for which calls (past and present) were queued
	callIndex = 1
	while !isempty(eventList)
		event = getNextEvent!(eventList)
		currentTime = event.time
		if event.form == callArrives
			event = addEvent!(eventList, form = considerDispatch, time = callArrivalTimes[callIndex] + coverBound.dispatchDelay)
			event.callIndex = callIndex
			callIndex += 1
			if callIndex <= numCalls
				event = addEvent!(eventList, form = callArrives, time = callArrivalTimes[callIndex])
			end

		elseif event.form == considerDispatch
			if numFreeAmbs == 0
				numQueuedCalls += 1
				if earliestQueuedCallIndex == nullIndex
					earliestQueuedCallIndex = event.callIndex
				end
			else
				addEvent!(eventList; form = ambDispatched, time = currentTime)
			end

		elseif event.form == ambDispatched
			recordStats = event.time >= cbs.warmUpDuration
			if earliestQueuedCallIndex == nullIndex # no queued calls
				recordStats && (numFreeAmbsCount[numFreeAmbs] += 1)
			else
				# assume that queued calls are responded to in first-in first-out order
				queuedDuration = event.time - (callArrivalTimes[earliestQueuedCallIndex] + coverBound.dispatchDelay)
				@assert(queuedDuration >= 0)
				recordStats && push!(queuedCallDurations, queuedDuration)
				earliestQueuedCallIndex = numQueuedCalls == 0 ? nullIndex : earliestQueuedCallIndex + 1
			end

			# generate amb busy duration
			t = rand(cbs.ambBusyDurationRng, coverBound.ambBusyDurationLowerBoundDistrs[min(numFreeAmbs, numStations)])
			addEvent!(eventList; form = ambBecomesFree, time = currentTime + t)

			numFreeAmbs -= 1
		elseif event.form == ambBecomesFree
			numFreeAmbs += 1
			if numQueuedCalls >= 1
				addEvent!(eventList; form = ambDispatched, time = currentTime)
				numQueuedCalls -= 1
			end
		else error("Unrecognised event.")
		end
	end

	numCallsRecorded = sum(numFreeAmbsCount) + length(queuedCallDurations)
	return numCallsRecorded, numFreeAmbsCount, queuedCallDurations
end

# mutates: coverBound.sim.reps, coverBound.sim.bound
function simulateCoverBound!(coverBound::CoverBound)
	coverBound.sim.reps = [CoverBoundSimRep() for i = 1:coverBound.sim.numReps]
	for rep in coverBound.sim.reps
		rep.numCalls, rep.numFreeAmbsCount, rep.queuedCallDurations = simulateRepCoverBound!(coverBound)
	end
	calcCoverBound!(coverBound)
	return coverBound.sim.bound
end

# calculate cover bound value after simulating
# mutates: coverBound.sim.bound, rep.bound for rep in coverBound.sim.reps
function calcCoverBound!(coverBound::CoverBound; accountForQueuedDurations::Bool = coverBound.accountForQueuedDurations)
	cb = coverBound # shorthand
	for (i,rep) in enumerate(cb.sim.reps)
		@unpack numCalls, numFreeAmbsCount, queuedCallDurations = rep
		n = min(length(cb.numAmbsMaxCoverageFrac), length(numFreeAmbsCount))
		coverageFracTotal = sum(cb.numAmbsMaxCoverageFrac[1:n] .* numFreeAmbsCount[1:n]) + cb.numAmbsMaxCoverageFrac[end] * sum(numFreeAmbsCount[n+1:end])
		if accountForQueuedDurations
			# account for duration that calls were queued
			for queuedCallDuration in queuedCallDurations
				j = binarySearch(cb.queuedDurationsToSample, queuedCallDuration)[1]
				if j == 0 # queuedCallDuration < queuedDurationsToSample[1], so need to assume call was queued for no time
					coverageFracTotal += cb.numAmbsMaxCoverageFrac[1] # maximum coverage from one ambulance
				else # j > 0
					coverageFracTotal += cb.queuedDurationsMaxCoverageFrac[j] # coverage accounting for queued time
				end
			end
		else
			# original cover bound, ignores duration for which calls were queued and gives these calls the maximum coverage from one ambulance
			coverageFracTotal += cb.numAmbsMaxCoverageFrac[1] * length(queuedCallDurations)
		end
		rep.bound = coverageFracTotal / numCalls
	end

	bounds = [rep.bound for rep in cb.sim.reps]
	cb.sim.bound = MeanAndHalfWidth(mean(bounds), tDistrHalfWidth(bounds))

	return cb.sim.bound
end

# Simulate to estimate bound on how much the cover bound value could be reduced if the distributions were not discrete (due to sampling) but continuous.
# This works by overestimating the ambulance busy duration and underestimating the coverage of queued calls.
# This requires setting (by knowledge or guess) the maximum duration that an ambulance can be busy (maxAmbBusyDuration).
# Note that this lower bound will change as coverBound.ambBusyDurationLowerBoundDistrs and coverBound.queuedDurationsMaxCoverageFrac are changed,
# e.g. from changing the values (durations) to sample or using the relaxed versions of sub-problems such as the p-median problem.
function simulateCoverBoundLowerBound!(coverBound::CoverBound; maxAmbBusyDuration::Float = coverBound.ambBusyDurationsToSample[end])
	cb = coverBound # shorthand

	@assert(maxAmbBusyDuration >= cb.ambBusyDurationsToSample[end])

	# create backup for some values
	coverBoundBak = CoverBound()
	fnames = (:ambBusyDurationsToSample, :ambBusyDurationProbUpperBounds, :ambBusyDurationLowerBoundDistrs, :queuedDurationsMaxCoverageFrac)
	for fname in fnames
		setfield!(coverBoundBak, fname, deepcopy(getfield(cb, fname)))
	end

	# shift ambBusyDurationProbUpperBounds right (and so change ambBusyDurationLowerBoundDistrs)
	n = size(cb.ambBusyDurationProbUpperBounds, 2)
	cb.ambBusyDurationProbUpperBounds = vcat(zeros(1,n), cb.ambBusyDurationProbUpperBounds)
	if maxAmbBusyDuration > cb.ambBusyDurationsToSample[end]
		push!(cb.ambBusyDurationsToSample, maxAmbBusyDuration)
	else
		@assert(maxAmbBusyDuration == cb.ambBusyDurationsToSample[end])
		cb.ambBusyDurationProbUpperBounds = cb.ambBusyDurationProbUpperBounds[1:end-1,:]
	end
	calcAmbBusyDurationLowerBoundDistrs!(cb) # sets coverBound.ambBusyDurationLowerBoundDistrs

	# shift queuedDurationsMaxCoverageFrac left
	cb.queuedDurationsMaxCoverageFrac = vcat(cb.queuedDurationsMaxCoverageFrac[2:end], 0)

	# simulate
	simulateCoverBound!(cb)

	# reset cover bound
	for fname in fnames
		setfield!(cb, fname, getfield(coverBoundBak, fname))
	end

	return cb.sim.bound
end

## read and write files
# have added the functions to this file instead of in write_sim_files.jl and read_sim_files.jl to make it clear that they are for the cover bound

function writeAmbBusyDurationProbUpperBoundsFile(filename::String, ambBusyDurationsToSample::Vector{Float}, ambBusyDurationProbUpperBounds::Array{Float,2})
	@assert(issorted(ambBusyDurationsToSample))
	n = size(ambBusyDurationProbUpperBounds, 2) # should equal numStations
	table = Table("ambBusyDurationProbUpperBounds", vcat("duration", ["numAmbs_$i" for i = 1:n]),
		rows = [vcat(ambBusyDurationsToSample[i], ambBusyDurationProbUpperBounds[i,:]) for i = 1:length(ambBusyDurationsToSample)])
	writeTablesToFile(filename, table)
end
function writeAmbBusyDurationProbUpperBoundsFile(filename::String, coverBound::CoverBound)
	writeAmbBusyDurationProbUpperBoundsFile(filename, coverBound.ambBusyDurationsToSample, coverBound.ambBusyDurationProbUpperBounds)
end

function readAmbBusyDurationProbUpperBoundsFile(filename::String)
	table = readTablesFromFile(filename)["ambBusyDurationProbUpperBounds"]
	ambBusyDurationsToSample = haskey(table.columns, "duration") ? table.columns["duration"] : table.columns["time"] # compat (time)
	ambBusyDurationProbUpperBounds = table.data[:,2:end]
	return ambBusyDurationsToSample, ambBusyDurationProbUpperBounds
end

function writeQueuedDurationsMaxCoverageFracFile(filename::String, queuedDurationsToSample::Vector{Float}, queuedDurationsMaxCoverageFrac::Vector{Float})
	table = Table("queuedDurationsMaxCoverageFrac", ["duration", "maxCoverageFrac"];
		rows = [[queuedDurationsToSample[i], queuedDurationsMaxCoverageFrac[i]] for i = 1:length(queuedDurationsToSample)])
	writeTablesToFile(filename, table)
end
function writeQueuedDurationsMaxCoverageFracFile(filename::String, coverBound::CoverBound)
	writeQueuedDurationsMaxCoverageFracFile(filename, coverBound.queuedDurationsToSample, coverBound.queuedDurationsMaxCoverageFrac)
end

function readQueuedDurationsMaxCoverageFracFile(filename::String)
	table = readTablesFromFile(filename)["queuedDurationsMaxCoverageFrac"]
	queuedDurationsToSample = haskey(table.columns, "duration") ? table.columns["duration"] : table.columns["queuedDuration"] # compat (queuedDuration)
	queuedDurationsMaxCoverageFrac = table.columns["maxCoverageFrac"]
	return queuedDurationsToSample, queuedDurationsMaxCoverageFrac
end


text:
"""
    neighbors(vec::AbstractVector)

Returns an iterator over neighboring pairs.
```julia
julia> neighbors([1,2,3]) |> collect
2-element Array{Tuple{Int64,Int64},1}:
 (1, 2)
 (2, 3)
 ```
 """
function neighbors(itr::AbstractVector)
    first  = @view itr[1:end-1]
    second = @view itr[2:end]
    zip(first, second)
end

ext: jl


text:
using Test
using BEncode

# Decode

@test bdecode("i42e") == 42
@test bdecode("i0e") == 0
@test bdecode("i-42e") == -42

@test bdecode("6:Trajan") == "Trajan"
@test bdecode("0:") == ""

@test bdecode("li4ei2ee") == [4, 2]
@test bdecode("li4e1:2e") == [4, "2"]
@test bdecode("le") == []

@test bdecode("d4:dictli1ei4eee") == Dict("dict" => [1, 4])
@test bdecode("de") == Dict()

# Encode

@test bencode(42) == "i42e"
@test bencode(0) == "i0e"
@test bencode(-42) == "i-42e"

@test bencode("Trajan") == "6:Trajan"
@test bencode("") == "0:"

@test bencode([4, 2]) == "li4ei2ee"
@test bencode([4, "2"]) == "li4e1:2e"
@test bencode([]) == "le"

@test bencode(Dict("dict" => [1, 4])) == "d4:dictli1ei4eee"
@test bencode(Dict()) == "de"


text:
module NDTensors

using Random
using LinearAlgebra
using StaticArrays
using HDF5
using Strided

#####################################
# Exports
#
include("exports.jl")

#####################################
# DenseTensor and DiagTensor
#
include("tupletools.jl")
include("dims.jl")
include("tensorstorage.jl")
include("tensor.jl")
include("contraction_logic.jl")
include("dense.jl")
include("symmetric.jl")
include("linearalgebra.jl")
include("diag.jl")
include("combiner.jl")
include("truncate.jl")
include("svd.jl")

#####################################
# BlockSparseTensor
#
include("blocksparse/blockdims.jl")
include("blocksparse/blockoffsets.jl")
include("blocksparse/blocksparse.jl")
include("blocksparse/blocksparsetensor.jl")
include("blocksparse/diagblocksparse.jl")
include("blocksparse/combiner.jl")
include("blocksparse/linearalgebra.jl")

#####################################
# Empty
#
include("empty.jl")

end # module NDTensors


text:
module UI

import GLFW
using ModernGL

include("gltools.jl")

GLFW.Init()

# OS X-specific GLFW hints to initialize the correct version of OpenGL
@osx_only begin
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
end

# Create a windowed mode window and its OpenGL context
window = GLFW.CreateWindow(600, 600, "OpenGL Example")

# Make the window's context current
GLFW.MakeContextCurrent(window)

# The data for our triangle
data = GLfloat[
    0.0, 0.5,
    0.5, -0.5,
    -0.5,-0.5
]

# Generate a vertex array and array buffer for our data
vao = glGenVertexArray()
glBindVertexArray(vao)

vbo = glGenBuffer()
glBindBuffer(GL_ARRAY_BUFFER, vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(data), data, GL_STATIC_DRAW)

# Create and initialize shaders
const vsh = """
    #version 330
    in vec2 position;

    void main() {
        gl_Position = vec4(position, 0.0, 1.0);
    }
"""

const fsh = """
    #version 330
    out vec4 outColor;

    void main() {
        outColor = vec4(1.0, 1.0, 1.0, 1.0);
    }
"""

vertexShader = createShader(vsh, GL_VERTEX_SHADER)
fragmentShader = createShader(fsh, GL_FRAGMENT_SHADER)
program = createShaderProgram(vertexShader, fragmentShader)
glUseProgram(program)

positionAttribute = glGetAttribLocation(program, "position");

glEnableVertexAttribArray(positionAttribute)
glVertexAttribPointer(positionAttribute, 2, GL_FLOAT, false, 0, 0)

t = 0

# Loop until the user closes the window
while !GLFW.WindowShouldClose(window)
    # Pulse the background blue
    t += 1
    glClearColor(0.0, 0.0, 0.5 * (1 + sin(t * 0.02)), 1.0)
    glClear(GL_COLOR_BUFFER_BIT)
    # Draw our triangle
    glDrawArrays(GL_TRIANGLES, 0, 3)

    # Swap front and back buffers
    GLFW.SwapBuffers(window)

    # Poll for and process events
    GLFW.PollEvents()
end

GLFW.Terminate()

end

max_stars_repo_name: yurivish/Julia-Playground

max_issues_repo_name: yurivish/Julia-Playground

max_forks_repo_name: yurivish/Julia-Playground


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## User Options

mutable struct Options
    hascolor::Bool
    extra_keymap::Union{Dict,Vector{<:Dict}}
    # controls the presumed tab width of code pasted into the REPL.
    # Must satisfy `0 < tabwidth <= 16`.
    tabwidth::Int
    # Maximum number of entries in the kill ring queue.
    # Beyond this number, oldest entries are discarded first.
    kill_ring_max::Int
    region_animation_duration::Float64
    beep_duration::Float64
    beep_blink::Float64
    beep_maxduration::Float64
    beep_colors::Vector{String}
    beep_use_current::Bool
    backspace_align::Bool
    backspace_adjust::Bool
    confirm_exit::Bool # ^D must be repeated to confirm exit
    auto_indent::Bool # indent a newline like line above
    auto_indent_tmp_off::Bool # switch auto_indent temporarily off if copy&paste
    auto_indent_bracketed_paste::Bool # set to true if terminal knows paste mode
    # cancel auto-indent when next character is entered within this time frame :
    auto_indent_time_threshold::Float64
    # refresh after time delay
    auto_refresh_time_delay::Float64
    # default IOContext settings at the REPL
    iocontext::Dict{Symbol,Any}
end

Options(;
        hascolor = true,
        extra_keymap = AnyDict[],
        tabwidth = 8,
        kill_ring_max = 100,
        region_animation_duration = 0.2,
        beep_duration = 0.2, beep_blink = 0.2, beep_maxduration = 1.0,
        beep_colors = ["\e[90m"], # gray (text_colors not yet available)
        beep_use_current = true,
        backspace_align = true, backspace_adjust = backspace_align,
        confirm_exit = false,
        auto_indent = true,
        auto_indent_tmp_off = false,
        auto_indent_bracketed_paste = false,
        auto_indent_time_threshold = 0.005,
        auto_refresh_time_delay = 0.05,
        iocontext = Dict{Symbol,Any}()) =
            Options(hascolor, extra_keymap, tabwidth,
                    kill_ring_max, region_animation_duration,
                    beep_duration, beep_blink, beep_maxduration,
                    beep_colors, beep_use_current,
                    backspace_align, backspace_adjust, confirm_exit,
                    auto_indent, auto_indent_tmp_off, auto_indent_bracketed_paste,
                    auto_indent_time_threshold, auto_refresh_time_delay,
                    iocontext)

# for use by REPLs not having an options field
const GlobalOptions = Options()


text:
using Test
include("../../utils/distancias.jl")

@testset "Distancias" begin
    @test EuclideanDistance([1,1,0],[2,1,2]) == sqrt(5)
    @test WeightedEuclideanDistance([1,1,1],[1,1,0],[5,5,4]) == 2.0
end


text:
# # Bacterial Growth

# ![](bacteria.gif)

# Bacterial colonies are a prime example for growing active matter, where systems
# are driven out of equilibrium by proliferation.
# This model is a simplified version of unpublished work by Yoav G. Pollack and
# [Philip Bittihn](https://www.ds.mpg.de/lmp/dyn-sys); similar models can be found in literature.
# Here, a bacterium is modelled by two soft disk "nodes" linked by a spring,
# whose rest length grows with a constant growth rate. When it has reached its
# full extension, the cell divides into two daughter cells with the same orientation.

# This example is a showcase of a complex continuous system.
# Agents will be splitting into more agents, thus having agent generation in continuous
# space. The model also uses advanced agent movement in continuous space, where a
# specialized "`move_agent`" function is created. Advanced plotting is also done,
# since each agent is a specialized shape.
# It is also available from the `Models` module as [`Models.growing_bacteria`](@ref).

using Agents, LinearAlgebra
using Random # hide
cd(@__DIR__) #src

mutable struct SimpleCell <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    length::Float64
    orientation::Float64
    growthprog::Float64
    growthrate::Float64

    ## node positions/forces
    p1::NTuple{2,Float64}
    p2::NTuple{2,Float64}
    f1::NTuple{2,Float64}
    f2::NTuple{2,Float64}
end

function SimpleCell(id, pos, l, φ, g, γ)
    a = SimpleCell(id, pos, l, φ, g, γ, (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0))
    update_nodes!(a)
    return a
end

# In this model, the agents have to store their state in two redundant ways:
# the cell coordinates (position, length, orientation) are required for the
# equations of motion, while the positions of the disk-shaped nodes are necessary
# for calculating mechanical forces between cells. To transform from one set of
# coordinates to the other, we need to write a function
function update_nodes!(a::SimpleCell)
    offset = 0.5 * a.length .* unitvector(a.orientation)
    a.p1 = a.pos .+ offset
    a.p2 = a.pos .- offset
end
nothing # hide

# Some geometry convenience functions
unitvector(φ) = reverse(sincos(φ))
cross2D(a, b) = a[1] * b[2] - a[2] * b[1]
nothing # hide

# ## Stepping functions

function model_step!(model)
    for a in allagents(model)
        if a.growthprog ≥ 1
            ## When a cell has matured, it divides into two daughter cells on the
            ## positions of its nodes.
            add_agent!(a.p1, model, 0.0, a.orientation, 0.0, 0.1 * rand() + 0.05)
            add_agent!(a.p2, model, 0.0, a.orientation, 0.0, 0.1 * rand() + 0.05)
            kill_agent!(a, model)
        else
            ## The rest lengh of the internal spring grows with time. This causes
            ## the nodes to physically separate.
            uv = unitvector(a.orientation)
            internalforce = model.hardness * (a.length - a.growthprog) .* uv
            a.f1 = -1 .* internalforce
            a.f2 = internalforce
        end
    end
    ## Bacteria can interact with more than on other cell at the same time, therefore,
    ## we need to specify the option `:all` in `interacting_pairs`
    for (a1, a2) in interacting_pairs(model, 2.0, :all)
        interact!(a1, a2, model)
    end
end
nothing # hide

# Here we use a custom [`move_agent!`](@ref) function,
# because the agents have several moving parts.
# Notice that the first derivatives of all degrees of freedom is directly
# proportional to the force applied to them. This overdamped approximation is
# valid for small length scales, where viscous forces dominate over inertia.
function agent_step!(agent::SimpleCell, model::ABM)
    fsym, compression, torque = transform_forces(agent)
    new_pos = agent.pos .+ model.dt * model.mobility .* fsym
    move_agent!(agent, new_pos, model)
    agent.length += model.dt * model.mobility .* compression
    agent.orientation += model.dt * model.mobility .* torque
    agent.growthprog += model.dt * agent.growthrate
    update_nodes!(agent)
    return agent.pos
end
nothing # hide

# ### Helper functions
function interact!(a1::SimpleCell, a2::SimpleCell, model)
    n11 = noderepulsion(a1.p1, a2.p1, model)
    n12 = noderepulsion(a1.p1, a2.p2, model)
    n21 = noderepulsion(a1.p2, a2.p1, model)
    n22 = noderepulsion(a1.p2, a2.p2, model)
    a1.f1 = @. a1.f1 + (n11 + n12)
    a1.f2 = @. a1.f2 + (n21 + n22)
    a2.f1 = @. a2.f1 - (n11 + n21)
    a2.f2 = @. a2.f2 - (n12 + n22)
end

function noderepulsion(p1::NTuple{2,Float64}, p2::NTuple{2,Float64}, model::ABM)
    delta = p1 .- p2
    distance = norm(delta)
    if distance ≤ 1
        uv = delta ./ distance
        return (model.hardness * (1 - distance)) .* uv
    end
    return (0, 0)
end

function transform_forces(agent::SimpleCell)
    ## symmetric forces (CM movement)
    fsym = agent.f1 .+ agent.f2
    ## antisymmetric forces (compression, torque)
    fasym = agent.f1 .- agent.f2
    uv = unitvector(agent.orientation)
    compression = dot(uv, fasym)
    torque = 0.5 * cross2D(uv, fasym)
    return fsym, compression, torque
end
nothing # hide

# ## Animating bacterial growth

# Okay, we can now initialize a model and see what it does.

space = ContinuousSpace((14, 9), 1.0; periodic = false)
model = ABM(
    SimpleCell,
    space,
    properties = Dict(:dt => 0.005, :hardness => 1e2, :mobility => 1.0),
)

# Let's start with just two agents.

add_agent!((6.5, 4.0), model, 0.0, 0.3, 0.0, 0.1)
add_agent!((7.5, 4.0), model, 0.0, 0.0, 0.0, 0.1)
nothing # hide

# The model has several parameters, and some of them are of interest.
# We could e.g. define

adata = [:pos, :length, :orientation, :growthprog, :p1, :p2, :f1, :f2]
nothing # hide

# and then [`run!`](@ref) the model. But we'll animate the model directly.

# Here we once again use the huge flexibility provided by [`plotabm`](@ref) to
# plot the becteria cells. We define a function that creates a custom `Shape` based
# on the agent:
using AgentsPlots
gr() # hide

function cassini_oval(agent)
    t = LinRange(0, 2π, 50)
    a = agent.growthprog
    b = 1
    m = @. 2 * sqrt((b^4 - a^4) + a^4 * cos(2 * t)^2) + 2 * a^2 * cos(2 * t)
    C = sqrt.(m / 2)

    x = C .* cos.(t)
    y = C .* sin.(t)

    uv = reverse(sincos(agent.orientation))
    θ = atan(uv[2], uv[1])
    R = [cos(θ) -sin(θ); sin(θ) cos(θ)]

    bacteria = R * permutedims([x y])
    Shape(bacteria[1, :], bacteria[2, :])
end
nothing # hide

# set up some nice colors

bacteria_colors(agent) =
    HSV.(agent.id * 2.718 .% 1, agent.id * 3.14 .% 1, agent.id * 1.618 .% 1)
nothing # hide

# and proceed with the animation

Random.seed!(1680)
e = model.space.extent
anim = @animate for i in 0:50:5000
    step!(model, agent_step!, model_step!, 100)
    p1 = plotabm(
        model,
        am = cassini_oval,
        as = 30,
        ac = bacteria_colors,
        showaxis = false,
        grid = false,
        xlims = (0, e[1]),
        ylims = (0, e[2]),
    )
    title!(p1, "n = $(i)")
end

gif(anim, "bacteria.gif", fps = 25)


text:
using Dash

app = dash()

app.layout = html_div() do
    dcc_input(id="input"),
    html_div(id="output-clientside"),
    html_div(id="output-serverside")
end

callback!(app, Output("output-serverside", "children"), Input("input","value")) do value
    return "Server says \"$(value)\""
end

callback!(
    ClientsideFunction("clientside", "display"),
    app,
    Output("output-clientside", "children"),
    Input("input", "value")
)

run_server(app)


text:
using Test, LinearAlgebra
using DiffEqOperators
using ModelingToolkit

@testset "utility functions" begin
    @test DiffEqOperators.unit_indices(2) == (CartesianIndex(1,0), CartesianIndex(0,1))
    @test DiffEqOperators.add_dims(zeros(2,2), ndims(zeros(2,2)) + 2) == [6. 6.; 0. 0.; 0. 0.]
    @test DiffEqOperators.perpindex(collect(1:5), 3) == [1, 2, 4, 5]
    @test DiffEqOperators.perpsize(zeros(2,2,3,2), 3) == (2, 2, 2)
end

@testset "count differentials 1D" begin
    @parameters t x
    @variables u(..)
    Dt = Differential(t)

    Dx = Differential(x)
    eq  = Dt(u(t,x)) ~ -Dx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))

    Dxx = Differential(x)^2
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))

    Dxxxx = Differential(x)^4
    eq  = Dt(u(t,x)) ~ -Dxxxx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 4
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
end

@testset "count differentials 2D" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dt = Differential(t)

    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
    @test first(DiffEqOperators.differential_order(eq.rhs, y.val)) == 2
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
    @test isempty(DiffEqOperators.differential_order(eq.lhs, y.val))
end

@testset "count with mixed terms" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dx = Differential(x)
    Dy = Differential(y)
    Dt = Differential(t)

    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + Dx(Dy(u(t,x,y)))
    @test DiffEqOperators.differential_order(eq.rhs, x.val) == Set([2, 1])
    @test DiffEqOperators.differential_order(eq.rhs, y.val) == Set([2, 1])
end

@testset "Kuramoto–Sivashinsky equation" begin
    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dx2 = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    α = 1
    β = 4
    γ = 1
    eq = Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0
    @test DiffEqOperators.differential_order(eq.lhs, x.val) == Set([4, 3, 2, 1])
end


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

__precompile__()

"""
This package implements the evaluation of the cumulative distribution function of the `Squares test statistic` originally defined in

Frederik Beaujean and Allen Caldwell. *A Test Statistic for Weighted Runs.*

https://arxiv.org/abs/1005.3233

The authors further derived an approximation to be able to compute the cumulative also for large numbers of observations in

Frederik Beaujean and Allen Caldwell. *Is the bump significant? An axion-search example*

https://arxiv.org/abs/1710.06642

Where they renamed the weighted-runs statistic to the `*Squares statistic*`.
This code is based on the original implementation by Frederik Beaujean in c++ and mathematica:

https://github.com/fredRos/runs
"""
module RunStatistics

using Distributions
using QuadGK
using ArgCheck
using Distributions
# using Threads

include("t_obs.jl")
include("partitions.jl")
include("squares.jl")
include("squares_approx.jl")

end # module


text:
# Rounded random numbers from discrete uniform distribution between [1,n]
unidrnd(n::Int64) = round(Int,rand(Uniform(1,n)))
# Sort rows based on the col reference
function sortrows!(A, col; order=false)
    if isa(col, Int)
        return A[sortperm(A[:,col], rev=order),:]
    else
        # Sorted along col[1]
        A_sorted = A[sortperm(A[:,col[1]], rev = order),:]
        # Break ties
    end
end
# Check convergence
isconverged!(new::Float64, old::Float64, reltol::Float64, iteration::Int64) = abs(new - old) / new <= reltol ? iteration += 1 : iteration = 0
# Check time
istimedout(t::Float64, t0::Float64, time_limit::Int64) = (t - t0) / 3600. > time_limit
# Status
function get_status(it::Int64, it_unchanged::Int64, timedout::Bool, options::Options)
    if it >= options.iterations
        status = "+ Status: Maximum iteration number is reached"
    elseif !timedout
        status = "+ Status: Timed out"
    elseif it_unchanged >= options.itmax_unchanged
        status = "+ Status: Converged"
    end
    return status
end
# Show results in REPL
function show_results(results::MetaheuristicResults)
    if results.options.log
        # Summary
        println("___")
        println()
        println("Optimization summary...")
        println("___")
        println()
        # Status
        println(results.status)
        # Minimizer
        println("+ Optimal decisions: ", results.minimizer)
        # Minimum
        println("+ Minimum: ", results.minimum)
        # Iteration
        println("+ Iterations: ", results.iterations)
    end
end
# Show verbose in REPL
function show_verbose(it::Int64, time::Float64, obj::Float64)
    # Header
    if it == 1
        println()
        println("___")
        println()
        println("Iteration     Objective        Time (s)")
        # Write to log
        open("metaheuristic_log.txt", "a+") do io
            write(io, "Iteration;Objective;Time (s) \n")
        end
    end
    # Iteration
    print(it)
    # Objective
    print("             ", round(obj, digits = 2))
    # Time
    print("        ", round(time, digits = 1))
    println()
    # Write to log
    open("metaheuristic_log.txt", "a+") do io
        write(io, string(it, ";", round(obj, digits = 2), ";", round(time, digits = 1), "\n"))
    end
end


text:
# This file is part of Kpax3. License is MIT.

"""
logp(C | R)
"""
function logpriorC(C::Matrix{UInt8},
                   cl::Vector{Int},
                   k::Int,
                   priorC::AminoAcidPriorCol)
  logp = 0.0

  for b in 1:size(C, 2)
    if C[cl[1], b] == UInt8(1)
      logp += priorC.logγ[1]
    elseif C[cl[1], b] == UInt8(2)
      logp += priorC.logγ[2]
    else
      logp += priorC.logγ[3]

      for l in 1:k
        if C[cl[l], b] == UInt8(3)
          logp += priorC.logω[k][1]
        else
          logp += priorC.logω[k][2]
        end
      end
    end
  end

  logp
end

function logpriorC(C::Matrix{UInt8},
                   k::Int,
                   priorC::AminoAcidPriorCol)
  logp = 0.0

  for b in 1:size(C, 2)
    if C[1, b] == UInt8(1)
      logp += priorC.logγ[1]
    elseif C[1, b] == UInt8(2)
      logp += priorC.logγ[2]
    else
      logp += priorC.logγ[3]

      for l in 1:k
        if C[l, b] == UInt8(3)
          logp += priorC.logω[k][1]
        else
          logp += priorC.logω[k][2]
        end
      end
    end
  end

  logp
end

"""
logp(S | R, X)
"""
function logcondpostS(S::Vector{UInt8},
                      cl::Vector{Int},
                      k::Int,
                      v::Vector{Int},
                      n1s::Matrix{Float64},
                      priorC::AminoAcidPriorCol)
  logp = 0.0

  lγ = zeros(Float64, 3)

  g = 0
  tmp1 = zeros(Float64, 2)
  tmp2 = 0.0

  for b in 1:length(S)
    lγ[1] = priorC.logγ[1]
    lγ[2] = priorC.logγ[2]
    lγ[3] = priorC.logγ[3]

    for l in 1:k
      g = cl[l]

      lγ[1] += logmarglik(n1s[g, b], v[g], priorC.A[1, b], priorC.B[1, b])
      lγ[2] += logmarglik(n1s[g, b], v[g], priorC.A[2, b], priorC.B[2, b])

      tmp1[1] = priorC.logω[k][1] +
                logmarglik(n1s[g, b], v[g], priorC.A[3, b], priorC.B[3, b])

      tmp1[2] = priorC.logω[k][2] +
                logmarglik(n1s[g, b], v[g], priorC.A[4, b], priorC.B[4, b])

      if tmp1[1] > tmp1[2]
        lγ[3] += tmp1[1] + log1p(exp(tmp1[2] - tmp1[1]))
      else
        lγ[3] += tmp1[2] + log1p(exp(tmp1[1] - tmp1[2]))
      end
    end

    if (lγ[1] >= lγ[2]) && (lγ[1] >= lγ[3])
      tmp2 = lγ[1] + log1p(exp(lγ[2] - lγ[1]) + exp(lγ[3] - lγ[1]))
    elseif (lγ[2] >= lγ[1]) && (lγ[2] >= lγ[3])
      tmp2 = lγ[2] + log1p(exp(lγ[1] - lγ[2]) + exp(lγ[3] - lγ[2]))
    else
      tmp2 = lγ[3] + log1p(exp(lγ[1] - lγ[3]) + exp(lγ[2] - lγ[3]))
    end

    logp += (lγ[S[b]] - tmp2)
  end

  logp
end

"""
Logp(C | R, X)
"""
function logcondpostC(C::Matrix{UInt8},
                      cl::Vector{Int},
                      k::Int,
                      v::Vector{Int},
                      n1s::Matrix{Float64},
                      priorC::AminoAcidPriorCol)
  logp = 0.0

  lγ = zeros(Float64, 3)

  logq = zeros(Float64, 4, length(cl))

  g = 0
  tmp = 0.0

  for b in 1:size(C, 2)
    lγ[1] = priorC.logγ[1]
    lγ[2] = priorC.logγ[2]
    lγ[3] = priorC.logγ[3]

    for l in 1:k
      g = cl[l]

      logq[1, l] = logmarglik(n1s[g, b], v[g], priorC.A[1, b], priorC.B[1, b])
      logq[2, l] = logmarglik(n1s[g, b], v[g], priorC.A[2, b], priorC.B[2, b])
      logq[3, l] = priorC.logω[k][1] +
                   logmarglik(n1s[g, b], v[g], priorC.A[3, b], priorC.B[3, b])
      logq[4, l] = priorC.logω[k][2] +
                   logmarglik(n1s[g, b], v[g], priorC.A[4, b], priorC.B[4, b])

      lγ[1] += logq[1, l]
      lγ[2] += logq[2, l]

      if logq[3, l] > logq[4, l]
        tmp = log1p(exp(logq[4, l] - logq[3, l]))

        lγ[3] += logq[3, l] + tmp

        if C[g, b] == 0x03
          logp -= tmp
        elseif C[g, b] == 0x04
          logp += logq[4, l] - logq[3, l] - tmp
        end
      else
        tmp = log1p(exp(logq[3, l] - logq[4, l]))

        lγ[3] += logq[4, l] + tmp

        if C[g, b] == 0x03
          logp += logq[3, l] - logq[4, l] - tmp
        elseif C[g, b] == 0x04
          logp -= tmp
        end
      end
    end

    if (lγ[1] >= lγ[2]) && (lγ[1] >= lγ[3])
      tmp = lγ[1] + log1p(exp(lγ[2] - lγ[1]) + exp(lγ[3] - lγ[1]))
    elseif (lγ[2] >= lγ[1]) && (lγ[2] >= lγ[3])
      tmp = lγ[2] + log1p(exp(lγ[1] - lγ[2]) + exp(lγ[3] - lγ[2]))
    else
      tmp = lγ[3] + log1p(exp(lγ[1] - lγ[3]) + exp(lγ[2] - lγ[3]))
    end

    if C[cl[1], b] == 0x01
      logp += lγ[1] - tmp
    elseif C[cl[1], b] == 0x02
      logp += lγ[2] - tmp
    else
      logp += lγ[3] - tmp
    end
  end

  logp
end


text:
module Corpora

require("LoadEnvironment.jl")
require("DataTypes.jl")

## Using
using ArrayViews
using DataTypes

## Exports
export Token, DataToken, DummyToken
export Sentence
export read_linecorpus
export text
export get_context, initialize_contexts
export replace_vocab!

abstract Token;

immutable DummyToken <: Token
    word::String;
    count::Vector{CInt};
end

DummyToken(word::String) = DummyToken(word, [0])

immutable DataToken <: Token
    word;
    idx::Int64;
    count::Vector{Int64};
    vector::DenseArray{Float,1}
    gradient::DenseArray{Float,1}
end

immutable Sentence
    text::Vector{Token}
end
Sentence() = Sentence(Token[])

text(s::Vector{Token}) = join(map(x->x.word, s), " ")
text(sentence::Sentence) = text(sentence.text)

function readuntil_ws_eol(s::IO)
    out = IOBuffer()
    lineend = false
    firstchar = true
    while !eof(s)
        c = read(s, Char)
        if c == ' ' || c == '\n'
            if !firstchar
                if c == '\n'
                    lineend = true
                end
                break
            end
        else
            firstchar = false
            write(out, c)
        end
    end
    return (lineend,takebuf_string(out))
end

function read_linecorpus(path::String, init_from_data, init_from_dummy,
                            init_params;
                            dummy_vocab=Dict{String,Token}(), unk="_UNK_",
                            limit_lines=0, skiprows=0, skipunk=true,
                            no_new_vocab=false, lock_tokens=false,
                            dummies_only=false)
    # Estimate number of lines
    lcounter = 0
    open(path) do fin
        for (iline, line) in enumerate(eachline(fin))
            if iline < skiprows
                continue
            end
            lcounter += 1
            if limit_lines > 0 && lcounter >= limit_lines
                break
            end
        end
    end
    data = Array(Sentence, lcounter)
    n_words_total = 0
    lcounter = 1
    max_lensen = -1
    open(path) do fin
        iline = 1
        sentence = Sentence()
        data[lcounter] = sentence
        lensen = 0
        while !eof(fin)
            (eol, token) = readuntil_ws_eol(fin)

            if iline < skiprows
                continue
            end

            ## For every token in the text
            # Check if it is already in the vocabulary
            # If not: add it with count 1
            # If it is: increase the counters
            if !haskey(dummy_vocab, token)
                if !no_new_vocab
                    #vocab[token] = inittoken(token, lock_tokens);
                    dummy_vocab[token] = DummyToken(token);
                else
                    if skipunk
                        n_words_total += 1
                        continue
                    else
                        token = unk
                    end
                end
            end
            dummytok = dummy_vocab[token]
            push!(sentence.text,dummytok)
            lensen += 1
            dummytok.count[1] += 1
            n_words_total += 1

            if eol
                if lensen > max_lensen
                    max_lensen = lensen
                end
            end

            if eol
                if (limit_lines > 0 && lcounter >= limit_lines)
                    break
                end
                if eof(fin)
                    break
                end
                lcounter += 1
                sentence = Sentence()
                data[lcounter] = sentence
                lensen = 0
            end

        end
    end

    vocab = if dummies_only
        dummy_vocab
    else
        v = Dict{String,Token}()
        # Convert DummyTokens to real tokens
        ## Signal the model that the vocabulary has been acquired

        init_params(length(dummy_vocab))
        ## Initialize token types
        for (iword,(word, tok)) in enumerate(dummy_vocab)
            if typeof(tok) == DataToken
                v[word] = init_from_data(iword,tok)
            else
                v[word] = init_from_dummy(iword,tok)
            end
        end

        ## Replace the dummies by data tokens
        for isentence in 1:size(data,1)
            for itoken in 1:size(data[isentence].text,1)
                data[isentence].text[itoken] = v[data[isentence].text[itoken].word]
            end
        end
        v
    end

    return data, vocab, n_words_total, max_lensen;
end

function replace_vocab!(data::Vector{Sentence}, vocab::Dict{String,Token})
    for isentence in 1:size(data,1)
        for itoken in 1:size(data[isentence].text,1)
            data[isentence].text[itoken] = vocab[data[isentence].text[itoken].word]
        end
    end
end

function get_context(sentence::Vector{Token}, pos::Int64, windowsize::Int64)
    lensen = size(sentence,1)
    windowsize_half = int64(windowsize/2)
    window_start = max(1, pos - windowsize_half)
    window_end = min(lensen, pos + windowsize_half)
    window = Token[]
    for i=window_start:window_end
        if i != pos
            push!(window, sentence[i])
        end
    end
    return window
end

function initialize_contexts(data::Vector{Sentence}, windowsize::Int64,
                                minlensen=5)
    for sentence in data
        if size(sentence.text,1) < 5
            continue
        end
        for pos=1:size(sentence.text,1)
            empty!(sentence.text[pos].context)
            append!(sentence.text[pos].context,
                    get_context(sentence.text, pos, windowsize))
        end
    end
end

end

max_stars_repo_name: ogh/binclusion

max_issues_repo_name: ogh/binclusion

max_forks_repo_name: ogh/binclusion


text:
ImPlotPoint() = ImPlotPoint(0, 0)
ImPlotPoint(p::ImVec2) = ImPlotPoint(p.x, p.y)
ImPlotRange() = ImPlotRange(0, 0)
ImPlotLimits() = ImPlotLimits(ImPlotRange(), ImPlotRange())


text:
"""
    disp(x)
    disp(x...)
    disp(io::IO,x)
    disp(io::IO,x...)
    disp(m::Array; format)
    disp(io::IO, m::Array; format)

Write an object `x` (or objects `x...`) to `io` (or `stdout` if `io` is not given).
The output is the same as `println` and `Base.print_matrix`, but the element type will not be shown.
If an array is given, you can specify the C-style format with `format`.

```juliadoctests
A = [1 2; 3 4]
b = [9,10]

disp(A)
disp(b)
disp(A,b)
```
"""
function disp(x)
   println(x)
end
function disp(args...)
   for x in args
      disp(x)
   end
end
function disp(io::IO,x)
   println(io,x)
end
function disp(io::IO,args...)
   for x in args
      disp(io,x)
   end
end
function disp(x::Array; format::Union{Nothing,String}=nothing)
   disp(stdout,x, format=format)
end
function disp(io::IO,x::Array; format::Union{Nothing,String}=nothing)
   if isnothing(format) || ndims(x)>2
      Base.print_matrix(io,x)
      print(io,"\n")
   elseif ndims(x)==1
      for i=1:length(x)
         print_formatted(io,format,x[i])
         print("\n")
      end
   elseif ndims(x)==2
      for i=1:size(x,1)
         for j=1:size(x,2)
            print_formatted(io,format,x[i,j])
         end
         print("\n")
      end
   end
end

# see https://discourse.julialang.org/t/printf-with-variable-format-string/3805/3
print_formatted(io, fmt, args...) = @eval @printf($io, $fmt, $(args...))


text:
const ManyMatrices = AbstractVecOrTuple{AbstractMatrix}

@inline inner(A,B) = dot(A,B)

@inline inner(A::NTuple{2,Operator}) = @inbounds inner(A[1],A[2])

@inline inner(A::Tuple{Vararg{Operator}},B::Operator) = inner_cs(data.(A),data(B))

function inner_cs(A,B::AbstractMatrix,s::BitSet=BitSet())
    # calculate |trace(A'*B)|cs (coherent subspaces)
    # this inner product cares only about relative phase between the subspaces contained in s
    x,y = _inner_cs_1(A,B,s)
    return _inner_cs_2(x,y)
end

function inner_cs_grad(Pj,JXj::AbstractMatrix,x,y,s::BitSet=BitSet())
    # calculate the partial derivative of |trace(A'*B)|cs (coherent subspaces)
    w,z = _inner_cs_1(Pj,JXj,s)
    res = real(w*x)
    @inbounds for i = 1:length(z)
        res += real(z[i]*y[i])
        res += real(w*normalize(x))*abs(y[i]) + real(z[i]*normalize(y[i]))*abs(x)
        for j = 1:i-1
            res += real(z[i]*normalize(y[i]))*abs(y[j]) + real(z[j]*normalize(y[j]))*abs(y[i])
        end
    end
    return res/_inner_cs_2(x,y)
end

# sub-routines

function _inner_cs_1(A::AbstractMatrix{T},B::AbstractMatrix{S},s::BitSet) where {T,S}
    m, n = size(A)
    size(B) == (m,n) || throw(DimensionMismatch("matrices must have the same dimensions"))
    x = zero(promote_type(T,S))
    y = Vector{typeof(x)}(undef,n-length(s))
    t = 1
    @inbounds for j = 1:n
        a = zero(x)
        for i = 1:m
            a += A[i,j]' * B[i,j]
        end
        if j in s
            x += a
        else
            y[t] = a
            t += 1
        end
    end
    return x, y
end

function _inner_cs_1(A::ManyMatrices,B::AbstractMatrix,::BitSet)
    # calculate |trace(A'*B)|cs (coherent subspaces)
    # this inner product cares only about relative phase between the subspaces contained in s
    x = inner.(A,(B,))
    return first(x), tail(x)
end

function _inner_cs_1(A::AbstractMatrix,B::ManyMatrices,::BitSet)
    # calculate |trace(A'*B)|cs (coherent subspaces)
    # this inner product cares only about relative phase between the subspaces contained in s
    x = inner.((A,),B)
    return first(x), tail(x)
end

function _inner_cs_2(x::Number,y::Union{Vector{<:Number},NTuple{M,<:Number} where M})
    res = abs2(x)
    @inbounds for i = 1:length(y)
        res += abs2(y[i])
        res += 2*abs(y[i])*abs(x)
        for j = 1:i-1
            res += 2*abs(y[i])*abs(y[j])
        end
    end
    return sqrt(res)
end

f_cs_2(x::Number,y::Union{Vector{<:Number},NTuple{M,<:Number} where M}) = sum(abs,y)+abs(x)


text:
module ImageBase

export
    # two-fold downsampling
    # originally from ImageTransformations.jl
    restrict,

    # basic image statistics, from Images.jl
    minimum_finite,
    maximum_finite,
    meanfinite,
    varfinite,
    sumfinite

# Introduced in ColorVectorSpace v0.9.3
# https://github.com/JuliaGraphics/ColorVectorSpace.jl/pull/172
using ImageCore.ColorVectorSpace.Future: abs2

using Reexport

using Base.Cartesian: @nloops
@reexport using ImageCore
using ImageCore.OffsetArrays

include("diff.jl")
include("restrict.jl")
include("utils.jl")
include("statistics.jl")
include("deprecated.jl")

if VERSION >= v"1.4.2" # work around https://github.com/JuliaLang/julia/issues/34121
    include("precompile.jl")
    _precompile_()
end

end


text:
function findfirst(f::File,r::Regex;from::Int=1,to::Int=length(f.data))
    for l in from:to
        if occursin(r,f.data[l])
            return l
        end
    end
    return nothing
end

findnext(f::File,r::Regex,ln::Int) = findfirst(f,r,from=ln)

function finduntil(f::File,r::Regex;from::Int=1,to::Int=length(f.data))
    for l in from:to
        if !occursin(r,f.data[l])
            return l-1
        end
    end
    return to
end

function findall(f::File,r::Regex;from::Int=1,to::Int=length(f.data))
    res = Vector{Int}()
    for l in from:to
        if occursin(r,f.data[l])
            push!(res,l)
        end
    end
    if length(res) == 0
        return nothing
    else
        return res
    end
end

parseSpecifier(f::File,s::Line,::Val{:after})= s.ln+1
parseSpecifier(f::File,s::Int,::Val{:after}) = s+1
parseSpecifier(f::File,s::Line,::Val{:before})= s.ln-1
parseSpecifier(f::File,s::Int,::Val{:before}) = s-1
parseSpecifier(f::File,s::Line,::Val{T}) where T= s.ln
parseSpecifier(f::File,s::Int,::Val{T}) where T= s
parseSpecifier(f::File,s::Nothing,::Val{:after}) = 1
parseSpecifier(f::File,s::Nothing,::Val{:before}) = length(f.data)

function Line(f::File,r::Regex;
              after::A=nothing,
              before::B=nothing
              ) where {A<:Union{Regex,Line,Int,Nothing}, B<:Union{Regex,Line,Int,Nothing}}
    if A == Regex
        _after = findfirst(f,after)
    else
        _after = parseSpecifier(f,after,Val{:after}())
    end

    _after == nothing && error("After ($after) was not found in the file")
    _after > length(f.data) && error("After is larger then the files length.")

    if B == Regex
        _before = findfirst(f,before,from=_after+1)
    else
        _before = parseSpecifier(f,before,Val{:before}())
    end

    _before == nothing && error("Before ($before) was not found in the file")
    _before > length(f.data) && error("Before is larger then the files length.")
    _after > _before && error("After must be smaller than before!")
    ind = findfirst(f,r,from=_after,to=_before)
    ind == nothing && error("$r was not found in the file")
    return Line(f,ind)
end

function Range(f::File;
               from::F=nothing,
               to::T=nothing,
               until::U=nothing,
               after::A=nothing,
               before::B=nothing) where {
                   F<:Union{Regex,Line,Int,Nothing},
                   T<:Union{Regex,Line,Int,Nothing},
                   U<:Union{Regex,Nothing},
                   A<:Union{Regex,Line,Int,Nothing},
                   B<:Union{Regex,Line,Int,Nothing}
               }
    (T != Nothing && U != Nothing) && error("A range cannot be specified by to AND until.")
    (T == Nothing && U == Nothing) && error("A range must be definied with to or until.")

    # after --> from ----> to ----> before
    # after --> from ---> until ---> before

    if A == Regex
        _after = findfirst(f,after)
    else
        _after = parseSpecifier(f,after,Val{:after}())
    end

    _after > length(f.data) && error("After is larger then the files length.")
    _after == nothing && error("After ($after) was not found in the file")

    if B == Regex
        _before = findfirst(f,before,from=_after+1)
    else
        _before = parseSpecifier(f,before,Val{:before}())
    end

    _before == nothing && error("Before ($before) was not found in the file")
    _before > length(f.data) && error("Before is larger then the files length.")
    _after > _before && error("After must be smaller than before!")

    if F == Regex
        _from = findfirst(f,from,from=_after,to=_before)
        _from == nothing && error("From ($from) was not found in the file")
    elseif from != nothing
        _from = parseSpecifier(f,from,Val{:from}())
    else
        _from = 1
    end

    _before < _from && error("From must be smaller than before!")

    if U == Regex
        _to = finduntil(f,until,from=max(_from,_after)+1,to=_before)
        _to == nothing && error("Until ($until) was not found in the file")
    elseif T == Regex
        _to = findfirst(f,to,from=max(_from,_after)+1,to=_before)
        _to == nothing && error("To ($to) was not found in the file")
    elseif isa(to,Int)
        _to = _from + to
    elseif isa(to,Line)
        _to = to.ln
    elseif to != nothing
        _to = parseSpecifier(f,from,Val{:to}())
    else
        _to = _before
    end

    _after > _to && error("To must be larger than after!")
    _before < _to && error("To must be smaller than before!")
    _from > _to && error("To must be larger than from!")

    Range(f,_from,_to)
end

function Lines(f::File,r::Regex;
              after::A=nothing,
              before::B=nothing
              ) where {A<:Union{Regex,Line,Int,Nothing}, B<:Union{Regex,Line,Int,Nothing}}
    if A == Regex
        _after = findfirst(f,after)
    else
        _after = parseSpecifier(f,after,Val{:after}())
    end

    _after == nothing && error("After ($after) was not found in the file")
    _after > length(f.data) && error("After is larger then the files length.")

    if B == Regex
        _before = findfirst(f,before,from=_after+1)
    else
        _before = parseSpecifier(f,before,Val{:before}())
    end

    _before == nothing && error("Before ($before) was not found in the file")
    _before > length(f.data) && error("Before is larger then the files length.")
    _after > _before && error("After must be smaller than before!")
    ind = findall(f,r,from=_after,to=_before)
    ind == nothing && error("$r was not found in the file")
    return Line.(Ref{File}(f),ind) |> Collection
end

function Ranges(f::File,r::Regex;
              after::A=nothing,
              before::B=nothing
              ) where {A<:Union{Regex,Line,Int,Nothing}, B<:Union{Regex,Line,Int,Nothing}}
    if A == Regex
        _after = findfirst(f,after)
    else
        _after = parseSpecifier(f,after,Val{:after}())
    end

    _after == nothing && error("After ($after) was not found in the file")
    _after > length(f.data) && error("After is larger then the files length.")

    if B == Regex
        _before = findfirst(f,before,from=_after+1)
    else
        _before = parseSpecifier(f,before,Val{:before}())
    end

    _before == nothing && error("Before ($before) was not found in the file")
    _before > length(f.data) && error("Before is larger then the files length.")
    _after > _before && error("After must be smaller than before!")
    ind = findall(f,r,from=_after,to=_before)
    ind == nothing && error("$r was not found in the file")

    ret = Collection(Reference[])

    start = ind[1]
    last = ind[1]
    for i in ind[2:end]
        if i != (last+1) || i == ind[end]
            if start < last
                push!(ret.refs,Range(f,start,i == ind[end] ? i : last))
            else
                push!(ret.refs,Line(f,start))
                if i == ind[end]
                    push!(ret.refs,Line(f,i))
                end
            end
            start = i
            last = i
        else
            last = i
        end
    end
    return ret
end


text:
# Reads in the string that describes the hierarchy,
# and sets up A matrix and masses for each Keplerian.
# There are N bodies, and N-1 Keplerians.

function get_indices(string)
# Takes a Kepler hierarchy vector, and returns indices of all planets
ind1 = readdlm(IOBuffer(hierarchy[1:i]))
# Converts strings to integers
nm = size(string)
for j=1:nm
  replace(string[j],"(","")
  replace(string[j],")","")
end
indices = zeros(Int64,nm)
for j=1:nm
  indice[j]=parse(string[j])
end
return indices
end

function split_kepler(hierarchy,mass)
# Splits the hierarchy into components
ndelim = count(c -> c == ',' , hierarchy)
if ndelim > 1
  # Find first complete Keplerian:
  i = 1
  nleft = 0
  nright = 0
  c = collect(hierarchy)
  nchar = endof(hierarchy)
  while (nleft == 0) || (nleft != nright) || i == nchar
    if c[i] == '('
      nleft += 1
    end
    if c[i] == ')'
      nright += 1
    end
    if c[i] == ','
      ndelim += 1
    end
  end
  @assert(nleft == nright)
  @assert(nleft == ndelim)
  # Compute sum of masses of both Kepler components:
  # First find the indices of planets:
  ind1 = get_indices(hierarchy[1:i])
  ind2 = get_indices(hierarchy[i+1:nchar])
  m1 = sum(mss[ind1])
  m2 = sum(mss[ind2])
  return m1,m2,hierarchy[1:i],hierarchy[i+1:nchar]
elseif ndelim == 1
  # Found a Keplerian of two bodies
  ind = get_indices(hierarchy)
  return mass[ind[1]],mass[ind[2]],split(hierarchy,",")
else
  # Down to a single body
  ind = get_indices(hierarchy)
  return mass[ind[1]],0.0,hierarchy,""
end
end

function strip_parentheses(hierarchy)
# Strips the leading and trailing parentheses (if they exist):
nchar = endof(hierarchy)
if hierarchy[1] == '(' && hierarchy[nchar] == ')'
  return hierarchy[2:nchar-1]
else
  return hierarchy
end
end

nleft  = count(c -> c == '(',hierarchy)
nright = count(c -> c == ')',hierarchy)
@assert(nleft == nright)
return
end

function setup_hierarchy(hierarchy,nbody,mass)
# Creates a routine which computes the matrix
# and total masses in each Keplerian.
nleft  = count(c -> c == '(',hierarchy)
nright = count(c -> c == ')',hierarchy)
ndelim = count(c -> c == ',',hierarchy)
@assert(nleft == nright)
@assert(nleft == ndelim)

return A,mtot
end


text:
# load requried packages
using IncrementalInference
using Test

## parameters

lm_prior_noise = 0.1
meas_noise = 0.15
odom_noise = 0.1
n_samples = 100

# initialize mean landmark locations
l0 = 0.0
l1 = 10.0
l2 = 40.0

# "Ground-truth" robot poses
x0 = 0.0
x1 = 10.0
x2 = 20.0
x3 = 40.0

## Initialize empty factor graph
fg = initfg()

# Place strong prior on locations of three "doors"
addVariable!(fg, :l0, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l0], Prior(Normal(l0, lm_prior_noise)))

addVariable!(fg, :l1, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l1], Prior(Normal(l1, lm_prior_noise)))

addVariable!(fg, :l2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l2], Prior(Normal(l2, lm_prior_noise)))

# Add first pose
addVariable!(fg, :x0, ContinuousScalar, N=n_samples)

# Make first "door" measurement
addFactor!(fg, [:x0; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0])

# Add second pose
addVariable!(fg, :x1, ContinuousScalar, N=n_samples)

# Gaussian transition model
addFactor!(fg, [:x0; :x1], LinearRelative(Normal(x1-x0, odom_noise)))

# Make second "door" measurement
addFactor!(fg, [:x1; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0])

## Add one more pose/odometry to invoke issue #236

# Add third pose
addVariable!(fg, :x2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(x2-x1, odom_noise)))

# Add fourth pose
# addVariable!(fg, :x3, ContinuousScalar, N=n_samples)

# Add odometry transition and new landmark sighting
# addFactor!(fg, [:x2, :x3], LinearRelative(Normal(2, odom_noise)))
# addFactor!(fg, [:x3; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0])

## Do some debugging

ensureAllInitialized!(fg)

drawGraph(fg, show=true)

tree = buildTreeReset!(fg, drawpdf=true, show=true)

## Solve graph
tree = solveTree!(fg)

## Plotting functions below

using RoMEPlotting

plotKDE(fg, [:x0])

plotKDE(fg, [:x0;:x1])

plotKDE(fg, [:x0;:x1;:x2])

plotKDE(fg, [:l0;:l1])

plotKDE(fg, [:l2])

spyCliqMat(tree, :l0)

spyCliqMat(tree, :x2)

## swap iteration order

#TODO
# getCliqueData(tree.cliques[2]).itervarIDs = [9;7;1;3;5]
getCliqueData(tree.cliques[1]).itervarIDs = [:x1, :x0, :l2, :l1, :l0]

solveTree!(fg, tree)

## Reconstruct individual steps for broader clique factor selection

# Cliq 2:
# L0,X0,L1,X1
# x ,  ,  ,
# x ,x ,x ,
# x ,  ,x ,x
#   ,x ,  ,x
#   ,  ,x ,

# choose iteration order (priors last): :x1, :x0, :l2, :l1, :l0

# get factors for :x1 in clique2:
# :x1l0l1l2f1, :x0x1f1
plotLocalProduct(fg, :x1, sidelength=20cm)
ptsX1, = predictbelief(fg, :x1, [:x1l0l1l2f1, :x0x1f1])
X1 = manikde!(ptsX1, (:Euclid,))
plotKDE([X1; getKDE(fg, :x1)], c=["red";"green"])
setValKDE!(fg, :x1, X1)

# get factors for :x0 in clique2:
# :x0l0l1l2f1, :x0x1f1
plotLocalProduct(fg, :x0, sidelength=20cm)
ptsX0, = predictbelief(fg, :x0, [:x0l0l1l2f1; :x0x1f1])
X0 = manikde!(ptsX0, (:Euclid,))
plotKDE([X0; getKDE(fg, :x0)], c=["red";"green"])
setValKDE!(fg, :x0, X0)

# get factors for :l2
# :x0l0l1l2f1, :x1l0l1l2f1, :l2f1
plotLocalProduct(fg, :l2, sidelength=20cm)
ptsL2, = predictbelief(fg, :l2, [:x0l0l1l2f1, :x1l0l1l2f1, :l2f1])
L2 = manikde!(ptsL2, (:Euclid,))
plotKDE([L2; getKDE(fg, :l2)], c=["red";"green"])
setValKDE!(fg, :l2, L2)

# get factors for :l1
# :x0l0l1l2f1, :x1l0l1l2f1, :l1f1
plotLocalProduct(fg, :l1, sidelength=20cm)
ptsL1, = predictbelief(fg, :l1, [:x0l0l1l2f1, :x1l0l1l2f1, :l1f1])
L1 = manikde!(ptsL1, (:Euclid,))
plotKDE([L1; getKDE(fg, :l1)], c=["red";"green"])
setValKDE!(fg, :l1, L1)

# get factors for :l0
# :x0l0l1l2f1, :x1l0l1l2f1, :l0f1
plotLocalProduct(fg, :l0, sidelength=20cm)
ptsL0, = predictbelief(fg, :l0, [:x0l0l1l2f1, :x1l0l1l2f1, :l0f1])
L0 = manikde!(ptsL0, (:Euclid,))
plotKDE([L0; getKDE(fg, :l0)], c=["red";"green"])
setValKDE!(fg, :l0, L0)

## Now do root clique manually too

ptsX2, = predictbelief(fg, :x2, [:x1x2f1;])
X2 = manikde!(ptsX2, (:Euclid,))
plotKDE([X2; getKDE(fg, :x2)], c=["red";"green"])
setValKDE!(fg, :x2, X2)

upmsgX1 = deepcopy(getKDE(fg, :x1))
ptsX1 = approxConv(fg, :x1x2f1, :x1)
pX1 = manikde!(ptsX1, (:Euclid,))
X1 = manifoldProduct([upmsgX1; pX1], (:Euclid,))
plotKDE([X1; getKDE(fg, :x1)], c=["red";"green"])
setValKDE!(fg, :x1, X1)

##  Complete downward pass

# :x0l0l1l2f1, :x1l0l1l2f1, :l1f1
ptsL1, = predictbelief(fg, :l1, [:x0l0l1l2f1, :x1l0l1l2f1, :l1f1])
L1 = manikde!(ptsL1, (:Euclid,))
plotKDE([L1; getKDE(fg, :l1)], c=["red";"green"])
setValKDE!(fg, :l1, L1)

## quick debug tests

##

pts = approxConv(fg, :x1l0l1l2f1, :l1)
plotKDE(manikde!(pts, (:Euclid,)))

##

plotKDE(fg, :x0)
pts = approxConv(fg, :x0l0l1l2f1, :l1)
plotKDE(manikde!(pts, (:Euclid,)))

## Debug draft-tree based autoinit:
# x0, l2, l1, l0, with one multihypo factor and 3 priors
# xxxx
#  x
#   x
#    x

fg = initfg()
addVariable!(fg, :l0, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l0], Prior(Normal(l0, lm_prior_noise)), graphinit=false)
addVariable!(fg, :l1, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l1], Prior(Normal(l1, lm_prior_noise)), graphinit=false)
addVariable!(fg, :l2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l2], Prior(Normal(l2, lm_prior_noise)), graphinit=false)
addVariable!(fg, :x0, ContinuousScalar, N=n_samples)
addFactor!(fg, [:x0; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0], graphinit=false)
addVariable!(fg, :x1, ContinuousScalar, N=n_samples)
addFactor!(fg, [:x0; :x1], LinearRelative(Normal(x1-x0, odom_noise)),graphinit=false)
addFactor!(fg, [:x1; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0],graphinit=false)
addVariable!(fg, :x2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(x2-x1, odom_noise)),graphinit=false)

drawGraph(fg, show=true)

tree = buildTreeReset!(fg, drawpdf=true, show=true)

spyCliqMat(tree, :l0)
spyCliqMat(tree, :x2)

cliq = tree.cliques[2]

## quick debug

##

initManual!(fg, :l2, [:l2f1])
initManual!(fg, :l1, [:l1f1])
initManual!(fg, :l0, [:l0f1])

# regular procedure
ptsX1, = predictbelief(fg, :x1, [:x1l0l1l2f1])
X1 = manikde!(ptsX1, (:Euclid,))
plotKDE(X1)
# just because the init flag has not been set yet
setValKDE!(fg, :x1, X1)

# regular procedure
ptsX0, = predictbelief(fg, :x0, [:x0l0l1l2f1])
X0 = manikde!(ptsX0, (:Euclid,))
plotKDE(X0)
# just because the init flag has not been set yet
setValKDE!(fg, :x0, X0)

# regular procedure
ptsX2, = predictbelief(fg, :x2, [:x1x2f1])
X2 = manikde!(ptsX2, (:Euclid,))
plotKDE(X2)
# just because the init flag has not been set yet
setValKDE!(fg, :x2, X2)

## other debugging

plotLocalProduct(fg, :x0)
plotLocalProduct(fg, :x1)
plotLocalProduct(fg, :l1, sidelength=15cm)

##

stuff = treeProductUp(fg, tree, :x2, :x2)

##

ensureAllInitialized!(fg)

tree = buildTreeReset!(fg, drawpdf=true, show=true)
cliqorder = getEliminationOrder(tree)

#FIXME spyCliqMat(cliqorder[end])

#


text:
"""
    blockmodel(inblock::ImageTensor{N}, outblock::OneHotTensor{0}, backbone)
    blockmodel(inblock::ImageTensor{N}, outblock::OneHotTensorMulti{0}, backbone)

Construct a model for N-dimensional image classification. `backbone` should
be a convolutional feature extractor taking in batches of image tensors with
`inblock.nch` color channels.
"""
function blockmodel(
        inblock::ImageTensor{N},
        outblock::Union{OneHotTensor{0}, OneHotTensorMulti{0}},
        backbone) where N
    outsz = Flux.outputsize(backbone, (ntuple(_ -> 256, N)..., inblock.nchannels, 1))
    outch = outsz[end-1]
    head = Models.visionhead(outch, length(outblock.classes), p = 0.)
    return Flux.Chain(backbone, head)
end

"""
    blockmodel(inblock::ImageTensor{N}, outblock::OneHotTensor{N}, backbone; kwargs...)

Construct a model for N-dimensional image segmentation. `backbone` should
be a convolutional feature extractor taking in batches of image tensors with
`inblock.nch` color channels. Keyword arguments are passed to [`UNetDynamic`](#).
"""
function blockmodel(inblock::ImageTensor{N}, outblock::OneHotTensor{N}, backbone; kwargs...) where N
    return Models.UNetDynamic(
        backbone,
        (ntuple(_ -> 256, N)..., inblock.nchannels, 1),
        length(outblock.classes);
        kwargs...)
end

"""
    blockmodel(inblock::ImageTensor{N}, outblock::Keypoints{N}, backbone)

Construct a model for image to keypoint regression. `backbone` should
be a convolutional feature extractor taking in batches of image tensors with
`inblock.nch` color channels.
"""
function blockmodel(inblock::ImageTensor{N}, outblock::KeypointTensor{N}, backbone) where N
    outsz = Flux.outputsize(backbone, (ntuple(_ -> 256, N)..., inblock.nchannels, 1))
    outch = outsz[end-1]
    head = Models.visionhead(outch, prod(outblock.sz)*N, p = 0.)
    return Flux.Chain(backbone, head)
end

"""
    blockbackbone(ImageTensor{2}(ch))

Construct a XResNet18 model that takes in an encoded image with `ch`
color channels.
"""
blockbackbone(inblock::ImageTensor{2}) = Models.xresnet18(c_in = inblock.nchannels)

# ## Tests

@testset "blockbackbone" begin
    @test_nowarn FastAI.blockbackbone(Vision.ImageTensor{2}(3))
end


text:
# functions that return information about h2o

function frames()
    ## Retrieve all frames in the h2o cluster
    results = h2o_get("Frames")
	result_json = JSON.parse(readall(results))

    @printf "%-34s" "frame_id"
    @printf "%15s" "rows"
    @printf "%15s" "columns"
    @printf "%16s" "size (MB)\n"

    for i in 1:length(result_json["frames"])
        @printf "%-34s" result_json["frames"][i]["frame_id"]["name"]
        @printf "%15s" result_json["frames"][i]["rows"]
        @printf "%15s" result_json["frames"][i]["columns"]
        @printf "%15.2f" result_json["frames"][i]["byte_size"]/1000000
        @printf "MB\n"
    end

# need to decide later, as a whole for H2O.jl when we want to return
# results vs. just print
# perhaps bifurcate into two separate types of functions:

# list_frames(), which would simply print
# get_frames(), which would return a JSON object

#    return result_json

end


text:
###############################################################################
#
# Flow mode, steady state or non-steady state
#
###############################################################################
"""
    abstract type AbstractFlowMode

Hierarchy of AbstractFlowMode
- [`SteadyStateMode`](@ref)
- [`NonSteadyStateMode`](@ref)
"""
abstract type AbstractFlowMode end

"""
    struct SteadyStateMode <: AbstractFlowMode
"""
struct SteadyStateMode <: AbstractFlowMode end

"""
    struct NonSteadyStateMode <: AbstractFlowMode
"""
struct NonSteadyStateMode <: AbstractFlowMode end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, InfiniteArrays, Test, Base64, RecipesBase
import ContinuumArrays: ℵ₁, materialize, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, ℵ₁,
                        MappedBasisLayout, AdjointMappedBasisLayout, MappedWeightedBasisLayout, TransformFactorization, Weight, WeightedBasisLayout, SubWeightedBasisLayout, WeightLayout,
                        basis, invmap, Map, checkpoints, _plotgrid, mul
import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout, LdivStyle, MulStyle

@testset "Inclusion" begin
    @testset "basics" begin
        x = Inclusion(-1..1)
        @test eltype(x) == Float64
        @test x[0.1] ≡ 0.1
        @test x[0] ≡ x[0.0] ≡ 0.0
    end

    @testset "broadcast" begin
        x = Inclusion(-1.0..1)
        X = QuasiDiagonal(x)
        @test X[-1:0.1:1,-1:0.1:1] == Diagonal(-1:0.1:1)
        @test Base.BroadcastStyle(typeof(x)) == LazyQuasiArrayStyle{1}()
        @test x .* x isa BroadcastQuasiArray
        @test (x.*x)[0.1] == 0.1^2

        @test exp.(x)[1.0] == exp(1.0)
        @test exp.(-x)[1.0] == exp(-1.0)
        @test exp.(x.^2)[1.0] == exp(1.0)
        @test exp.(-x.^2/(20/2))[1.0] == exp(-1.0^2/(20/2))

        @test dot(x,x) ≈ 2/3
        @test norm(x) ≈ sqrt(2/3)
    end

    @testset "Derivative" begin
        x = Inclusion(-1..1)
        D = Derivative(x)
        @test D == Derivative{Float64}(x) == Derivative{Float64}(-1..1)
        @test D*x ≡ QuasiOnes(x)
        @test D^2 * x ≡ QuasiZeros(x)
        @test D*[x D*x] == [D*x D^2*x]
        @test stringmime("text/plain", D) == "Derivative(Inclusion(-1..1))"
    end
end

include("test_maps.jl")

@testset "DiracDelta" begin
    δ = DiracDelta()
    @test δ == DiracDelta(0, Inclusion(0))
    @test axes(δ) ≡ (axes(δ,1),) ≡ (Inclusion(0.0),)
    @test size(δ) ≡ (length(δ),) ≡ (1,)
    @test_throws BoundsError δ[1.1]
    @test δ[0.0] ≡ Inf
    @test Base.IndexStyle(δ) ≡ Base.IndexLinear()

    δ = DiracDelta(0, -1..3)
    @test δ == DiracDelta{Float64}(0, -1..3)
    @test axes(δ) ≡ (axes(δ,1),) ≡ (Inclusion(-1..3),)
    @test size(δ) ≡ (length(δ),) ≡ (ℵ₁,)
    @test δ[1.1] ≡ 0.0
    @test δ[0.0] ≡ Inf
    @test Base.IndexStyle(δ) ≡ Base.IndexLinear()

    @test stringmime("text/plain", δ) == "δ at 0.0 over Inclusion(-1..3)"
    x = Inclusion(-1..1)
    @test stringmime("text/plain", δ[2x .+ 1]) == "δ at 0.0 over Inclusion(-1..3) affine mapped to -1..1"
end

@testset "Kernels" begin
    x = Inclusion(0..1)
    K = x .- x'
    @test K[0.1,0.2] == K[Inclusion(0..0.5), Inclusion(0..0.5)][0.1,0.2] == 0.1 - 0.2
    @test_throws BoundsError K[Inclusion(0..0.5), Inclusion(0..0.5)][1,1]

    L = LinearSpline(0:0.1:1)
    @test_broken L\K isa ApplyQuasiArray # broken since it assumes non-broadcasting
    @test L \ exp.(K) isa ApplyQuasiArray
end

include("test_splines.jl")
include("test_chebyshev.jl")
include("test_basisconcat.jl")

@testset "Plotting" begin
    L = LinearSpline(0:5)
    rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), L)
    @test rep[1].args == (L.points,L[L.points,:])

    u = L*randn(6)
    rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
    @test rep[1].args == (L.points,u[L.points])

    @testset "padded" begin
        u = L * Vcat(rand(3), Zeros(3))
        rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
        @test rep[1].args == (L.points,u[L.points])
    end

    @testset "Chebyshev and weighted Chebyshev" begin
        T =  Chebyshev(10)
        w =  ChebyshevWeight()
        wT = w .* T
        x =  axes(T, 1)

        u = T * Vcat(rand(3), Zeros(7))
        v = wT * Vcat(rand(3), Zeros(7))

        rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
        @test rep[1].args == (grid(T), u[grid(T)])
        wrep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), v)
        @test wrep[1].args == (grid(wT), v[grid(wT)])

        @test plotgrid(v) == plotgrid(u) == grid(T) == grid(wT) == _plotgrid(MemoryLayout(v), v) == _plotgrid(MemoryLayout(u), u)
        y = affine(0..1, x)
        @test plotgrid(T[y,:]) == (plotgrid(T) .+ 1)/2
    end
end


text:
using Test
using Provenance