text:
using MacroTools: isdef

@testset "utils" begin
    ex1 = :(function foo(a) return a; end)
    @test isdef(ex1)
    ex2 = :(function bar(a)::Int return 1; end)
    @test isdef(ex2)
    ex3 = :(function foo(a::T) where T return a; end)
    @test isdef(ex3)
    ex4 = :(function bar(a::T)::Int where T return 1; end)
    @test isdef(ex4)
    ex5 = :(function bar(a::S, b::T)::Union{S,T} where {S,T} if rand() < 0.5 return a; end; return b; end)
    @test isdef(ex5)

    ex6 = :(f(a) = a)
    @test isdef(ex6)
    ex7 = :(f(a)::Int == 1)
    @test isdef(ex7)
    ex8 = :(f(a::T) where T = a)
    @test isdef(ex8)
    ex9 = :(f(a::T)::Int where T = 1)
    @test isdef(ex9)
    ex10 = :(f(a::S, b::T)::Union{S,T} where {S,T} = rand() < 0.5 ? a : b)
    @test isdef(ex10)
end


text:
"""
    Scan(dataSource::DataSource, intervals::Vector{<:Interval}; <keyword arguments>)

A method with all arguments as keyword arguments is also provided.

# Arguments
- columns::Vector{String} = nothing
- filter::Filter = nothing
- virtualColumns::Vector{<:VirtualColumn} = nothing
- order::String = nothing
- limit::Integer = nothing
- offset::Integer = nothing
- resultFormat::String = nothing
- batchSize::Integer = nothing
- context::Dict = nothing
- legacy::Bool = nothing
"""
mutable struct Scan <: Query
    queryType::String
    dataSource::DataSource
    intervals::Vector{<:Interval}
    columns
    filter
    virtualColumns
    order
    limit
    offset
    resultFormat
    batchSize
    context
    legacy
    function Scan(
        dataSource, intervals;
        columns=nothing, filter=nothing, virtualColumns=nothing, order=nothing,
        limit=nothing, offset=nothing, resultFormat=nothing,
        batchSize=nothing, context=nothing, legacy=nothing
    )
        nothing_or_type(columns, Vector{String})
        nothing_or_type(filter, Filter)
        nothing_or_type(virtualColumns, Vector{<:VirtualColumn})
        order === nothing || (order = lowercase(order)) ∈ ["ascending", "descending", "none"] || error("Invalid order")
        limit === nothing || (isa(limit, Integer) && limit >= 0) || error("limit must be a non-negative integer")
        offset === nothing || (isa(offset, Integer) && offset >= 0) || error("offset must be a non-negative integer")
        resultFormat === nothing || resultFormat ∈ ["list", "compactedList"] || error("Invalid resultFormat")
        batchSize === nothing || (isa(batchSize, Integer) && batchSize >= 0) || error("batchSize must be a non-negative integer")
        nothing_or_type(context, Dict)
        nothing_or_type(legacy, Bool)
        new("scan", dataSource, intervals, columns, filter, virtualColumns, order, limit, offset, resultFormat, batchSize, context, legacy)
    end
end
Scan(
    ; dataSource, intervals,
    columns=nothing, filter=nothing, virtualColumns=nothing, order=nothing,
    limit=nothing, offset=nothing, resultFormat=nothing,
    batchSize=nothing, context=nothing, legacy=nothing
) = Scan(
    dataSource, intervals;
    columns, filter, virtualColumns, order, limit, offset, resultFormat,
    batchSize, context, legacy
)

function execute(client::Client, query::Scan; pretty=false)
    res = convert(Vector{Dict}, JSON.parse(execute_native_query(client, query; pretty)))
    if length(res) != 0
        names = Symbol.(res[1]["columns"])
        num_rows = sum(length(subres["events"]) for subres ∈ res)
    else
        names = Symbol[]
        num_rows = 0
    end
    ScanResult(names, num_rows, res)
end

struct ScanResult <: QueryResult
    names::Vector{Symbol}
    num_rows::Int
    inner_array::Vector{Dict}
end

names(sr::ScanResult) = getfield(sr, :names)

Base.summary(io::IO, sr::ScanResult) = print(io, string(length(sr)) * "-element " * string(typeof(sr)))
Base.getindex(sr::ScanResult, i::Int) = (i <= length(sr) || throw(BoundsError(sr, i))) && ScanRow(i, sr)

Tables.rowaccess(::ScanResult) = true
Tables.rows(sr::ScanResult) = sr

Base.eltype(::ScanResult) = ScanRow
Base.length(sr::ScanResult) = getfield(sr, :num_rows)
Base.iterate(sr::ScanResult, state=1) = state > length(sr) ? nothing : (ScanRow(state, sr), state + 1)

function Tables.schema(sr::ScanResult)
    if length(sr) != 0
        types = [typeof(cl) for cl in sr[1]]
    else
        types = DataType[]
    end
    Tables.Schema(names(sr), types)
end

struct ScanRow <: Tables.AbstractRow
    row::Int
    source::ScanResult
end

function Tables.getcolumn(sr::ScanRow, name::Symbol)
    elems = 0
    i = getfield(sr, :row)
    for subres ∈ getfield(getfield(sr, :source), :inner_array)
        subelems = length(subres["events"])
        if i <= elems + subelems
            return subres["events"][i - elems]
        end
        elems += subelems
    end
end

Tables.getcolumn(sr::ScanRow, i::Int) = Tables.getcolumn(sr, names(getfield(sr, :source))[i])
Tables.columnnames(sr::ScanRow) = names(getfield(sr, :source))

ext: jl


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

# given a route, a priority for travel, a start time, and an end location,
# (and an end node which should be the node nearest to end location),
# modify current route
function changeRoute!(sim::Simulation, route::Route, priority::Priority, startTime::Float, endLoc::Location, endFNode::Int)

	# shorthand:
	@unpack map, net, travel = sim

	# get data on current route before changing
	startLoc = getRouteCurrentLocation!(net, route, startTime)

	travelMode = getTravelMode!(travel, priority, startTime)

	(startFNode, startFNodeTravelTime) = getRouteNextNode!(sim, route, travelMode.index, startTime)
	startFNodeTime = startTime + startFNodeTravelTime
	startFNodeDist = getRouteNextNodeDist!(sim, route, startTime)

	(pathTravelTime, rNodes) = shortestPathData(net, travelMode.index, startFNode, endFNode)

	# shorthand:
	fNetTravel = travelMode.fNetTravel
	@unpack fNodeFromRNodeTime, fNodeToRNodeTime = fNetTravel

	## change route

	route.priority = priority
	route.travelModeIndex = travelMode.index

	# start and end fNodes, times, and distances
	route.startFNode = startFNode
	route.startFNodeTime = startFNodeTime
	route.startFNodeDist = startFNodeDist
	route.endFNode = endFNode
	route.endFNodeTime = startFNodeTime + pathTravelTime

	# start and end rNodes and times
	route.startRNode = rNodes[1]
	route.endRNode = rNodes[2]
	if route.startRNode != nullIndex
		@assert(route.endRNode != nullIndex)
		route.startRNodeTime = startFNodeTime + fNodeToRNodeTime[startFNode][route.startRNode]
		route.endRNodeTime = route.startRNode == route.endRNode ? route.startRNodeTime : route.endFNodeTime - fNodeFromRNodeTime[endFNode][route.endRNode]
	else
		route.startRNodeTime = nullTime
		route.endRNodeTime = nullTime
	end

	# start and end locations and times
	copy!(route.startLoc, startLoc)
	route.startTime = startTime
	copy!(route.endLoc, endLoc)
	route.endTime = route.endFNodeTime + offRoadTravelTime(travelMode, map, net.fGraph.nodes[endFNode].location, endLoc)

	# recent rArc, recent fNode, next fNode, status
	setRouteStateBeforeStartFNode!(route, startTime)

	# first rArc
	setRouteFirstRArc!(net, route)
end

# Initialise an empty route to be at a given location,
# along with the start node to travel to and distance to that node.
function initRoute!(sim::Simulation, route::Route;
	startLoc::Location = Location(), startFNode::Int = nullIndex, startFNodeDist::Float = nullDist)

	@assert(route.status == routeNullStatus)
	@assert(startLoc != Location())
	@assert(startFNode != nullIndex)
	@assert(startFNodeDist >= 0.0)

	# make route that starts at time = Inf
	route.startTime = Inf
	copy!(route.startLoc, startLoc)
	route.startFNode = startFNode
	route.startFNodeDist = startFNodeDist
	route.startFNodeTime = Inf
	# route.endTime = Inf # leave as nullTime, for getRouteNextNode!
	copy!(route.endLoc, startLoc) # needed for animation
	route.endFNode = startFNode
	route.endFNodeTime = Inf
	route.nextFNode = startFNode
	setRouteStateBeforeStartFNode!(route, Inf)
	@assert(route.status == routeBeforeStartNode)

	# check that route functions return expected values
	if checkMode
		t = route.recentUpdateTime = 0.0
		travelModeIndex = 1
		@assert(isRouteUpToDate(route, t))
		@assert(getRouteCurrentLocation!(sim.net, route, t) == route.startLoc)
		@assert(getRouteNextNode!(sim, route, travelModeIndex, t)[1] == startFNode)
		@assert(getRouteNextNodeDist!(sim, route, t) == startFNodeDist)
		@assert(calcRouteDistance!(sim, route, t) == 0)
		route.recentUpdateTime = nullTime # reset
	end
end

# given a route and time, get current location
# Note that this may return the instance of `route.startLoc` or `route.endLoc` instead of a copy.
function getRouteCurrentLocation!(net::Network, route::Route, time::Float)
	updateRouteToTime!(net, route, time)

	fNodes = net.fGraph.nodes # shorthand
	if time <= route.startTime
		return route.startLoc

	elseif time >= route.endTime
		return route.endLoc

	elseif time <= route.startFNodeTime # and time > route.startTime
		# currently between startLoc and startFNode
		return linearInterpLocation(route.startLoc, fNodes[route.startFNode].location, route.startTime, route.startFNodeTime, time)

	elseif time >= route.endFNodeTime # and time < route.endTime
		# currently between endFNode and endLoc
		return linearInterpLocation(fNodes[route.endFNode].location, route.endLoc, route.endFNodeTime, route.endTime, time)

	else
		# currently somewhere on network
		return linearInterpLocation(fNodes[route.recentFNode].location, fNodes[route.nextFNode].location, route.recentFNodeTime, route.nextFNodeTime, time)
	end
end

# given route and current time,
# return next fNode index and time remaining to reach this node
# if already past last node in route, return that fNode and time to return to it
# testing: in certain cases, use given travel mode to determine travel time to the node
function getRouteNextNode!(sim::Simulation, route::Route, travelModeIndex::Int, time::Float)

	# shorthand:
	@unpack map, net = sim
	travelModes = sim.travel.modes

	# first need to update route
	updateRouteToTime!(net, route, time)

	# if route already finished, return nearest node
	# also return travel time to node based on route.endLoc
	if route.endTime <= time
		nearestFNode = route.endFNode
		travelTime = offRoadTravelTime(travelModes[travelModeIndex], map, route.endLoc, net.fGraph.nodes[nearestFNode].location)

		return nearestFNode, travelTime
	end

	nextFNode = nullIndex # init
	travelTime = nullTime # init
	if route.nextFNode != nullIndex
		# between startLoc and endFNode, go to nextFNode
		nextFNode = route.nextFNode
		travelTime = route.nextFNodeTime - time
	else
		# between endFNode and endLoc, return to recentFNode
		nextFNode = route.recentFNode
		travelTime = time - route.recentFNodeTime
	end

	# testing: changing travel time to match bartsim code...
	# scale travel time according to any change in travel mode
	if route.travelModeIndex != travelModeIndex
		if route.nextFNode == nullIndex
			# currently somewhere between route.endFNode and route.endLoc
			travelTime *= travelModes[route.travelModeIndex].offRoadSpeed / travelModes[travelModeIndex].offRoadSpeed
		end
	end

	if checkMode
		@assert(nextFNode != nullIndex)
		@assert(travelTime >= 0)
	end

	return nextFNode, travelTime
end

# Get the distance of the route to the next node at the given time.
# If already past last node in route, return the distance to return to it.
# See also: getRouteNextNode!
function getRouteNextNodeDist!(sim::Simulation, route::Route, time::Float)
	# shorthand:
	@unpack map, net = sim
	fGraph = net.fGraph

	# first need to update route
	updateRouteToTime!(net, route, time)

	if route.nextFNodeDist != nullDist return route.nextFNodeDist end

	dist = nullDist # init
	if route.status == routeAfterEndNode
		# off-road, after last node
		dist = normDist(map, route.endLoc, fGraph.nodes[route.endFNode].location)
		if time < route.endTime
			dist *= (time - route.endFNodeTime) / (route.endTime - route.endFNodeTime)
		end
	elseif route.status == routeBeforeStartNode
		# before first node, but may be off-road or partway along first arc
		dist = route.startFNodeDist # can be different from normDist between startLoc and startFNode if starting along arc
		if time > route.startTime
			dist *= (route.startFNodeTime - time) / (route.startFNodeTime - route.startTime)
		end
	else
		@assert(route.status == routeOnPath)
		arcIndex = fGraph.nodePairArcIndex[route.recentFNode, route.nextFNode]
		@assert(arcIndex != 0)
		arc = fGraph.arcs[arcIndex]
		dist = arc.distance * (route.nextFNodeTime - time) / (route.nextFNodeTime - route.recentFNodeTime)
	end
	@assert(dist >= 0)

	route.nextFNodeDist = dist

	return dist
end

function isRouteUpToDate(route::Route, time::Float)
	@assert(time != nullTime)
	@assert(route.recentUpdateTime <= time) # route should not be ahead of time

	# return time == route.recentUpdateTime # basic check, assuming that other fields of route are correct
	if time != route.recentUpdateTime return false end

	# further checks
	result = true
	if time >= route.endFNodeTime
		result &= (route.status == routeAfterEndNode)
		result &= (route.recentFNode == route.endFNode && route.nextFNode == nullIndex)
	elseif time <= route.startFNodeTime
		result &= (route.status == routeBeforeStartNode)
		result &= (route.recentFNode == nullIndex && route.nextFNode == route.startFNode)
	else # route.startFNodeTime <= time < route.endFNodeTime
		result &= (route.status == routeOnPath)
		result &= (route.recentFNode != nullIndex && route.nextFNode != nullIndex)
	end
	return result
end

# updates route fields for given time
function updateRouteToTime!(net::Network, route::Route, time::Float)
	@assert(time != nullTime)
	@assert(route.recentUpdateTime <= time)

	if isRouteUpToDate(route, time) return end # already up to date
	route.recentUpdateTime = time
	route.nextFNodeDist = nullDist # route is not up to date, so value should have changed

	if time <= route.startFNodeTime # time < route.startFNodeTime does not always work; see setRouteStateBeforeStartFNode!()
		@assert(route.status == routeBeforeStartNode && route.nextFNode == route.startFNode) # see setRouteStateBeforeStartFNode!()
	elseif time >= route.endFNodeTime
		setRouteStateAfterEndFNode!(route, time)
		# @assert(route.status == routeAfterEndNode && route.recentFNode == route.endFNode) # see setRouteStateAfterEndFNode!()
	else
		# currently somewhere on network
		updateRouteRecentRArc!(net, route, time)
		updateRouteRecentRArcFNode!(net, route, time)
		@assert(route.status == routeOnPath)
	end
end

# update route.recentRArc and similar fields for given time
# assumes that route is still on network
function updateRouteRecentRArc!(net::Network, route::Route, time::Float)

	if route.firstRArc == nullIndex
		return # do nothing
	end

	# should be on network
	@assert(route.startFNodeTime <= time < route.endFNodeTime)

	if route.status != routeOnPath
		# previously not on path, set to be at startFNode
		setRouteStateAfterStartFNode!(net, route, time)
	end

	if time < route.recentRArcEndTime # or: time < min(route.recentRArcEndTime, route.endFNodeTime)
		return # still on recentRArc, do nothing
	end

	# currently have: min(route.recentRArcEndTime, route.endFNodeTime) <= time
	# also know now that recentRArc is no longer current and so needs updating

	# shorthand:
	rNetTravel = net.rNetTravels[route.travelModeIndex]
	rNodeFNode = net.rNodeFNode

	@assert(route.startRNode != nullIndex != route.endRNode)

	if route.endRNodeTime <= time
		# on last rArc
		route.recentRArc = findRArcFromFNodeToFNode(net, rNodeFNode[route.endRNode], route.endFNode)
		route.recentRArcRecentFNode = 1 # remaining fNode data will be set in updateRouteRecentRArcFNode!()
		route.recentRArcStartTime = route.endRNodeTime
		route.recentRArcEndTime = route.endRNodeTime + rNetTravel.arcTimes[route.recentRArc]
		@assert(route.endFNodeTime < route.recentRArcEndTime) # should leave network before rArc ends
		return
	end

	# somewhere between startRNode and endRNode
	if route.startRNode != route.endRNode
		rArcs = net.rGraph.arcs # shorthand
		rArc = rArcs[route.recentRArc]
		recentRNode = rArc.fromNodeIndex
		nextRNode = rArc.toNodeIndex
		while route.recentRArcEndTime <= time && nextRNode != route.endRNode
			recentRNode = nextRNode

			rArcIndex = shortestPathNextRArc(net, route.travelModeIndex, nextRNode, route.endRNode)
			nextRNode = rArcs[rArcIndex].toNodeIndex

			route.recentRArc = rArcIndex
			route.recentRArcStartTime = route.recentRArcEndTime
			route.recentRArcEndTime += rNetTravel.arcTimes[route.recentRArc] # do not use route.recentRArcEndTime += rNetTravel.spTimes[recentRNode, nextRNode], as this may have low precision
		end
		if rNodeFNode[nextRNode] == route.endFNode
			@assert(isapprox(route.recentRArcEndTime, route.endFNodeTime))
			route.recentRArcEndTime = route.endFNodeTime # adjust value, to compensate for problems with numerical precision
		elseif nextRNode == route.endRNode
			@assert(isapprox(route.recentRArcEndTime, route.endRNodeTime))
			route.recentRArcEndTime = route.endRNodeTime # adjust value, to compensate for problems with numerical precision
		end
		route.recentRArcRecentFNode = 1 # remaining fNode data will be set in updateRouteRecentRArcFNode!()
	end

end

# for a given route and current time,
# updates 'recentRArcRecentFNode', 'recentRArcNextFNode', and other data in route
function updateRouteRecentRArcFNode!(net::Network, route::Route, time::Float)
	# assumes that route is currently still on the recent rArc

	if route.firstRArc == nullIndex
		return # do nothing
	end

	# check that recent rArc is still current
	@assert(route.recentFNodeTime != nullTime)
	@assert(route.recentFNodeTime <= time < min(route.recentRArcEndTime, route.endFNodeTime))

	if time < route.nextFNodeTime
		return # do nothing
	end

	# shorthand:
	fNetTravel = net.fNetTravels[route.travelModeIndex]
	rArcFNodes = net.rArcFNodes[route.recentRArc]
	recentRArcFNodesTimes = fNetTravel.rArcFNodesTimes[route.recentRArc]

	i = findMaxIndexLeqTime(recentRArcFNodesTimes, time - route.recentRArcStartTime, route.recentRArcRecentFNode)
	route.recentRArcRecentFNode = min(i, length(rArcFNodes) - 1) # adjust value, to compensate for problems with numerical precision in time values
	route.recentFNode = rArcFNodes[route.recentRArcRecentFNode]
	route.recentFNodeTime = route.recentRArcStartTime + recentRArcFNodesTimes[route.recentRArcRecentFNode]

	route.recentRArcNextFNode = route.recentRArcRecentFNode + 1
	route.nextFNode = rArcFNodes[route.recentRArcNextFNode]
	route.nextFNodeTime = route.recentRArcStartTime + recentRArcFNodesTimes[route.recentRArcNextFNode]
	if route.nextFNode == route.endFNode
		@assert(isapprox(route.nextFNodeTime, route.endFNodeTime))
		route.nextFNodeTime = route.endFNodeTime # adjust value, to compensate for problems with numerical precision
	end
end

# for a vector of strictly increasing time values,
# return the index of the largest value less than or equal (leq) to the given time
# minIndex is an optional lower bound on the index
function findMaxIndexLeqTime(times::Vector{Float}, time::Float, minIndex::Int = 1)

	if checkMode && false # skip this check, it is slow
		@assert(issorted(times, lt=<=)) # values should be strictly increasing
	end

	n = length(times)
	@assert(1 <= minIndex < n)
	if times[n] <= time
		return n
	end
	@assert(times[minIndex] <= time < times[n])

	# find range that contains given time
	# start at minIndex and step with increasing step sizes
	i = minIndex
	j = i + 1
	step = 2
	while times[j] <= time
		i = j
		j += step
		j >= n ? (j = n; break) : step *= 2
	end
	# should now have: times[i] <= time < times[j]

	# # or, instead of finding an upper bound for j, just set j to n
	# i = minIndex
	# j = n

	# binary search between i and j
	# maintain the property: times[i] <= time < times[j]
	while i < j-1
		k = div(i+j,2)
		times[k] <= time ? i = k : j = k
	end

	return i
end

# given a route with the start and end fNodes and rNodes already set,
# set the firstRArc of the route
function setRouteFirstRArc!(net::Network, route::Route)
	@assert(route.startFNode != nullIndex != route.endFNode)

	rNodeFNode = net.rNodeFNode # shorthand

	if route.startFNode == route.endFNode
		# no firstRArc needed
		route.firstRArc = nullIndex
		return

	elseif route.startRNode == nullIndex # and route.startFNode != route.endFNode
		@assert(route.endRNode == nullIndex)
		# find single rArc which contains path from startFNode to endFNode
		route.firstRArc = findRArcFromFNodeToFNode(net, route.startFNode, route.endFNode)

	elseif route.startFNode != rNodeFNode[route.startRNode] # and neither nodes are null
		# find arc from startFNode to startRNode
		route.firstRArc = findRArcFromFNodeToFNode(net, route.startFNode, rNodeFNode[route.startRNode])

	elseif route.startRNode == route.endRNode # and route.startFNode == rNodeFNode[route.startRNode]
		@assert(rNodeFNode[route.endRNode] != route.endFNode)
		# find arc from startRNode to endFNode
		route.firstRArc = findRArcFromFNodeToFNode(net, rNodeFNode[route.startRNode], route.endFNode)

	else # nullIndex != route.startRNode != route.endRNode != nullIndex
		# find first arc on path from startRNode to endRNode
		route.firstRArc = shortestPathNextRArc(net, route.travelModeIndex, route.startRNode, route.endRNode)
	end

	@assert(route.firstRArc != nullIndex)
end

# set temporally varying fields of route to represent state before reaching startFNode
function setRouteStateBeforeStartFNode!(route::Route, time::Float)
	@assert(time <= route.startFNodeTime) # @assert(time < route.startFNodeTime) does not work if starting directly at route.startFNode

	route.status = routeBeforeStartNode

	# recent rArc
	route.recentRArc = nullIndex
	route.recentRArcStartTime = nullTime
	route.recentRArcEndTime = nullTime

	# recent fNode
	route.recentRArcRecentFNode = nullIndex
	route.recentFNode = nullIndex
	route.recentFNodeTime = nullTime

	# next fNode
	route.recentRArcNextFNode = nullIndex
	route.nextFNode = route.startFNode
	route.nextFNodeTime = route.startFNodeTime
end

# set temporally varying fields of route to represent state just after leaving startFNode
function setRouteStateAfterStartFNode!(net::Network, route::Route, time::Float)
	@assert(route.startFNodeTime <= time < route.endFNodeTime)
	@assert(route.firstRArc != nullIndex)
	# or: if (route.firstRArc == nullIndex) setRouteStateBeforeStartFNode!(route, time); return; end

	# shorthand:
	rNetTravel = net.rNetTravels[route.travelModeIndex]
	fNetTravel = net.fNetTravels[route.travelModeIndex]
	rArc = net.rGraph.arcs[route.firstRArc]
	firstRArcFNodesTimes = fNetTravel.rArcFNodesTimes[route.firstRArc]
	firstRArcFNodes = net.rArcFNodes[route.firstRArc]

	route.status = routeOnPath

	# recent rArc
	route.recentRArc = route.firstRArc
	route.recentRArcStartTime = route.startFNodeTime - fNetTravel.fNodeFromRNodeTime[route.startFNode][rArc.fromNodeIndex]
	route.recentRArcEndTime = route.recentRArcStartTime + rNetTravel.arcTimes[route.firstRArc]
	if firstRArcFNodes[end] == route.endFNode
		@assert(isapprox(route.recentRArcEndTime, route.endFNodeTime))
		route.recentRArcEndTime = route.endFNodeTime # adjust value, to compensate for problems with numerical precision
	end

	# recent fNode
	route.recentRArcRecentFNode = net.rArcFNodeIndex[route.firstRArc][route.startFNode]
	route.recentFNode = route.startFNode
	route.recentFNodeTime = route.startFNodeTime
	if checkMode
		@assert(firstRArcFNodes[route.recentRArcRecentFNode] == route.startFNode)
	end

	# next fNode
	route.recentRArcNextFNode = route.recentRArcRecentFNode + 1
	route.nextFNode = firstRArcFNodes[route.recentRArcNextFNode]
	route.nextFNodeTime = route.recentRArcStartTime + firstRArcFNodesTimes[route.recentRArcNextFNode]
	if route.nextFNode == route.endFNode
		@assert(isapprox(route.nextFNodeTime, route.endFNodeTime))
		route.nextFNodeTime = route.endFNodeTime # adjust value, to compensate for problems with numerical precision
	end
end

# set temporally varying fields of route to represent state after leaving endFNode
function setRouteStateAfterEndFNode!(route::Route, time::Float)
	@assert(route.endFNodeTime <= time)

	route.status = routeAfterEndNode

	# recent rArc
	route.recentRArc = nullIndex
	route.recentRArcStartTime = nullTime
	route.recentRArcEndTime = nullTime

	# recent fNode
	route.recentRArcRecentFNode = nullIndex
	route.recentFNode = route.endFNode
	route.recentFNodeTime = route.endFNodeTime

	# next fNode
	route.recentRArcNextFNode = nullIndex
	route.nextFNode = nullIndex
	route.nextFNodeTime = nullTime
end

"""
	shortestRouteTravelTime!(sim::Simulation;
		startLoc::Location, firstNode::Int, dist1::Float, time1::Float, route::Route,
		endLoc::Location, lastNode::Int, dist2::Float, time2::Float,
		travelMode::TravelMode, travelPriority::Priority, startTime::Float)
Returns the travel time of the shortest route, given information about the start and end.

# Keyword arguments
Requires one of (in order of preference):
- `route` and `startTime` (use this if considering starting on `route`)
- `firstNode` (index of first node in new route) and `time1` (travel time from `startLoc` to `firstNode`)
- `firstNode` and `dist1` (distance between `startLoc` and `firstNode`)
- `startLoc` (can additionally pass in `firstNode` or `time1`)
Also requires one of (in order of preference):
- `lastNode` (index of last node in new route) and `time2` (travel time from `lastNode` to `endLoc`)
- `lastNode` and `dist2` (distance between `lastNode` and `endLoc`)
- `endLoc` (can additionally pass in `lastNode` or `time2`)
Also requires one of (in order of preference):
- `travelMode`
- `travelPriority` and `startTime`

Mutates: `sim.travel`, `route`

See also: [`shortestPathTravelTime`](@ref)
"""
function shortestRouteTravelTime!(sim::Simulation;
	route::Union{Route, Nothing} = nothing, startLoc::Union{Location, Nothing} = nothing, firstNode::Int = nullIndex, dist1::Float = nullDist, time1::Float = nullTime,
	endLoc::Union{Location, Nothing} = nothing, lastNode::Int = nullIndex, dist2::Float = nullDist, time2::Float = nullTime,
	travelMode::Union{TravelMode, Nothing} = nothing, travelPriority::Priority = nullPriority, startTime::Float = nullTime)

	# determine travel mode
	if travelMode == nothing
		@assert(startTime != nullTime)
		@assert(travelPriority != nullPriority)
		travelMode = getTravelMode!(sim.travel, travelPriority, startTime)
	end

	# find firstNode, and time to reach it from startLoc
	fNodes = sim.net.fGraph.nodes # shorthand
	if route == nothing
		if firstNode == nullIndex
			(firstNode, dist1) = findNearestNode(sim.map, sim.grid, fNodes, startLoc)
		end
		if time1 == nullTime
			if dist1 == nullDist
				@assert(startLoc != nothing)
				dist1 = normDist(sim.map, startLoc, fNodes[firstNode].location)
			end
			time1 = offRoadTravelTime(travelMode, dist1)
		end
	else
		@assert(startTime != nullTime)
		(firstNode, time1) = getRouteNextNode!(sim, route, travelMode.index, startTime)
	end

	# find lastNode, and time from it to endLoc
	if lastNode == nullIndex
		(lastNode, dist2) = findNearestNode(sim.map, sim.grid, fNodes, endLoc)
	end
	if time2 == nullTime
		if dist2 == nullDist
			@assert(endLoc != nothing)
			dist2 = normDist(sim.map, endLoc, fNodes[lastNode].location)
		end
		time2 = offRoadTravelTime(travelMode, dist2)
	end

	# path (on-road) travel time
	pathTravelTime = shortestPathTravelTime(sim.net, travelMode.index, firstNode, lastNode)

	# total travel time
	travelTime = time1 + pathTravelTime + time2

	return travelTime
end

# Return the distance travelled along a route,
# from the route start time to the given time.
function calcRouteDistance!(sim::Simulation, route::Route, time::Float)::Float
	@assert(route.status != routeNullStatus)

	# shorthand
	net = sim.net
	fGraph = net.fGraph

	# first need to update route
	updateRouteToTime!(net, route, time)

	dist = 0.0
	nextFNodeDist = getRouteNextNodeDist!(sim, route, time)
	if route.status == routeBeforeStartNode
		dist = route.startFNodeDist - nextFNodeDist
	elseif route.status == routeAfterEndNode
		dist += route.startFNodeDist
		dist += shortestPathDistance(net, route.travelModeIndex, route.startFNode, route.endFNode, startRNode = route.startRNode, endRNode = route.endRNode)
		dist += nextFNodeDist # distance to return back to route.endFNode
	else
		@assert(route.status == routeOnPath)
		dist += route.startFNodeDist
		dist += shortestPathDistance(net, route.travelModeIndex, route.startFNode, route.recentFNode)
		fArcIndex = fGraph.nodePairArcIndex[route.recentFNode, route.nextFNode]
		fArc = fGraph.arcs[fArcIndex]
		dist += fArc.distance - nextFNodeDist
	end
	@assert(dist >= 0)

	return dist
end


text:
module MieSeries
using StaticArrays
import SpecialFunctions
import AssociatedLegendrePolynomials as Legendre

export mieseries, incident_planewave

include("legendre.jl")
include("bessel.jl")
include("scattering.jl")

end # module

max_stars_repo_name: Riarrieta/MieSeries

max_issues_repo_name: Riarrieta/MieSeries

max_forks_repo_name: Riarrieta/MieSeries


text:
using GigaSOM, DataFrames, XLSX, CSV, Test, Random, Distributed, DistributedData
using FileIO, DataFrames, Distances
using JSON, SHA
import LinearAlgebra

owd = pwd()

"""
Check if the `pwd()` is the `/test` directory, and if not it changes to it.
"""
function checkDir()
    files = readdir()
    if !in("runtests.jl", files)
        cd(dirname(dirname(pathof(GigaSOM))))
    end
end

checkDir()

@testset "GigaSOM test suite" begin
    include("testDataOps.jl")
    include("testTrainutils.jl")
    include("testSplitting.jl")
    include("testInput.jl")

    #this loads the PBMC dataset required for the batch/parallel tests
    include("testLoadPBMC8.jl")
    include("testBatch.jl")
    include("testParallel.jl")

    #misc tests that require some of the above data too
    include("testInputCSV.jl")
    include("testFileSplitting.jl")
end

cd(owd)


text:
##############################################################################
#
# All Document types share a common metadata profile using DocumentMetadata
#
##############################################################################

type DocumentMetadata
    language::DataType
    name::Compat.UTF8String
    author::Compat.UTF8String
    timestamp::Compat.UTF8String
end
DocumentMetadata() = DocumentMetadata(
    EnglishLanguage,
    "Unnamed Document",
    "Unknown Author",
    "Unknown Time"
)

##############################################################################
#
# The abstract Document type
#
##############################################################################

abstract AbstractDocument

##############################################################################
#
# FileDocument type and constructors
#
##############################################################################

type FileDocument <: AbstractDocument
    filename::Compat.UTF8String
    metadata::DocumentMetadata
end

function FileDocument(f::AbstractString)
    d = FileDocument(Compat.UTF8String(f), DocumentMetadata())
    d.metadata.name = f
    return d
end

##############################################################################
#
# StringDocument type and constructors
#
##############################################################################

type StringDocument <: AbstractDocument
    text::Compat.UTF8String
    metadata::DocumentMetadata
end

StringDocument(txt::AbstractString) = StringDocument(Compat.UTF8String(txt), DocumentMetadata())

##############################################################################
#
# TokenDocument type and constructors
#
##############################################################################

type TokenDocument <: AbstractDocument
    tokens::Vector{Compat.UTF8String}
    metadata::DocumentMetadata
end
function TokenDocument(txt::AbstractString, dm::DocumentMetadata)
    TokenDocument(tokenize(dm.language, Compat.String(txt)), dm)
end
function TokenDocument{T <: AbstractString}(tkns::Vector{T})
    TokenDocument(tkns, DocumentMetadata())
end
TokenDocument(txt::AbstractString) = TokenDocument(Compat.UTF8String(txt), DocumentMetadata())

##############################################################################
#
# NGramDocument type and constructors
#
##############################################################################

type NGramDocument <: AbstractDocument
    ngrams::Dict{Compat.UTF8String,Int}
    n::Int
    metadata::DocumentMetadata
end
function NGramDocument(txt::AbstractString, dm::DocumentMetadata, n::Integer=1)
    NGramDocument(ngramize(dm.language, tokenize(dm.language, Compat.UTF8String(txt)), n),
        n, dm)
end
function NGramDocument(txt::AbstractString, n::Integer=1)
    NGramDocument(txt, DocumentMetadata(), n)
end
function NGramDocument{T <: AbstractString}(ng::Dict{T, Int}, n::Integer=1)
    NGramDocument(merge(Dict{AbstractString,Int}(), ng), n, DocumentMetadata())
end

##############################################################################
#
# text() / text!(): Access to document text as a string
#
##############################################################################

function text(fd::FileDocument)
    !isfile(fd.filename) && error("Can't find file: $(fd.filename)")
    readstring(fd.filename)
end

text(sd::StringDocument) = sd.text
function text(td::TokenDocument)
    warn("TokenDocument's can only approximate the original text")
    join(td.tokens, " ")
end
function text(ngd::NGramDocument)
    error("The text of an NGramDocument cannot be reconstructed")
end

text!(sd::StringDocument, new_text::AbstractString) = (sd.text = new_text)
function text!(d::AbstractDocument, new_text::AbstractString)
    error("The text of a $(typeof(d)) cannot be edited")
end

##############################################################################
#
# tokens() / tokens!(): Access to document text as a token array
#
##############################################################################

tokens(d::(@compat Union{FileDocument, StringDocument})) = tokenize(language(d), text(d))
tokens(d::TokenDocument) = d.tokens
function tokens(d::NGramDocument)
    error("The tokens of an NGramDocument cannot be reconstructed")
end

tokens!{T <: AbstractString}(d::TokenDocument, new_tokens::Vector{T}) = (d.tokens = new_tokens)
function tokens!{T <: AbstractString}(d::AbstractDocument, new_tokens::Vector{T})
    error("The tokens of a $(typeof(d)) cannot be directly edited")
end

##############################################################################
#
# ngrams() / ngrams!(): Access to document text as n-gram counts
#
##############################################################################

function ngrams(d::NGramDocument, n::Integer)
    error("The n-gram complexity of an NGramDocument cannot be increased")
end
ngrams(d::AbstractDocument, n::Integer) = ngramize(language(d), tokens(d), n)
ngrams(d::NGramDocument) = d.ngrams
ngrams(d::AbstractDocument) = ngrams(d, 1)

ngrams!(d::NGramDocument, new_ngrams::Dict{AbstractString, Int}) = (d.ngrams = new_ngrams)
function ngrams!(d::AbstractDocument, new_ngrams::Dict)
    error("The n-grams of $(typeof(d)) cannot be directly edited")
end

##############################################################################
#
# Length describes length of document in characters
#
##############################################################################

function Base.length(d::NGramDocument)
    error("NGramDocument's do not have a well-defined length")
end

Base.length(d::AbstractDocument) = length(text(d))

##############################################################################
#
# Length describes length of document in characters
#
##############################################################################

ngram_complexity(ngd::NGramDocument) = ngd.n

function ngram_complexity(d::AbstractDocument)
    error("$(typeof(d))'s have no n-gram complexity")
end

##############################################################################
#
# Union type that refers to a generic, non-abstract document type
#
##############################################################################

typealias GenericDocument @compat Union{
    FileDocument,
    StringDocument,
    TokenDocument,
    NGramDocument
}

##############################################################################
#
# Easier Document() constructor that decides types based on inputs
#
##############################################################################

Document(str::AbstractString) = isfile(str) ? FileDocument(str) : StringDocument(str)
Document{T <: AbstractString}(tkns::Vector{T}) = TokenDocument(tkns)
Document(ng::Dict{Compat.UTF8String, Int}) = NGramDocument(ng)

##############################################################################
#
# Conversion rules
#
##############################################################################

function Base.convert(::Type{StringDocument}, d::FileDocument)
    StringDocument(text(d), d.metadata)
end
function Base.convert(::Type{TokenDocument}, d::(@compat Union{FileDocument, StringDocument}))
    TokenDocument(tokens(d), d.metadata)
end
function Base.convert(::Type{NGramDocument},
            d::(@compat Union{FileDocument, StringDocument, TokenDocument}))
    NGramDocument(ngrams(d), 1, d.metadata)
end
Base.convert(::Type{TokenDocument}, d::TokenDocument) = d
Base.convert(::Type{NGramDocument}, d::NGramDocument) = d

##############################################################################
#
# getindex() methods: StringDocument("This is text and that is not")["is"]
#
##############################################################################

Base.getindex(d::AbstractDocument, term::AbstractString) = ngrams(d)[term]


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DiscretizeRelax/utilities/discretize_relax.jl
# Defines the structure used to compute bounds/relaxations via a discretize
# and relax approach.
#############################################################################

"""
DiscretizeRelax

An integrator for discretize and relaxation techniques.

$(TYPEDFIELDS)
"""
mutable struct DiscretizeRelax{M <: AbstractStateContractor, T <: Number, S <: Real, F, K, X, NY, JX, JP, INT, N} <: AbstractODERelaxIntegrator

    # Problem description
    "Initial Condition for pODEs"
    x0f::X
    "Jacobian w.r.t x"
    Jx!::JX
    "Jacobian w.r.t p"
    Jp!::JP
    "Parameter value for pODEs"
    p::Vector{Float64}
    "Lower Parameter Bounds for pODEs"
    pL::Vector{Float64}
    "Upper Parameter Bounds for pODEs"
    pU::Vector{Float64}
    "Number of state variables"
    nx::Int
    "Number of decision variables"
    np::Int
    "Time span to integrate over"
    tspan::Tuple{Float64, Float64}
    "Individual time points to evaluate"
    tsupports::Vector{Float64}

    # Options and internal storage
    "Maximum number of integration steps"
    step_limit::Int
    "Steps taken"
    step_count::Int
    "Stores solution X (from step 2) for each time"
    storage::Vector{Vector{T}}
    "Stores solution X (from step 1) for each time"
    storage_apriori::Vector{Vector{T}}
    "Stores each time t"
    time::Vector{Float64}
    "Support index to storage dictory"
    support_dict::Dict{Int,Int}
    "Holds data for numeric error encountered in integration step"
    error_code::TerminationStatusCode
    "Storage for bounds/relaxation of P"
    P::Vector{T}
    "Storage for bounds/relaxation of P - p"
    rP::Vector{T}
    "Relaxation Type"
    style::T
    "Flag indicating that only apriori bounds should be computed"
    skip_step2::Bool

    # Main functions used in routines
    "Functor for evaluating Taylor coefficients over a set"
    set_tf!::TaylorFunctor!{F,K,S,T}
    method_f!::M

    exist_result::ExistStorage{F,K,S,T}
    contractor_result::ContractorStorage{T}
    step_result::StepResult{T}
    step_params::StepParams

    new_decision_pnt::Bool
    new_decision_box::Bool

    relax_t_dict_indx::Dict{Int,Int}
    relax_t_dict_flt::Dict{Float64,Int}

    calculate_local_sensitivity::Bool
    local_problem_storage

    constant_state_bounds::Union{Nothing,ConstantStateBounds}
    polyhedral_constraint::Union{Nothing,PolyhedralConstraint}
    prob
end

function DiscretizeRelax(d::ODERelaxProb, m::SCN; repeat_limit = 50, step_limit = 1000, tol = 1E-4, hmin = 1E-13,
                         relax = false, h = 0.0, skip_step2 = false, Jx! = nothing, Jp! = nothing) where SCN <: AbstractStateContractorName

    γ = state_contractor_γ(m)::Float64
    k = state_contractor_k(m)::Int
    method_steps = state_contractor_steps(m)::Int

    tsupports = d.support_set.s
    support_dict = Dict{Int,Int}()
    error_code = RELAXATION_NOT_CALLED

    T = relax ? MC{d.np,NS} : Interval{Float64}
    style = zero(T)
    time = zeros(Float64, 1000)
    storage = Vector{T}[]
    storage_apriori = Vector{T}[]
    for i = 1:1000
        push!(storage, zeros(T, d.nx))
        push!(storage_apriori, zeros(T, d.nx))
    end
    P = zeros(T, d.np)
    rP = zeros(T, d.np)

    Δ = FixedCircularBuffer{Vector{T}}(method_steps)
    A_Q = FixedCircularBuffer{Matrix{Float64}}(method_steps)
    A_inv = FixedCircularBuffer{Matrix{Float64}}(method_steps)
    for i = 1:method_steps
        push!(Δ, zeros(T, d.nx))
        push!(A_Q, Float64.(Matrix(I, d.nx, d.nx)))
        push!(A_inv, Float64.(Matrix(I, d.nx, d.nx)))
    end

    state_method = state_contractor(m, d.f, Jx!, Jp!, d.nx, d.np, style, zero(Float64), h)
    is_adaptive = (h <= 0.0)

    set_tf! = TaylorFunctor!(d.f, d.nx, d.np, Val(k), style, zero(Float64))
    exist_storage = ExistStorage(set_tf!, style, P, d.nx, d.np, k, h, method_steps)
    contractor_result = ContractorStorage(style, d.nx, d.np, k, h, method_steps)
    contractor_result.γ = γ
    contractor_result.P = P
    contractor_result.rP = rP

    contractor_result.is_adaptive = is_adaptive
    step_params = StepParams(tol, hmin, repeat_limit, is_adaptive, skip_step2)
    step_result = StepResult{typeof(style)}(zeros(d.nx), zeros(typeof(style), d.nx),
                                            A_Q, A_inv, Δ, 0.0, 0.0)

    relax_t_dict_indx = Dict{Int,Int}()
    relax_t_dict_flt = Dict{Float64,Int}()

    calculate_local_sensitivity = false
    local_integrator = state_contractor_integrator(m)
    local_problem_storage = ODELocalIntegrator(d, local_integrator)

    return DiscretizeRelax{typeof(state_method), T, Float64, typeof(d.f), k+1,
                           typeof(d.x0), d.nx+d.np, typeof(Jx!), typeof(Jp!),
                           typeof(local_integrator), d.np}(d.x0, Jx!, Jp!, d.p,
                           d.pL, d.pU, d.nx, d.np, d.tspan, tsupports,
                           step_limit, 0, storage, storage_apriori, time,
                           support_dict, error_code, P, rP, skip_step2, style,
                           set_tf!, state_method, exist_storage, contractor_result,
                           step_result, step_params, true, true,
                           relax_t_dict_indx, relax_t_dict_flt,calculate_local_sensitivity,
                           local_problem_storage, d.constant_state_bounds,
                           d.polyhedral_constraint, d)
end

function DiscretizeRelax(d::ODERelaxProb; kwargs...)
    DiscretizeRelax(d, LohnerContractor{4}(); kwargs...)
end

"""
set_P!(d::DiscretizeRelax)

Initializes the `P` and `rP` (P - p) fields of `d`.
"""
function set_P!(d::DiscretizeRelax{M,Interval{Float64},S,F,K,X,NY}) where {M<:AbstractStateContractor, S, F, K, X, NY}

    @__dot__ d.P = Interval(d.pL, d.pU)
    @__dot__ d.rP = d.P - d.p

    return nothing
end

function set_P!(d::DiscretizeRelax{M,MC{N,T},S,F,K,X,NY}) where {M<:AbstractStateContractor, T<:RelaxTag, S <: Real, F, K, X, N, NY}

    @__dot__ d.P = MC{N,NS}.(d.p, Interval(d.pL, d.pU), 1:d.np)
    @__dot__ d.rP = d.P - d.p

    return nothing
end

"""
compute_X0!(d::DiscretizeRelax)

Initializes the circular buffer that holds `Δ` with the `out - mid(out)` at
index 1 and a zero vector at all other indices.
"""
function compute_X0!(d::DiscretizeRelax)

    d.storage[1] .= d.x0f(d.P)
    d.storage_apriori[1] .= d.storage[1]

    if in(d.tspan[1], d.tsupports)
        d.relax_t_dict_indx[1] = 1
        d.relax_t_dict_flt[d.tspan[1]] = 1
    end

    d.step_result.Xⱼ .= d.storage[1]
    d.step_result.xⱼ .= mid.(d.step_result.Xⱼ)

    d.exist_result.Xj_0 .= d.step_result.Xⱼ
    d.exist_result.predicted_hj = d.tspan[2] - d.tspan[1]

    d.contractor_result.Xj_0 .= d.step_result.Xⱼ
    d.contractor_result.xval .= mid.(d.exist_result.Xj_0)
    d.contractor_result.pval .= d.p

    return nothing
end

"""
set_Δ!

Initializes the circular buffer, `Δ`, that holds `Δ_i` with the `out - mid(out)` at
index 1 and a zero vector at all other indices.
"""
function set_Δ!(Δ::FixedCircularBuffer{Vector{T}}, out::Vector{Vector{T}}) where T

    Δ[1] .= out[1] .- mid.(out[1])
    for i = 2:length(Δ)
        fill!(Δ[i], zero(T))
    end

    return nothing
end


text:
"""
    InterpolationLineSearch{C, T}

`InterpolationLineSearch` uses quadratic and cubic interpolations to compute the step length
satisfying the Armijo sufficient decrease condition.

# Fields
- `c::C`: the constant for the Armijo condition in the open interval (0, 1)
- `ϵ::T`: the minimum step length decrease from the initial step length
"""
struct InterpolationLineSearch{C, T}
    c::C
    ϵ::T

    function InterpolationLineSearch(c, ϵ)
        0 < c < 1 || error("`c` must be in the open interval (0, 1)")
        return new{typeof(c), typeof(ϵ)}(c, ϵ)
    end
end

"""
    InterpolationLineSearch(; kwargs...)

Initiate `InterpolationLineSearch`.

# Keywords
- `c=1e-4`: the constant for the Armijo condition in the open interval (0, 1)
- `ϵ=1e-6`: the minimum step length decrease from the initial step length
"""
InterpolationLineSearch(; c=1e-4, epsilon=1e-6) = InterpolationLineSearch(c, epsilon)

function (ils::InterpolationLineSearch)(f, state, p, α₀::T) where {T}
    x = state.x
    ϕ₀ = state.f
    ∇y = state.∇f
    dϕ₀ = ∇y⋅p
    αₖ = α₀
    c = ils.c
    ϵ = ils.ϵ
    ϕ(α) = f(x + α * p)
    ϕαₖ = ϕ(αₖ)
    αₖ₊₁ = T(-(dϕ₀ * αₖ^2) / (2 * (ϕαₖ - ϕ₀ - dϕ₀ * αₖ)))
    while f(x + αₖ₊₁*p) > ϕ₀ + c*αₖ₊₁*dϕ₀
        ϕαₖ₊₁ = ϕ(αₖ₊₁)
        a, b = [αₖ^2 -αₖ₊₁^2; -αₖ^3 αₖ₊₁^3] *
            [ϕαₖ₊₁ - ϕ₀ - dϕ₀*αₖ₊₁; ϕαₖ - ϕ₀ - dϕ₀*αₖ] /
            (αₖ^2*αₖ₊₁^2*(αₖ₊₁-αₖ))
        αₖ₊₂ = (-b + sqrt(b^2 - 3a*dϕ₀)) / (3a)
        αₖ₊₂ = max(min(αₖ₊₂, αₖ₊₁ - ϵ), αₖ₊₁/2)
        αₖ = αₖ₊₁
        ϕαₖ = ϕαₖ₊₁
        αₖ₊₁ = T(αₖ₊₂)
    end
    return αₖ₊₁
end


text:
# This code should perhaps be moved into a separate package

specializations(f::Function) =
    MethodInstance[spec for m in methods(f).ms for spec in specializations(m)]
specializations(m::Method) = gather_specializations(m.specializations)
gather_specializations(::Void) = MethodInstance[]  # needs special-casing...
function gather_specializations(specs::Union{Base.TypeMapEntry, Base.TypeMapLevel})
    res = MethodInstance[]
    Base.visit(specs) do spec
        push!(res, spec)
    end
    res
end

is_compiled(spec::MethodInstance) = spec.fptr!=C_NULL

"""    number_of_specializations(obj::Module/Function)

For each function/method in `obj`, return the number of specializations, and the
number of _compiled_ specializations (excluding those that were merely inferred). """
function number_of_specializations(mod::Module)
    res = []
    for name in names(mod, true) # see ?names
        r = @ignore_errors nothing eval(mod, name)
        if r isa Function
            specs = specializations(r)
            if length(specs) > 0
                push!(res, (r, length(specs), count(is_compiled, specs)))
            end
        end
    end
    sort(res, by=last, rev=true)
end
number_of_specializations(fun::Function) =
    sort([(specs = specializations(m);
           (m, length(specs), count(is_compiled, specs)))
          for m in methods(fun)], by=last, rev=true)

is_inferred(spec::MethodInstance) =
    isleaftype(spec.rettype)   # seems reasonable, but not 100% sure


text:
using Documenter, HTTP

"""
Loops through the files in the examples folder and adds them (with any header comments) to the examples.md
markdown file.
"""
function generateExamples()
    f = open(joinpath(@__DIR__, "src/examples.md"), "w")
    write(
        f,
        "```@meta
    # NOTE this file is autogenerated, do not edit examples.md directly. To make an example, upload the .jl file to the examples folder. Header comments may be included at the top of the file using \"\"\" syntax
``` ",
    )
    write(f, "\n")
    write(f, "# Examples")
    write(
        f,
        "\nSome examples that may prove potentially useful for those using
`HTTP.jl`. The code for these examples can also be found on Github
 in the `docs/examples` folder.",
    )
    for (root, dirs, files) in walkdir(joinpath(@__DIR__, "examples"))
        #set order of files so simple is first, and Readme examples are last
        temp = files[1]
        files[findfirst(isequal("simple_server.jl"), files)] = temp
        files[1] = "simple_server.jl"
        temp = files[2]
        files[findfirst(isequal("cors_server.jl"), files)] = temp
        files[2] = "cors_server.jl"
        temp = files[length(files)]
        files[findfirst(isequal("readme_examples.jl"), files)] = temp
        files[length(files)] = "readme_examples.jl"

        for file in files
            println(file)
            #extract title from example
            write(f, "\n")
            title = file
            title = replace(title, "_" => " ")
            title = replace(title, ".jl" => "")
            title = titlecase(title)
            title = "## " * title * "\n"
            write(f, title)
            #open each file and read contents
            opened = open(joinpath(@__DIR__, "examples/") * file)
            lines = readlines(opened, keep = true)
            index = 1
            #find doc string intro if exists
            if "\"\"\"\n" in lines
                index = findall(isequal("\"\"\"\n"), lines)[2]
                print(index)
                for i = 2:index-1
                    write(f, lines[i])
                end
                lines = lines[index+1:end]
            end

            write(f, "```julia")
            write(f, "\n")
            for line in lines
                write(f, line)
            end
            write(f, "\n")
            write(f, "```")
            close(opened)
        end
    end
    close(f)
end

generateExamples()

makedocs(
    modules = [HTTP],
    sitename = "HTTP.jl",
    pages = [
        "Home" => "index.md",
        "public_interface.md",
        "internal_architecture.md",
        "internal_interface.md",
        "examples.md",
    ],
)

deploydocs(repo = "github.com/JuliaWeb/HTTP.jl.git")


text:
using Test
using SetProg
const MOIT = MOI.Test

function bridged_mock(mock_optimize!::Function...;
                      # We use a UF to add support for SOSCone, ... so that
                      # we don't have to set the variables created by the SOS
                      # bridges
                      model = MOIU.UniversalFallback(MOI.Utilities.Model{Float64}()))
    mock = MOI.Utilities.MockOptimizer(model)
    bridged = MOI.Bridges.full_bridge_optimizer(mock, Float64)
    MOI.Utilities.set_mock_optimize!(mock, mock_optimize!...)
    return bridged
end


text:
# Tests for helper functions in src/high_level.jl.
# Tests for get_* functions are in test/get_*.jl.

using SeisRequests, Test

import StationXML

@testset "High-level" begin
    @testset "filter_xml" begin
        sxml = StationXML.FDSNStationXML(source="Me", created=DateTime("2000-01-01"),
        schema_version="1.1")
        net1 = StationXML.Network(code="AN")
        net2 = StationXML.Network(code="BO")
        sta1 = StationXML.Station(code="FAKE", longitude=1, latitude=2, elevation=3,
            site=StationXML.Site(name="Fake site"))
        sta2 = StationXML.Station(code="FAKE2", longitude=2, latitude=3, elevation=4,
            site=StationXML.Site(name="Fake site 2"))
        cha1 = StationXML.Channel(code="UHT", longitude=10, latitude=20, elevation=30,
            depth=40, location_code="")
        cha2 = StationXML.Channel(code="XYZ", longitude=20, latitude=30, elevation=40,
            depth=50, location_code="00")

        # sxml′ is the filtered version, sxml″ is the correct one
        @testset "Network" begin
            sxml′ = deepcopy(sxml)
            sxml″ = deepcopy(sxml)
            append!(sxml′.network, deepcopy([net1, net2]))
            push!(sxml″.network, deepcopy(net1))
            @test SeisRequests.filter_stationxml(sxml′, net1) == sxml″
        end

        @testset "Station" begin
            sxml′ = deepcopy(sxml)
            sxml″ = deepcopy(sxml)
            append!(sxml′.network, deepcopy([net1, net2]))
            for n in sxml′.network
                append!(n.station, deepcopy([sta1, sta2]))
            end
            push!(sxml″.network, deepcopy(net1))
            push!(sxml″.network[1].station, deepcopy(sta1))
            @test SeisRequests.filter_stationxml(sxml′, net1, sta1) == sxml″
        end

        @testset "Channel" begin
            sxml′ = deepcopy(sxml)
            sxml″ = deepcopy(sxml)
            append!(sxml′.network, deepcopy([net1, net2]))
            for n in sxml′.network
                append!(n.station, deepcopy([sta1, sta2]))
                for s in n.station
                    append!(s.channel, deepcopy([cha1, cha2]))
                end
            end
            push!(sxml″.network, deepcopy(net1))
            push!(sxml″.network[1].station, deepcopy(sta1))
            push!(sxml″.network[1].station[1].channel, deepcopy(cha1))
            @test SeisRequests.filter_stationxml(sxml′, net1, sta1, cha1) == sxml″
        end
    end

    @testset "_getifnotmissing" begin
        @test SeisRequests._getifnotmissing((a=1, b=missing), :a) == 1
        @test SeisRequests._getifnotmissing(missing, :a) === missing
    end
end


text:
module Losses

using Statistics
import Knet
using Photon: ϵ
import Photon.Metrics: Metric

export  Loss, L1Loss, MAELoss, L2Loss, MSELoss, LNLoss, PseudoHuberLoss, HingeLoss,
        SquaredHingeLoss, BCELoss, CrossEntropyLoss, FocalLoss

# Simple wrapper since Knet nll loss function expects regular array.
# Knet.nll(y, a::Knet.KnetArray{<:Integer}; dims=1, average=true) = nll(y, Array(a); dims=dims, average=average)

"""
Abstract type for the loss functions. However Photon accepts
any function as a loss function as long as it is callable and
returns the loss value as a scalar type.
"""
abstract type Loss <: Metric end

"""
Calculates the mean absolute error between `label` and `pred`.

``L = \\sum_i \\vert {label}_i - {pred}_i \\vert``

`label` and `pred` can have arbitrary shape as long as they have the same
number of elements.
"""
struct L1Loss <: Loss
  reduction

  L1Loss(;reduction=mean) = new(reduction)
end

function (l::L1Loss)(ŷ, y)
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)
  l.reduction(abs.(ŷ - y))
end

MAELoss = L1Loss

"""
Calculates the mean squared error between `label` and `pred`.

``L = \\frac{1}{2} \\sum_i ({label}_i - {pred}_i)^2``

`label` and `pred` can have arbitrary shape as long as they have the same
number of elements.
"""
struct L2Loss <: Loss
  reduction

  L2Loss(;reduction=mean) = new(reduction)
end

function (l::L2Loss)(ŷ, y)
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)
  l.reduction((ŷ - y).^2)
end

MSELoss = L2Loss

"""
More generic version of L1loss and L2loss that allows L{N}Loss.

``L = \\frac{1}{2} \\sum_i \\vert {label}_i - {pred}_i \\vert^{N}``

`label` and `pred` can have arbitrary shape as long as they have the same
number of elements.

"""
struct LNLoss <: Loss
  n
  reduction

  LNLoss(n::Int;reduction=mean) = new(n, reduction)
end

function (l::LNLoss)(ŷ, y)
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)
  l.reduction(abs.(ŷ - y).^l.n)
end

"""
Pseudo Huber Loss implementation, somewhere between a L1 and L2 loss.
"""
struct PseudoHuberLoss <: Loss
  delta
  reduction
  PseudoHuberLoss(delta=1; reduction=mean) = new(delta, reduction)
end

function (l::PseudoHuberLoss)(ŷ, y)
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)
  δ = l.delta
  r = δ^2 .* (sqrt.(1 .+ ((ŷ - y)/δ).^2) .- 1)
  l.reduction(r)
end

"""
BCE (Binary Cross Entropy) loss is useful when training logistic regression.

``L = - \\sum_i {label}_i * \\log({pred}_i) +
            (1 - {label}_i) * \\log(1 - {pred}_i)``

If `use_sigmoid` is true, first a sigmoid activation will be applied before
calcuating the BCE loss.
"""
struct BCELoss <: Loss
  reduction
  use_sigmoid
  BCELoss(;reduction=mean, use_sigmoid=false) = new(reduction, use_sigmoid)
end

function (l::BCELoss)(ŷ, y)
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)

  ŷ = l.use_sigmoid ? Knet.sigm(ŷ) : ŷ

  r = -y .* log.(ŷ .+ ϵ) .- (1 .- y) .* log.(1 .- ŷ .+ ϵ)
  l.reduction(r)
end

"""
CrossEntropy loss function with support for an optional weight parameter.
The weight parameter can be static (for example to handle class inbalances)
or dynamic, so passed every time when the lost function is invoked.

Cross Entropy Function:

``L = - \\sum_i {label}_i \\log({output}_i)``

If use_softmax is true, first the softmax(output) will be performed before
cross entropy.

# Usage

```julia
workout = Workout(model, CrossEntropyLoss(), SGD())
```

"""
struct CrossEntropyLoss <: Loss
  weight
  reduction
  use_softmax
  CrossEntropyLoss(weight=1; reduction=mean, use_softmax=false) =
    new(weight, reduction, use_softmax)
end

function (l::CrossEntropyLoss)(ŷ, y, weight=nothing)
  w = weight !== nothing ? weight : l.weight
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)

  ŷ = l.use_softmax ? Knet.softmax(ŷ) : ŷ

  r = sum(y .* log.(ŷ .+ ϵ), dims=1) .* w
  l.reduction(-r)
end

"""
Calculates the hinge loss function often used in SVMs:

``L = \\sum_i max(0, {margin} - {pred}_i \\cdot {label}_i)``

where `pred` is the classifier prediction and `label` is the target tensor
containing values -1 or 1. `label` and `pred` must have the same number of
elements.

If autofix is true, will convert label from {0,1} to {-1,1}

"""
struct HingeLoss <: Loss
  reduction
  margin
  autofix
  HingeLoss(;margin=1.0, reduction=mean, autofix=true) = new(reduction,margin, autofix)
end

function (l::HingeLoss)(ŷ, y)
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)

  # Hinge lost requires y to be {-1,1}. So the following converts {0,1} to
  # the right values
  if l.autofix
    y[Array(y) .< 1.0] .= -1.0
  end

  r = max.(0, l.margin .- (ŷ .* y))
  l.reduction(r)
end

"""
Calculates the soft-margin loss function used in SVMs:

``L = \\sum_i max(0, {margin} - {pred}_i \\cdot {label}_i)^2``

where `pred` is the classifier prediction and `label` is the target tensor
containing values -1 or 1. `label` and `pred` can have arbitrary shape as
long as they have the same number of elements.

"""
struct SquaredHingeLoss <: Loss
  reduction
  margin
  autofix
  SquaredHingeLoss(;margin=1.0, reduction=mean, autofix=true) = new(reduction,margin, autofix)
end

function (l::SquaredHingeLoss)(ŷ, y)
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)

  # Hinge lost requires y to be {-1,1}. So the following converts {0,1} to
  # the right values
  if l.autofix
    y[Array(y) .< 1.0] .= -1.0
  end

  r = max.(0, l.margin .- (ŷ .* y)).^2
  l.reduction(r)
end

"""
Focal loss for multi-classification

``L=-alpha(1-p_t)^{gamma}ln(p_t)``

See also Focal Loss for Dense Object Detection
    https://arxiv.org/abs/1708.02002
"""
struct FocalLoss <: Loss
  reduction
  gamma
  alpha
  FocalLoss(;gamma=2, alpha=2, reduction=mean) = new(reduction, gamma, alpha)
end

function (l::FocalLoss)(ŷ, y)
  ŷ, y = Knet.mat(ŷ), Knet.mat(y)
  γ = l.gamma

  loss =   @.      y  * ( - (1 - ŷ)^γ  * log(    ŷ + ϵ)) # if y = 1
  loss .+= @. (1 - y) * ( - (    ŷ)^γ  * log(1 - ŷ + ϵ)) # if y = 0

  l.reduction(loss)
end

end


text:
##############################################################################
##
## DynamicSolution
##
##############################################################################

type DynamicAiyagariSolution{T}
    as::AiyagariSolution{T}
    G1::Matrix{T}
    impact::Matrix{T}
    eu::Vector{Int}
end

##############################################################################
##
## F
##
##############################################################################

function F(ap::AiyagariProblem, ρπ, σπ, V, g, K, r, w, π, Vpost, gpost, Kpost, rpost, wpost, πpost, Verrors, ηπ)
    newg = vcat(one(eltype(g)) - sum(g), g)
    as = AiyagariSolution(V, newg, K, r, w)
    aa = AiyagariArrays(ap, as)
    update_Au!(ap, aa, as)
    update_B!(ap, aa, as)
    hjbResidual = aa.u - aa.B' * V .+ Vpost .- V .- Verrors
    newgIntermediate = aa.A * as.g
    gResidual = gpost .- g -  newgIntermediate[2:end]
    KResidual = K - sum_capital(ap.a, newg)
    rResidual = exp(π) * ap.α * K^(ap.α-1) - ap.δ - r
    wResidual = exp(π) * (1-ap.α) * K^ap.α  - w
    πResidual = πpost - π - (-(1-ρπ) * π + σπ * ηπ)
    return hjbResidual, gResidual, KResidual, rResidual, wResidual, πResidual
end

function F(ap, ρπ, σπ, x)
    lV = length(ap.a) * length(ap.z)
    lg = length(ap.a) * length(ap.z) - 1
    start = 1
    V = x[start:(start + lV - 1)]
    start += lV
    g = x[start:(start + lg - 1)]
    start += lg
    K = x[start]
    start += 1
    r = x[start]
    start += 1
    w = x[start]
    start += 1
    π = x[start]
    start += 1
    Vpost = x[start:(start + lV - 1)]
    start += lV
    gpost = x[start:(start + lg - 1)]
    start += lg
    Kpost = x[start]
    start += 1
    rpost = x[start]
    start += 1
    wpost = x[start]
    start += 1
    πpost = x[start]
    start += 1
    Verrors = x[start:(start + lV - 1)]
    start += lV
    ηπ = x[start]
    @assert start == length(x)
    return vcat(F(ap, ρπ, σπ, V, g, K, r, w, π, Vpost, gpost, Kpost, rpost, wpost, πpost, Verrors, ηπ)...)
end

function solve(ap::AiyagariProblem, ρπ, σπ)
    println("Compute Steady State")
    as = solve(ap, verbose = false)
    g = as.g[2:end]
    x0 = vcat(as.V, g, as.K, as.r, as.w, ap.π, deepcopy(as.V), deepcopy(g), as.K, as.r, as.w, ap.π, zeros(as.V), 0.0)
    ap.invΔ = 0.0

    # check function gives 0.0 at steady state
    @assert maxabs(F(ap, ρπ, σπ, x0)) < 1e-4

    # compute jacobian around 0
    println("Compute Jacobian")
    out = Calculus.jacobian(x -> F(ap, ρπ, σπ, x))(x0)

    # unpack derivatives
    Γ1 = - out[:, 1:(length(as.V) + length(g) + 4)]
    Γ0 = out[:, (length(as.V) + length(g) + 4 + 1):(2 * length(as.V) + 2 * length(g) + 2 * 4)]
    Π = - out[:, (2 * length(as.V) + 2 * length(g) + 2 * 4 + 1):(2 * length(as.V) + 2 * length(g) + 2 * 4 + length(as.V))]
    Ψ = - out[:, end:end]

    # solve system
    println("Solve Linearized Model")
    c = zeros(length(as.V) + length(g) + 4, 1)
    G1, C, impact, fmat, fwt, ywt, gev, eu, loose = gensysdt(Γ0, Γ1, c, Ψ, Π)
    return DynamicAiyagariSolution(as, G1, impact, eu)
end

function simulate(ap::AiyagariProblem, as::DynamicAiyagariSolution)
    T = 200
    N = 2000
    T1 = div(N, T)
    dt = T / N
    lV = length(ap.a) * length(ap.z)
    lg = length(ap.a) * length(ap.z) - 1
    g0 = as.as.g[2:end]
    x = zeros(Float64, lV + lg + 4)
    V = Array(Float64, lV, N)
    g = Array(Float64, lg, N)
    K = Array(Float64, N)
    r = Array(Float64, N)
    w = Array(Float64, N)
    A = as.G1 * dt
    B = as.impact * sqrt(dt)
    for n in 1:N
        # unpack solution
        start = 1
        V[:, n] = x[start:(start + lV - 1)]
        start += lV
        g[:, n] = x[start:(start + lg - 1)]
        start += lg
        K[n] = x[start]
        start += 1
        r[n] = x[start]
        start += 1
        w[n] = x[start]
        if n <= T1
            newx = A * x + B
        else
            newx = A * x
        end

        # update
        newx, x = x, newx
    end
    return linspace(1, T, N), V, g, K/as.as.K, r/as.as.r, w/as.as.w
end


text:
using OrdinaryDiffEq, OpenQuantumTools, Plots

coupling_1 = ConstantCouplings(["X"])
bath_1 = Ohmic(1e-4, 4, 16)
interaction_1 = Interaction(coupling_1, bath_1)

coupling_2 = ConstantCouplings(["Z"])
bath_2 = Ohmic(1e-4, 0.1, 16)
interaction_2 = Interaction(coupling_2, bath_2)

interaction_set = InteractionSet(interaction_1, interaction_2);

H = DenseHamiltonian([(s) -> 1.0], -[σz], unit = :ħ)
u0 = PauliVec[1][1]
annealing_1 = Annealing(H, u0, coupling=coupling_1, bath = bath_1)
annealing_2 = Annealing(H, u0, coupling=coupling_2, bath = bath_2)
annealing = Annealing(H, u0, interactions=interaction_set)

tf = 10
# Generate the unitary first
U = solve_unitary(annealing, tf, alg=Tsit5(), reltol=1e-6)
# tag the unitary so the solver know it has inplace update method
# this will speed up the calculation of integral
U = InplaceUnitary(U)
 # Solve the Redfield equation
sol_1 = solve_redfield(annealing_1, tf, U, alg = Tsit5(), abstol = 1e-6, reltol = 1e-6)
sol_2 = solve_redfield(annealing_2, tf, U, alg = Tsit5(), abstol = 1e-6, reltol = 1e-6)
sol = solve_redfield(annealing, tf, U, alg = Tsit5(), abstol = 1e-6, reltol = 1e-6)

t_list = range(0,tf,length=200)
x_nopulse = []
z_nopulse = []
xz_nopulse = []
for s in t_list
    push!(x_nopulse, real(tr(σx*sol_1(s))))
    push!(z_nopulse, real(tr(σx*sol_2(s))))
    push!(xz_nopulse, real(tr(σx*sol(s))))
end
plot(t_list, x_nopulse, linewidth=2, label="X coupling")
plot!(t_list, z_nopulse, linewidth=2, label="Z coupling")
plot!(t_list, xz_nopulse, linewidth=2, label="X+Z coupling")
xlabel!("s")
ylabel!("<X>")

# in this example, we apply a Z pulse in the middle of the evolution
# for the InstPulseCallback constructor
# the first argument is a list of times where the pulses are applied
# the second argument is a function to update the state update!(c, pulse_index
# the function will update the state c with give pulse_index
cbu = InstPulseCallback([0.5 * tf], (c, x) -> c .= σx * c)
cb = InstPulseCallback([0.5 * tf], (c, x) -> c .= σx * c * σx)
annealing_1 = Annealing(H, u0, coupling = coupling_1, bath = bath_1)
annealing_2 = Annealing(H, u0, coupling=coupling_2, bath = bath_2)
annealing = Annealing(H, u0, interactions=interaction_set)

tf = 10
U = solve_unitary(annealing, tf, alg=Tsit5(), reltol=1e-6, callback = cbu);
U = InplaceUnitary(U)

sol_1 = solve_redfield(annealing_1, tf, U, alg = Tsit5(), reltol = 1e-6, callback=cb)
sol_2 = solve_redfield(annealing_2, tf, U, alg = Tsit5(), reltol = 1e-6, callback=cb)
sol = solve_redfield(annealing, tf, U, alg = Tsit5(), reltol = 1e-6, callback=cb);
t_list = range(0,tf,length=200)
x_pulse = []
z_pulse = []
xz_pulse = []
for s in t_list
    push!(x_pulse, real(tr(σx*sol_1(s))))
    push!(z_pulse, real(tr(σx*sol_2(s))))
    push!(xz_pulse, real(tr(σx*sol(s))))
end
plot(t_list, x_pulse, linewidth=2, label="X coupling", legend=:topleft, color=1)
plot!(t_list, z_pulse, linewidth=2, label="Z coupling", color=2)
plot!(t_list, xz_pulse, linewidth=2, label="X+Z coupling", color=3)
plot!(t_list, x_nopulse, linewidth=2, linestyle=:dash, label="X coupling with no pulse", color=1)
plot!(t_list, z_nopulse, linewidth=2, linestyle=:dash, label="Z coupling with no pulse", color=2)
plot!(t_list, xz_nopulse, linewidth=2, linestyle=:dash, label="X+Z coupling with no pulse", color=3)
xlabel!("s")
ylabel!("<X>")


text:
using LinearAlgebra
using MixedModels
using Random
using SparseArrays
using Test

using MixedModels: allequal, average, densify, dataset

@testset "average" begin
	@test average(1.1, 1.2) == 1.15
end

@testset "densify" begin
	@test densify(sparse(1:5, 1:5, ones(5))) == Diagonal(ones(5))
	rsparsev = SparseVector(float.(rand(MersenneTwister(123454321), Bool, 20)))
	@test densify(rsparsev) == Vector(rsparsev)
	@test densify(Diagonal(rsparsev)) == Diagonal(Vector(rsparsev))
end

@testset "allequal" begin
	@test allequal((true, true, true))
	@test allequal([true, true, true])
	@test !allequal((true, false, true))
	@test !allequal([true, false, true])
	@test !allequal(collect(1:4))
	@test allequal((false, false, false))
	@test allequal([false, false, false])
	@test allequal(ones(3))
	@test allequal(1, 1, 1)
end

@testset "threaded_replicate" begin
	rng = MersenneTwister(42);
	single_thread = replicate(10,use_threads=false) do; randn(rng, 1)[1] ; end
	rng = MersenneTwister(42);
	multi_thread = replicate(10,use_threads=true) do
		if Threads.threadid() % 2 == 0
			sleep(0.001)
		end
		r = randn(rng, 1)[1];
	end

	@test all(sort!(single_thread) .== sort!(multi_thread))
end

@testset "datasets" begin
	@test isa(MixedModels.datasets(), Vector{String})
	@test size(MixedModels.dataset(:dyestuff)) == (30, 2)
	@test size(MixedModels.dataset("dyestuff")) == (30, 2)
	@test_throws ArgumentError MixedModels.dataset(:foo)
end


text:
module Neighbors

using ..Dorothy.Geometry
using ..Dorothy.PBC

export
		NeighborGridBuffer, NeighborGrid, findnear, findnear!, findwithin,
		findwithin!, NeighborDistance, neighborlist, neighborlist!,
		NeighborTableBuffer, NeighborTable

mutable struct NeighborGridBuffer
	cells::Array{Vector{Int},3}
	I::Vector{Tuple{Int,Int,Int}}
end

NeighborGridBuffer() = NeighborGridBuffer(Array{Vector{Int},3}(undef, (0,0,0)),
		Tuple{Int,Int,Int}[])

function prepare!(buffer::NeighborGridBuffer, N::Tuple{Integer,Integer,Integer})
	nx, ny, nz = N
	mx, my, mz = size(buffer.cells)
	if nx > mx || ny > my || nz > mz
		buffer.cells = newcells = Array{Vector{Int},3}(undef, N .+ 1)
		for i in eachindex(newcells)
			newcells[i] = Int[]
		end
	else
		for oldcell in buffer.cells
			empty!(oldcell)
		end
	end
	empty!(buffer.I)
	buffer
end

abstract type NeighborGrid end

struct NonperiodicNeighborGrid <: NeighborGrid
	cells::Array{Vector{Int},3}
	I::Vector{Tuple{Int,Int,Int}}
	N::Tuple{Int,Int,Int}
	O::Vector3D
	D::Vector3D
end

struct PeriodicNeighborGrid <: NeighborGrid
	cells::Array{Vector{Int},3}
	I::Vector{Tuple{Int,Int,Int}}
	N::Tuple{Int,Int,Int}
	D::Vector3D
end

NeighborGrid(N::Tuple{Integer,Integer,Integer}, O::AbstractVector{<:Real},
		D::AbstractVector{<:Real}, buffer::NeighborGridBuffer) =
		NonperiodicNeighborGrid(buffer.cells, buffer.I, N, O, D)

NeighborGrid(N::Tuple{Integer,Integer,Integer}, D::AbstractVector{<:Real},
		buffer::NeighborGridBuffer) =
		PeriodicNeighborGrid(buffer.cells, buffer.I, N, D)

function NeighborGrid(Rw::AbstractVector{Vector3D},
		Kw::AbstractVector{Vector3D}, cell::Union{TriclinicPBC,Nothing},
		dmax::Real, buffer = NeighborGridBuffer())
	N, = params = neighborgridparams(Rw, Kw, cell, dmax)
	prepare!(buffer, N)
	nbrgrid = NeighborGrid(params..., buffer)
	append!(nbrgrid, Kw)
end

function neighborgridparams(Rw::AbstractVector{Vector3D},
		Kw::AbstractVector{Vector3D}, ::Nothing, dmax::Real)
	D = Vector3D(dmax, dmax, dmax)
	bb = extent(Rw)
	O = minimum(bb) - 2*D
	Dbb = dims(bb)
	nx = ceil(Int, Dbb.x / D.x) + 4
	ny = ceil(Int, Dbb.y / D.y) + 4
	nz = ceil(Int, Dbb.z / D.z) + 4
	(nx,ny,nz), O, D
end

function neighborgridparams(Rw::AbstractVector{Vector3D},
		Kw::AbstractVector{Vector3D}, cell::TriclinicPBC, dmax::Real)
	D = neighborgriddims(cell, dmax)
	nx = max(1, floor(Int, 1 / D.x))
	ny = max(1, floor(Int, 1 / D.y))
	nz = max(1, floor(Int, 1 / D.z))
	D = Vector3D(1.0/nx, 1.0/ny, 1.0/nz)
	(nx,ny,nz), D
end

function neighborgriddims(cell::TriclinicPBC, dmax::Real)
	rmax = norm(inv(cell) * Vector3D(dmax, dmax, dmax))
	Vector3D(rmax, rmax, rmax)
end

function neighborgriddims(cell::RhombododecahedralPBC, dmax::Real)
	rmax = √(2*dmax^2)
	gridcell = RhombododecahedralCell(rmax)
	inv(cell) * dims(gridcell)
end

neighborgriddims(cell::OrthorhombicPBC, dmax::Real) =
		inv(cell) * Vector3D(dmax, dmax, dmax)

Base.length(nbrgrid::NeighborGrid) = length(nbrgrid.I)

Base.size(nbrgrid::NeighborGrid) = nbrgrid.N

function Base.sizehint!(nbrgrid::NeighborGrid, n::Integer,
		ncell::Integer = cld(n, prod(size(nbrgrid))))
	@boundscheck n >= 0 && ncell >= 0 || error("expected positive length")
	sizehint!(nbrgrid.I, n)
	for cell in nbrgrid.cells
		sizehint!(cell, ncell)
	end
	nbrgrid
end

findcell(nbrgrid::NeighborGrid, i::Integer) = nbrgrid.I[i]

function findcell(nbrgrid::NonperiodicNeighborGrid, R::Vector3D)
	x = floor(Int, (R.x - nbrgrid.O.x) / nbrgrid.D.x) + 1
	y = floor(Int, (R.y - nbrgrid.O.y) / nbrgrid.D.y) + 1
	z = floor(Int, (R.z - nbrgrid.O.z) / nbrgrid.D.z) + 1
	x, y, z
end

function findcell(nbrgrid::PeriodicNeighborGrid, Rw::Vector3D)
	x = floor(Int, Rw.x / nbrgrid.D.x) + 1
	y = floor(Int, Rw.y / nbrgrid.D.y) + 1
	z = floor(Int, Rw.z / nbrgrid.D.z) + 1
	x, y, z
end

function Base.push!(nbrgrid::NeighborGrid,
		(x,y,z)::Tuple{Integer,Integer,Integer})
	nx, ny, nz = size(nbrgrid)
	!(x in 1:nx) || !(y in 1:ny) || !(z in 1:nz) &&
			error("position is outside the grid")
	push!(nbrgrid.cells[x,y,z], length(nbrgrid.I) + 1)
	push!(nbrgrid.I, (x,y,z))
	nbrgrid
end

Base.push!(nbrgrid::NeighborGrid, R::Vector3D) =
		push!(nbrgrid, findcell(nbrgrid, R))

function Base.append!(nbrgrid::NeighborGrid, R::AbstractVector{Vector3D})
	sizehint!(nbrgrid, length(nbrgrid) + length(R))
	for Ri in R
		push!(nbrgrid, Ri)
	end
	nbrgrid
end

findnear(R::AbstractVector{<:Vector3D}, i::Integer, dmax::Real) =
		findnear!(Int[], R, i, dmax)

findnear(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findnear!(Int[], R, i, cell, dmax)

findnear(Rw::AbstractVector{<:Vector3D}, Kw::AbstractVector{<:Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findnear!(Int[], Rw, Kw, i, cell, dmax)

findnear!(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		i::Integer, dmax::Real) = findnear!(dest, R, i, nothing, dmax)

findnear!(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findnear!(dest, pbcpos(R, cell)..., i, cell, dmax)

function findnear!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{<:Vector3D}, Kw::AbstractVector{<:Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real)
	nbrgrid = NeighborGrid(Rw, Kw, cell, dmax)
	findnear!(dest, nbrgrid, i)
end

findnear(R::AbstractVector{<:Vector3D}, Rref::AbstractVector{<:Real},
		dmax::Real) = findnear!(Int[], R, Rref, dmax)

findnear(R::AbstractVector{<:Vector3D}, Rref::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findnear!(Int[], R, Rref, cell, dmax)

findnear(Rw::AbstractVector{<:Vector3D}, Kw::AbstractVector{<:Vector3D},
		Rrefw::AbstractVector{<:Real}, Krefw::AbstractVector{<:Vector3D},
		cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findnear!(Int[], Rw, Kw, Rrefw, Krefw, cell, dmax)

findnear!(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		Rref::AbstractVector{<:Real}, dmax::Real) =
		findnear!(dest, R, Rref, nothing, dmax)

findnear!(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		Rref::AbstractVector{<:Real}, cell::Union{TriclinicPBC,Nothing},
		dmax::Real) =
		findnear!(dest, pbcpos(R, cell)..., pbcpos(Rref, cell)..., cell, dmax)

findnear!(dest::AbstractVector{<:Integer}, Rw::AbstractVector{<:Vector3D},
		Kw::AbstractVector{<:Vector3D}, Rrefw::AbstractVector{<:Real},
		Krefw::AbstractVector{<:Vector3D}, cell::Union{TriclinicPBC,Nothing},
		dmax::Real) =
		findnear!(dest, Rw, Kw, Vector3D(Rrefw), Vector3D(Krefw), cell, dmax)

function findnear!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{<:Vector3D}, Kw::AbstractVector{<:Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real)
	nbrgrid = NeighborGrid(Rw, Kw, cell, dmax)
	findnear!(dest, nbrgrid, Krefw)
end

findnear!(dest::AbstractVector{<:Integer}, nbrgrid::NeighborGrid, i::Integer) =
		findnear!(dest, nbrgrid, findcell(nbrgrid, i))

findnear!(dest::AbstractVector{<:Integer}, nbrgrid::NeighborGrid, R::Vector3D) =
		findnear!(dest, nbrgrid, findcell(nbrgrid, R))

function findnear!(dest::AbstractVector{<:Integer},
		nbrgrid::NonperiodicNeighborGrid,
		(x,y,z)::Tuple{Integer,Integer,Integer})
	empty!(dest)
	nx, ny, nz = size(nbrgrid)
	for xi = x-1:x+1
		for yi = y-1:y+1
			for zi = z-1:z+1
				if 1 <= xi <= nx && 1 <= yi <= ny && 1 <= zi <= nz
					for i in nbrgrid.cells[xi,yi,zi]
						push!(dest, i)
					end
				end
			end
		end
	end
	dest
end

function findnear!(dest::AbstractVector{<:Integer},
		nbrgrid::PeriodicNeighborGrid, (x,y,z)::Tuple{Integer,Integer,Integer})
	empty!(dest)
	nx, ny, nz = size(nbrgrid)
	for xi = x-1:x+1
		for yi = y-1:y+1
			for zi = z-1:z+1
				if xi < 1
					xi = nx
				elseif xi > nx
					xi = 1
				end
				if yi < 1
					yi = ny
				elseif yi > ny
					yi = 1
				end
				if zi < 1
					zi = nz
				elseif zi > nz
					zi = 1
				end
				for i in nbrgrid.cells[xi,yi,zi]
					push!(dest, i)
				end
			end
		end
	end
	dest
end

findwithin(R::AbstractVector{Vector3D}, i::Integer, dmax::Real) =
		findwithin!(Int[], R, nothing, dmax)

findwithin(R::AbstractVector{Vector3D}, i::Integer,
		cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findwithin!(Int[], pbcpos(R, cell)..., i, cell, dmax)

findwithin(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findwithin!(Int[], Rw, Kw, i, cell, dmax)

findwithin!(dest::AbstractVector{<:Integer}, R::AbstractVector{Vector3D},
		i::Integer, dmax::Real) = findwithin!(dest, R, i, nothing, dmax)

findwithin!(dest::AbstractVector{<:Integer}, R::AbstractVector{Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findwithin!(dest, pbcpos(R, cell)..., i, cell, dmax)

function findwithin!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D}, i::Integer,
		cell::Union{TriclinicPBC,Nothing}, dmax::Real)
	dmax2 = dmax^2
	Rwi, Kwi = Rw[i], Kw[i]
	empty!(dest)
	for j in findnear(Rw, Kw, i, cell, dmax)
		if sqmindist(Rw[j], Kw[j], Rwi, Kwi, cell) <= dmax2
			push!(dest, j)
		end
	end
	dest
end

findwithin(R::AbstractVector{Vector3D}, Rref::AbstractVector{<:Real},
		dmax::Real) = findwithin!(Int[], R, Rref, nothing, dmax)

findwithin(R::AbstractVector{Vector3D}, Rref::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findwithin!(Int[], R, Rref, cell, dmax)

findwithin(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::AbstractVector{<:Real}, Krefw::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		findwithin!(Int[], Rw, Kw, Rrefw, Krefw, cell, dmax)

findwithin!(dest::AbstractVector{<:Integer}, R::AbstractVector{Vector3D},
		Rref::AbstractVector{<:Real}, dmax::Real) =
		findwithin!(dest, R, Rref, nothing, dmax)

findwithin!(dest::AbstractVector{<:Integer}, R::AbstractVector{Vector3D},
		Rref::AbstractVector{<:Real}, cell::Union{TriclinicPBC,Nothing},
		dmax::Real) =
		findwithin!(dest, pbcpos(R, cell)..., pbcpos(Rref, cell)..., cell, dmax)

findwithin!(dest::AbstractVector{<:Integer}, Rw::AbstractVector{Vector3D},
		Kw::AbstractVector{Vector3D}, Rrefw::AbstractVector{<:Real},
		Krefw::AbstractVector{<:Real}, cell::Union{TriclinicPBC,Nothing},
		dmax::Real) =
		findwithin!(dest, Rw, Kw, Vector3D(Rrefw), Vector3D(Krefw), cell, dmax)

function findwithin!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real)
	dmax2 = dmax^2
	empty!(dest)
	for i in findnear(Rw, Kw, Rrefw, Krefw, cell, dmax)
		if sqmindist(Rw[i], Kw[i], Rrefw, Krefw, cell) <= dmax2
			push!(dest, i)
		end
	end
	dest
end

struct NeighborDistance
	i::Int
	d2::Real
end

function Base.iterate(nd::NeighborDistance, i::Integer = 1)
	if i > 2
		nothing
	else
		nd[i], i + 1
	end
end

Base.length(::NeighborDistance) = 2

function Base.getindex(nd::NeighborDistance, i::Integer)
	if i == 1
		nd.i
	elseif i == 2
		nd.d2
	else
		throw(BoundsError(nd, i))
	end
end

Base.firstindex(::NeighborDistance) = 1

Base.lastindex(nd::NeighborDistance) = length(nd)

neighborlist(R::AbstractVector{<:Vector3D}, i::Integer, dmax::Real) =
		neighborlist!(NeighborDistance[], R, i, dmax)

neighborlist(R::AbstractVector{<:Vector3D}, i::Integer,
		cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		neighborlist!(NeighborDistance[], R, i, cell, dmax)

neighborlist(Rw::AbstractVector{<:Vector3D}, Kw::AbstractVector{<:Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		neighborlist!(NeighborDistance[], Rw, Kw, i, cell, dmax)

neighborlist!(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		i::Integer, dmax::Real) = neighborlist!(dest, R, i, nothing, dmax)

neighborlist!(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		neighborlist!(dest, pbcpos(R, cell)..., i, cell, dmax)

function neighborlist!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{<:Vector3D}, Kw::AbstractVector{<:Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real)
	reachable = findnear(Rw, Kw, i, cell, dmax)
	neighborlist!(dest, Rw, Kw, i, cell, dmax, reachable)
end

neighborlist(R::AbstractVector{<:Vector3D}, Rref::AbstractVector{<:Real},
		dmax::Real) = neighborlist!(NeighborDistance[], R, Rref, dmax)

neighborlist(R::AbstractVector{<:Vector3D}, Rref::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		neighborlist!(NeighborDistance[], R, Rref, cell, dmax)

neighborlist(Rw::AbstractVector{<:Vector3D}, Kw::AbstractVector{<:Vector3D},
		Rrefw::AbstractVector{<:Real}, Krefw::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}, dmax::Real) =
		neighborlist!(NeighborDistance[], Rw, Kw, Rrefw, Krefw, cell, dmax)

neighborlist!(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		Rref::AbstractVector{<:Real}, dmax::Real) =
		neighborlist!(dest, R, Rref, nothing, dmax)

neighborlist!(dest::AbstractVector{<:Integer}, R::AbstractVector{<:Vector3D},
		Rref::AbstractVector{<:Real}, cell::Union{TriclinicPBC,Nothing},
		dmax::Real) = neighborlist!(dest, pbcpos(R, cell)...,
		pbcpos(Rref, cell)..., cell, dmax)

neighborlist!(dest::AbstractVector{<:Integer}, Rw::AbstractVector{<:Vector3D},
		Kw::AbstractVector{<:Vector3D}, Rrefw::AbstractVector{<:Real},
		Krefw::AbstractVector{<:Real}, cell::Union{TriclinicPBC,Nothing},
		dmax::Real) = neighborlist!(dest, Rw, Kw, Vector3D(Rrefw),
		Vector3D(Krefw), cell, dmax)

function neighborlist!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{<:Vector3D}, Kw::AbstractVector{<:Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real)
	reachable = findnear(Rw, Kw, Rrefw, Krefw, cell, dmax)
	neighborlist!(dest, Rw, Kw, Rrefw, Krefw, cell, dmax, reachable)
end

function neighborlist!(dest::AbstractVector{NeighborDistance},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real,
		reachable::AbstractVector{<:Integer})
	dmax2 = dmax^2
	Rwi, Kwi = Rw[i], Kw[i]
	empty!(dest)
	for j in reachable
		if j != i
			d2 = sqmindist(Rw[j], Kw[j], Rwi, Kwi, cell)
			if d2 <= dmax2
				push!(dest, NeighborDistance(j, d2))
			end
		end
	end
	sort!(dest, by = neighbor -> neighbor.d)
end

function neighborlist!(dest::AbstractVector{NeighborDistance},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real, reachable::AbstractVector{<:Integer})
	dmax2 = dmax^2
	empty!(dest)
	for i in reachable
		d2 = sqmindist(Rw[i], Kw[i], Rrefw, Krefw, cell)
		if d2 <= dmax2
			push!(dest, NeighborDistance(i, d2))
		end
	end
	sort!(dest, by = neighbor -> neighbor.d)
end

neighborlist(nbrlist::AbstractVector{NeighborDistance}, dmax::Real) =
		neighborlist!(NeighborDistance[], nbrlist, dmax)

function neighborlist!(dest::AbstractVector{NeighborDistance},
		nbrlist::AbstractVector{NeighborDistance}, dmax::Real)
	empty!(dest)
	dmax2 = dmax^2
	i = length(nbrlist)
	while i > 0
		if nbrlist[i].d2 > dmax2
			i -= 1
		else
			break
		end
	end
	sizehint!(dest, length(1:i))
	for j = 1:i
		push!(dest, nbrlist[j])
	end
	dest
end

struct NeighborTableBuffer
	grids::Dict{Int,NeighborGridBuffer}
	lists::Vector{Vector{NeighborDistance}}
	dmaxima::Vector{Float64}
	dest::Vector{Int}
end

NeighborTableBuffer() = NeighborTableBuffer(Dict{Int,NeighborGridBuffer}(),
		Vector{Vector{NeighborDistance}}(), Float64[], Int[])

function prepare!(buffer::NeighborTableBuffer, n::Integer)
	lists = buffer.lists
	dmaxima = buffer.dmaxima
	m = length(dmaxima)
	for list in lists
		empty!(list)
	end
	if n > m
		resize!(lists, n)
		for i = (m+1):n
			lists[i] = NeighborDistance[]
		end
	end
	resize!(dmaxima, n)
	fill!(dmaxima, -1.0)
	buffer
end

struct NeighborTable{T<:NeighborGrid}
	grids::Dict{Int,T}
	gridbuffers::Dict{Int,NeighborGridBuffer}
	lists::Vector{Vector{NeighborDistance}}
	dmaxima::Vector{Float64}
	dest::Vector{Int}

	NeighborTable{T}(grids::Dict{Int,T},
			gridbuffers::Dict{Int,NeighborGridBuffer},
			lists::Vector{Vector{NeighborDistance}}, dmaxima::Vector{Float64},
			dest::Vector{Int}) where {T<:NeighborGrid} =
			new(grids, gridbuffers, lists, dmaxima, dest)
end

function NeighborTable{T}(n::Integer, buffer::NeighborTableBuffer) where
		{T<:NeighborGrid}
	prepare!(buffer, n)
	NeighborTable{T}(Dict{Int,T}(), buffer.grids, buffer.lists, buffer.dmaxima,
			buffer.dest)
end

NeighborTable(n::Integer, cell::TriclinicPBC, buffer::NeighborTableBuffer) =
		NeighborTable{PeriodicNeighborGrid}(n, buffer)

NeighborTable(n::Integer, ::Nothing, buffer::NeighborTableBuffer) =
		NeighborTable{NonperiodicNeighborGrid}(n, buffer)

NeighborTable(n::Integer, cell::Union{TriclinicPBC,Nothing}) =
		NeighborTable(n, cell, NeighborTableBuffer())

Base.length(nbrtable::NeighborTable) = length(nbrtable.dmaxima)

function getnbrgrid(nbrtable::NeighborTable, Rw::AbstractVector{Vector3D},
		Kw::AbstractVector{Vector3D}, cell::Union{TriclinicPBC,Nothing},
		dmax::Real)
	dgrid = ceil(Int, dmax)
	get!(nbrtable.grids, dgrid) do
		buffer = get!(nbrtable.gridbuffers, dgrid) do
			NeighborGridBuffer()
		end
		NeighborGrid(Rw, Kw, cell, dmax, buffer)
	end
end

function getnbrlist(nbrtable::NeighborTable,  i::Integer,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		cell::Union{TriclinicPBC,Nothing}, dmax::Real)
	nbrlist = nbrtable.lists[i]
	if nbrtable.dmaxima[i] < dmax
		nbrgrid = getnbrgrid(nbrtable, Rw, Kw, cell, dmax)
		reachable = findnear!(nbrtable.dest, nbrgrid, i)
		neighborlist!(nbrlist, Rw, Kw, i, cell, dmax, reachable)
	end
	nbrlist
end

findnear(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D}, i::Integer,
		cell::Union{TriclinicPBC,Nothing}, dmax::Real,
		nbrtable::NeighborTable) =
		findnear!(Int[], Rw, Kw, i, cell, dmax, nbrtable)

function findnear!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D}, i::Integer,
		cell::Union{TriclinicPBC,Nothing}, dmax::Real,
		nbrtable::NeighborTable)
	if nbrtable.dmaxima[i] >= dmax
		dmax2 = dmax^2
		empty!(dest)
		push!(dest, i)
		for (j, d2) in nbrtable.lists[i]
			if d2 <= dmax2
				push!(dest, j)
			else
				break
			end
		end
		dest
	else
		nbrgrid = getnbrgrid(nbrtable, Rw, Kw, cell, dmax)
		findnear!(dest, nbrgrid, i)
	end
end

findnear(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real, nbrtable::NeighborTable) =
		findnear!(Int[], Rw, Kw, Rrefw, Krefw, cell, dmax, nbrtable)

function findnear!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real, nbrtable::NeighborTable)
	nbrgrid = getnbrgrid(nbrtable, Rw, Kw, cell, dmax)
	findnear!(dest, nbrgrid, Krefw)
end

findwithin(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real,
		nbrtable::NeighborTable) =
		findwithin!(Int[], Rw, Kw, i, cell, dmax, nbrtable)

function findwithin!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real,
		nbrtable::NeighborTable)
	dmax2 = dmax^2
	empty!(dest)
	push!(dest, i)
	for (j, d2) in getnbrlist(nbrtable, i, Rw, Kw, cell, dmax)
		if d2 <= dmax2
			push!(dest, i)
		else
			break
		end
	end
	dest
end

findwithin(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real, nbrtable::NeighborTable) =
		findwithin!(Int[], Rw, Kw, Rrefw, Krefw, cell, dmax, nbrtable)

function findwithin!(dest::AbstractVector{<:Integer},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real, nbrtable::NeighborTable)
	dmax2 = dmax^2
	empty!(dest)
	reachable =
			findnear!(nbrtable.dest, Rw, Kw, Rrefw, Krefw, cell, dmax, nbrtable)
	for i in reachable
		if sqmindist(Rw[i], Kw[i], Rrefw, Krefw, cell) <= dmax2
			push!(dest, i)
		end
	end
	dest
end

neighborlist(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real,
		nbrtable::NeighborTable) =
		neighborlist!(NeighborDistance[], Rw, Kw, i, cell, dmax, nbrtable)

function neighborlist!(dest::AbstractVector{NeighborDistance},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		i::Integer, cell::Union{TriclinicPBC,Nothing}, dmax::Real,
		nbrtable::NeighborTable)
	nbrlist = getnbrlist(nbrtable, i, Rw, Kw, cell, dmax)
	neighborlist!(dest, nbrlist, dmax)
end

neighborlist(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real, nbrtable::NeighborTable) = neighborlist!(NeighborDistance[],
		Rw, Kw, Rrefw, Krefw, cell, dmax, nbrtable)

function neighborlist!(dest::AbstractVector{NeighborDistance},
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rrefw::Vector3D, Krefw::Vector3D, cell::Union{TriclinicPBC,Nothing},
		dmax::Real, nbrtable::NeighborTable)
	reachable =
			findnear!(nbrtable.dest, Rw, Kw, Rrefw, Krefw, cell, dmax, nbrtable)
	neighborlist!(dest, Rw, Kw, Rrefw, Krefw, cell, dmax, reachable)
end

end # module


text:
#=

  Maximum flow problem in ConstraintSolver.jl

  From Winston "Operations Research", page 420f, 423f
  Sunco Oil example.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function max_flow_winston1(problem, print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    cap  = problem[:cap]
    arcs = problem[:arcs]
    n    = problem[:n]
    nodes = problem[:nodes]
    max_val = maximum(cap)
    num_arcs = length(cap)

    @variable(model, 0 <= x[1:n,1:n] <= max_val, Int)
    @variable(model, 0 <= z <= max_val*n, Int)

    @constraint(model, z == x[n,1])
    for i in 1:num_arcs
        @constraint(model,x[arcs[i,1], arcs[i,2]] <= cap[i])
    end
    for i in nodes
        @constraint(model,
          sum([x[arcs[k,1], arcs[k,2]] for k in 1:num_arcs if arcs[k,1] == i])
          ==
          sum([x[arcs[k,1], arcs[k,2]] for k in 1:num_arcs if arcs[k,2] == i])
        )
    end

    @objective(model,Max,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                println("x:")
                for row in eachrow(x_val)
                    println(row)
                end
                println("z:$z_val")
                println()
            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

problem = Dict(
    :cap => [2,3,3,4,2,1,100],
    :arcs =>  resize_matrix([[1, 2],
                            [1, 3],
                            [2, 3],
                            [2, 4],
                            [3, 5],
                            [4, 5],
                            [5, 1]]),
    :n => 5,
    :nodes => 1:5,
)

@time max_flow_winston1(problem,true,false)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank