text:
const scikit_common = Dict(
    # Gauge bosons
    "photon"=>22,
    "gamma"=>22,
    "Z_0"=>23,
    "W_plus"=>24,
    "W_minus"=>-24,
    # Leptons
    "e_minus"=>11,
    "e_plus"=>-11,
    "mu_minus"=>13,
    "mu_plus"=>-13,
    "tau_minus"=>15,
    "tau_plus"=>-15,
    "nu_e"=>12,
    "nu_e_bar"=>-12,
    "nu_mu"=>14,
    "nu_mu_bar"=>-14,
    "nu_tau"=>16,
    "nu_tau_bar"=>-16,
    # Quarkonia
    "Jpsi_1S"=>443,
    "psi_2S"=>100443,
    "eta_c_1S"=>441,
    "eta_c_2S"=>100441,
    "chi_c0_1P"=>10441,
    "chi_c1_1P"=>20443,
    "chi_c2_1P"=>445,
    "Upsilon_1S"=>553,
    "Upsilon_2S"=>100553,
    "Upsilon_3S"=>200553,
    "Upsilon_4S"=>300553,
    "eta_b_1S"=>551,
    "eta_b_2S"=>100551,
    "chi_b0_1P"=>10551,
    "chi_b1_1P"=>20553,
    "chi_b2_1P"=>555,
    # Light mesons
    "pi_0"=>111,
    "pi_plus"=>211,
    "pi_minus"=>-211,
    "a_0_980_0"=>9000111,
    "a_0_980_plus"=>9000211,
    "rho_770_0"=>113,
    "rho_770_plus"=>213,
    "rho_770_minus"=>-213,
    "a_1_1260_0"=>20113,
    "a_1_1260_plus"=>20213,
    "eta"=>221,
    "etap_958"=>331,
    "f_0_500"=>9000221,
    "f_0_980"=>9010221,
    "omega_782"=>223,
    "phi_1020"=>333,
    "K_L_0"=>130,
    "K_S_0"=>310,
    "K_plus"=>321,
    "K_minus"=>-321,
    "Kst_892_0"=>313,
    "Kst_892_0_bar"=>-313,
    "Kst_892_plus"=>323,
    "Kst_892_minus"=>-323,
    "K_0st_1430_0"=>10311,
    "K_0st_1430_0_bar"=>-10311,
    "K_0st_1430_plus"=>10321,
    "K_0st_1430_minus"=>-10321,
    # Charm mesons
    "D_0"=>421,
    "D_0_bar"=>-421,
    "D_plus"=>411,
    "D_minus"=>-411,
    "Dst_2007_0"=>423,
    "Dst_2007_0_bar"=>-423,
    "Dst_2010_plus"=>413,
    "Dst_2010_minus"=>-413,
    "D_s_plus"=>431,
    "D_s_minus"=>-431,
    "D_sst_plus"=>433,
    "D_sst_minus"=>-433,
    "D_s0st_2317_plus"=>10431,
    "D_s0st_2317_minus"=>-10431,
    "D_s1_2460_plus"=>20433,
    "D_s1_2460_minus"=>-20433,
    "D_s1_2536_plus"=>10433,
    "D_s1_2536_minus"=>-10433,
    # Beauty mesons
    "B_0"=>511,
    "B_0_bar"=>-511,
    "B_plus"=>521,
    "B_minus"=>-521,
    "B_s_0"=>531,
    "B_s_0_bar"=>-531,
    # Beauty, charmed mesons
    "B_c_plus"=>541,
    "B_c_minus"=>-541,
    # Light baryons
    "p"=>2212,
    "proton"=>2212,
    "p_bar"=>-2212,
    "antiproton"=>-2212,
    "n"=>2112,
    "neutron"=>2112,
    "n_bar"=>-2112,
    "antineutron"=>-2112,
    "Lambda"=>3122,
    "Lambda_bar"=>-3122,
    "Sigma_minus"=>3112,
    "Sigma_plus_bar"=>-3112,
    "Sigma_0"=>3212,
    "Sigma_0_bar"=>-3212,
    "Sigma_plus"=>3222,
    "Sigma_minus_bar"=>-3222,
    "Xi_minus"=>3312,
    "Xi_plus_bar"=>-3312,
    "Xi_0"=>3322,
    "Xi_0_bar"=>-3322,
    "Omega_minus"=>3334,
    "Omega_plus_bar"=>-3334,
    "Delta_1232_minus"=>1114,
    "Delta_1232_plus_bar"=>-1114,
    "Delta_1232_0"=>2114,
    "Delta_1232_0_bar"=>-2114,
    "Delta_1232_plus"=>2214,
    "Delta_1232_minus_bar"=>-2214,
    "Delta_1232_pp"=>2224,
    "Delta_1232_mm_bar"=>-2224,
    "Lambda_1520"=>3124,
    "Lambda_1520_bar"=>-3124,
    # Charm baryons
    "Lambda_c_plus"=>4122,
    "Lambda_c_minus_bar"=>-4122,
    # Beauty baryons
    "Lambda_b_0"=>5122,
    "Lambda_b_0_bar"=>-5122,
    "Xi_b_0"=>5232,
    "Xi_b_0_bar"=>-5232,
    "Xi_b_minus"=>5132,
    "Xi_b_plus_bar"=>-5132,
    "Omega_b_minus"=>5332,
    "Omega_b_plus_bar"=>-5332,
)
const more_common = Dict(
    "electron" => 11,
    "e-" => 11,
    "e+" => -11,
    "muon" => 13,
    "mu-" => 13,
    "mu+" => -13,
    "tauon" => 15,
    "tau-" => 15,
    "tau+" => -15,
    "H" => 25,
    "H0" => 25,
    "higgs" => 25,
)

"""
    Particle(name::String)

Construct particles from PDG name.

Or commonly alias, alias do not cover all particles.
"""
function Particle(n::String, anti=false)
    # determine if anti particle
    if startswith(n, "anti-")
        bn = n[6:end]
        anti = true
    elseif endswith(n, "~")
        bn = n[1:end-1]
        anti=true
    else
        bn = n
    end

    # find match to base name
    for S in (scikit_common, more_common, inv_catalog)
        if bn ∈ keys(S)
            p = Particle(S[bn])
            return anti ? -p : p
        end
    end
    error("$n not found.")
end

ext: jl


text:
module FilterTest
using SimplePNGs
using Test

include("common.jl")
using .TestCommon: load_json

pl(name) = SimplePNGs.load(joinpath("PngSuite", name*".png"))

@testset "Filtering test files" begin
    @testset "grayscale, no interlacing, filter-type 0" begin
        img1 = load_json("f00n0g08")
        img2 = pl("f00n0g08")
        @test img1 == img2
    end

    @testset "grayscale, no interlacing, filter-type 1" begin
        img1 = load_json("f01n0g08")
        img2 = pl("f01n0g08")
        @test img1 == img2
    end

    @testset "grayscale, no interlacing, filter-type 2" begin
        img1 = load_json("f02n0g08")
        img2 = pl("f02n0g08")
        @test img1 == img2
    end

    @testset "grayscale, no interlacing, filter-type 3" begin
        img1 = load_json("f03n0g08")
        img2 = pl("f03n0g08")
        @test img1 == img2
    end

    @testset "colour, no interlacing, filter-type 0" begin
        img1 = load_json("f00n2c08")
        img2 = pl("f00n2c08")
        @test img1 == img2
    end

    @testset "colour, no interlacing, filter-type 1" begin
        img1 = load_json("f01n2c08")
        img2 = pl("f01n2c08")
        @test img1 == img2
    end

    @testset "colour, no interlacing, filter-type 2" begin
        img1 = load_json("f02n2c08")
        img2 = pl("f02n2c08")
        @test img1 == img2
    end

    @testset "colour, no interlacing, filter-type 3" begin
        img1 = load_json("f03n2c08")
        img2 = pl("f03n2c08")
        @test img1 == img2
    end
end

end # module


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

# Run a local search algorithm to optimise a priority list given some objective function (objFn).
# The local search is simple hill climbing, and the neighbourhood of a priority list
# involves swapping or inserting items (station indices) within/into the list.
# The local search stops when a local optimum (as evaluated by simulating) is found.
# Ambulances are assumed to be identical.

using JEMSS
using Random
using Statistics
using StatsFuns
using Base.Iterators: Stateful

global ObjVal = Float

function priorityListLocalSearch!(sim::Simulation, priorityLists::Vector{PriorityList}; outputFolder::String = "")
	isdir(outputFolder) || mkdir(outputFolder)
	@assert(isdir(outputFolder))
	@assert(!isempty(priorityLists))
	@assert(all(priorityList -> checkPriorityList(priorityList, sim), priorityLists))

	# parameters
	global solFilename = "$outputFolder/solutions.csv" # final solutions from each local search
	global priorityListsOutputFilename = "$outputFolder/priority_lists.csv" # save final priority lists from each search
	global logFilename = "$outputFolder/log.csv" # log search progress to file
	global sense = :max # :min or :max; direction of optimisation for objective function
	global conf = 0.95 # statistical confidence level
	global doPrint = true
	# reps:
	global minNumReps = 2 # min number of reps to simulate each priority list
	global maxNumRepsList = [2,5] # = vcat(10 * 2 .^ [0:10;], Inf) # can start with a small maximum, and increase each time local optimum found
	global maxNumRepsIncrease = 1 # = 20 # max number of reps to increase by at a time when comparing two priority lists
	# parallel (multi-threading):
	global parallel = false
	global numThreads = 1

	# some parameter checks
	@assert(sense == :min || sense == :max)
	@assert(0 <= conf < 1) # should be between 0 and 1, but not equal to 1 otherwise confidence interval is infinite
	@assert(minNumReps >= 2) # cannot calculate confidence interval for less than 2 reps
	@assert(minNumReps <= sim.numReps)
	@assert(all(maxNumRepsList .>= minNumReps))
	@assert(issorted(maxNumRepsList, lt=<=)) # maxNumRepsList should be strictly increasing
	@assert(maxNumRepsIncrease >= 1)
	@assert(isa(parallel, Bool))
	@assert(numThreads >= 1)

	global nullObjVal = sense == :max ? -Inf : Inf

	# keep track of priority lists tried, and their objective values per replication
	global priorityListsObjVals = Dict{PriorityList, Vector{ObjVal}}() # priorityListsObjVals[priorityList][repIndex] gives sim replication objective value

	global priorityListsPeriodStatsList = Dict{PriorityList, Vector{SimPeriodStats}}() # priorityListsPeriodStatsList[priorityList][repIndex] gives sim replication period stats

	# save the locally optimal priority lists found for each search
	global priorityListSols = Vector{PriorityList}() # priorityListSols[i] is solution for priorityLists[i]

	global stationCapacities = [station.capacity for station in sim.stations] # stationCapacities[i] gives ambulance holding capacity of station i

	repeatedLocalSearch!(sim, priorityLists)
end

function objValLookup(priorityList::PriorityList, repIndex::Int)::ObjVal
	global priorityListsObjVals
	if haskey(priorityListsObjVals, priorityList) && length(priorityListsObjVals[priorityList]) >= repIndex
		return priorityListsObjVals[priorityList][repIndex]
	end
	return nullObjVal
end
function objValsLookup(priorityList::PriorityList, repIndices::Vector{Int})::Vector{ObjVal}
	return [objValLookup(priorityList, repIndex) for repIndex in repIndices]
end
function objValsLookup(priorityList::PriorityList)::Vector{ObjVal}
	global priorityListsObjVals
	return get(priorityListsObjVals, priorityList, ObjVal[])
end

function setPriorityListObjVal!(priorityList::PriorityList, repIndex::Int, objVal::ObjVal)
	global priorityListsObjVals
	objVals = get!(priorityListsObjVals, priorityList, ObjVal[])
	@assert(length(objVals) == repIndex - 1) # otherwise will overwrite existing value / will have missing values in objVals
	push!(objVals, objVal)
end

function periodStatsListLookup(priorityList::PriorityList)::Vector{SimPeriodStats}
	global priorityListsPeriodStatsList
	return get(priorityListsPeriodStatsList, priorityList, SimPeriodStats[])
end

function setPriorityListPeriodStats!(priorityList::PriorityList, repIndex::Int, period::SimPeriodStats)
	global priorityListsPeriodStatsList
	periods = get!(priorityListsPeriodStatsList, priorityList, SimPeriodStats[])
	@assert(length(periods) == repIndex - 1) # otherwise will overwrite existing value / will have missing values in periods
	push!(periods, period)
end
function setPriorityListPeriodStatsList!(priorityList::PriorityList, repIndices::UnitRange{Int}, periods::Vector{SimPeriodStats})
	for (repIndex, period) in zip(repIndices, periods)
		setPriorityListPeriodStats!(priorityList, repIndex, period)
	end
end

function resetPriorityListsPeriodStatsList!()
	global priorityListsPeriodStatsList
	empty!(priorityListsPeriodStatsList)
end

# Return mean and half-width of mean of x.
# Assumes that values in x are from population with normal distribution with unknown standard deviation.
function calcMeanAndHalfWidth(x::Vector{T}; conf::Float = conf)::Tuple{Float,Float} where T <: Real
	x = convert(Vector{Float}, x)
	return mean(x), tDistrHalfWidth(x; conf = conf)
end

# Objective value for a single period
function objFn(period::SimPeriodStats)::Float
	return period.call.numResponsesInTime / period.call.numCalls # fraction of calls reached in time
end

# Return objective values from multiple replications (repIndices).
# Uses look-up if possible.
# Mutates: sim.reps, priorityListsObjVals
function simObjVals!(sim::Simulation, priorityList::PriorityList, repIndices::UnitRange{Int})::Tuple{Vector{ObjVal}, Vector{SimPeriodStats}, Int}
	# look-up any saved data
	objValsFound = objValsLookup(priorityList)
	periodsFound = periodStatsListLookup(priorityList)
	n = min(length(objValsFound), length(periodsFound))
	lookupIndices= intersect(1:n, repIndices)
	simRepIndices = setdiff(repIndices, lookupIndices)

	objVals = ObjVal[]
	periods = SimPeriodStats[]
	if !isempty(simRepIndices)
		# need to simulate for simRepIndices
		global parallel, numThreads
		reps = sim.reps[simRepIndices]
		resetReps!(reps; parallel = parallel, numThreads = numThreads)
		for rep in reps applyPriorityList!(rep, priorityList) end
		simulateReps!(reps; parallel = parallel, numThreads = numThreads)
		periods = getRepsPeriodStatsList(reps)
		objVals = [objFn(period) for period in periods]
		for (i,repIndex) in enumerate(simRepIndices)
			if length(objValsFound) < repIndex setPriorityListObjVal!(priorityList, repIndex, objVals[i]) end
			if length(periodsFound) < repIndex setPriorityListPeriodStats!(priorityList, repIndex, periods[i]) end
		end
	end

	objVals = vcat(objValsFound[lookupIndices], objVals)
	periods = vcat(periodsFound[lookupIndices], periods)
	@assert(length(objVals) == length(repIndices))
	@assert(length(periods) == length(repIndices))

	return objVals, periods, length(lookupIndices)
end

function applyPriorityList!(sim::Simulation, priorityList::PriorityList)
	@assert(!sim.used)
	for i = 1:sim.numAmbs
		setAmbStation!(sim, sim.ambulances[i], sim.stations[priorityList[i]])
	end
	initPriorityList!(sim, priorityList)
end

function printPriorityList(priorityList::PriorityList)
	for i = 1:length(priorityList)
		println("item ", i, ": station ", priorityList[i])
	end
end

function getPeriodsStatsDict(periods::Vector{SimPeriodStats})::Dict{String,Any}
	global periodDuration, simStatsKeys, simStatsEmpty
	@assert(all(p -> isapprox(p.duration, periodDuration), periods))
	stats = flatten(statsDictFromPeriodStatsList(periods))
	return merge(Dict(["$(key)_mean" => stats[key].mean for key in simStatsKeys]), Dict(["$(key)_halfWidth" => stats[key].halfWidth for key in simStatsKeys]))
end

global simStatsEmpty = flatten(statsDictFromPeriodStatsList(SimPeriodStats[]))
global simStatsKeys = collect(keys(simStatsEmpty))
global simStatsEmpty = merge(Dict(["$(key)_mean" => NaN for key in simStatsKeys]), Dict(["$(key)_halfWidth" => NaN for key in simStatsKeys]))

global solFileHeader = vcat(["search", "numReps", "objVal", "objValHalfWidth"], sort(collect(keys(simStatsEmpty)))) # ... and priorityList
global logFileHeader = vcat(["search", "maxNumRepsListIndex", "maxNumReps", "iter", "move", "i", "j", "usedMove", "numReps", "numLookups", "searchDurationSeconds", "objVal", "objValHalfWidth", "bestObjVal", "objValDiff", "objValDiffHalfWidth"], sort(collect(keys(simStatsEmpty)))) # ... and priorityList
global logFileDict = Dict{String,Any}([s => "" for s in logFileHeader])

# mutates: file
function fileWriteDlmLine!(file::IOStream, fileHeader::Vector{String}, data::Dict{String,T}, priorityList::PriorityList) where T <: Any
	line = []
	for header in fileHeader
		push!(line, data[header])
	end
	line = vcat(line, priorityList)
	writeDlmLine!(file, line...)
	flush(file)
end

# perform local search, starting at each of the priority lists provided
# mutates: sim
function repeatedLocalSearch!(sim::Simulation, priorityLists::Vector{PriorityList})
	global doPrint

	# check that sim is using priority list for move-up
	@assert(sim.moveUpData.useMoveUp == true && sim.moveUpData.moveUpModule == priorityListModule)

	# open files for writing solution
	global solFilename, logFilename
	solFile = open(solFilename, "w")
	logFile = open(logFilename, "w")

	# write misc data to files
	global warmUpDuration, periodDuration, conf
	numSearches = length(priorityLists)
	miscTable = Table("miscData", ["numAmbs", "numStations", "numSearches", "conf"];
		rows = [[sim.numAmbs, sim.numStations, numSearches, conf]]) # note that numCalls may vary between sim replications
	repsTable = Table("repsData", ["numReps", "warmUpDuration", "periodDuration"];
		rows = [[sim.numReps, warmUpDuration, periodDuration]])
	writeTablesToFile!(logFile, [miscTable, repsTable])
	writeTablesToFile!(solFile, [miscTable, repsTable])

	# write file headers
	writeDlmLine!(solFile, "solutions")
	writeDlmLine!(logFile, "log")
	writeDlmLine!(solFile, solFileHeader..., ["item_$i stationIndex" for i = 1:sim.numAmbs]...)
	writeDlmLine!(logFile, logFileHeader..., ["item_$i stationIndex" for i = 1:sim.numAmbs]...)
	flush(solFile)
	flush(logFile)

	for i = 1:numSearches
		logFileDict["search"] = i

		# initial priority list
		priorityList = priorityLists[i]
		doPrint && println("Initial priority list $i (of $numSearches):")
		doPrint && printPriorityList(priorityList)

		# perform local search until no improvements can be made, then increase max number of reps allowed
		doPrint && println()
		global minNumReps, maxNumRepsList
		numReps = minNumReps # init
		resetPriorityListsPeriodStatsList!()
		for j = 1:length(maxNumRepsList)
			maxNumReps = maxNumRepsList[j] < sim.numReps ? Int(maxNumRepsList[j]) : sim.numReps
			logFileDict["maxNumRepsListIndex"] = j
			logFileDict["maxNumReps"] = maxNumReps

			doPrint && println("Search $i (of $numSearches)")
			doPrint && println("Setting maximum number of reps to ", maxNumReps)
			priorityList, numReps = localSearch!(sim, priorityList, numReps, maxNumReps, logFile)
			if maxNumReps == sim.numReps break end # increasing max number of reps will not help
		end
		doPrint && println("Finished local search $i (of $numSearches); solution:")
		doPrint && printPriorityList(priorityList)

		# save best priority list found for this search
		global priorityListSols
		push!(priorityListSols, priorityList)

		# write solution to file
		objVal, objValHalfWidth = calcMeanAndHalfWidth(objValsLookup(priorityList, [1:numReps;]))
		solFileDict = Dict("search" => i, "numReps" => numReps, "objVal" => objVal, "objValHalfWidth" => objValHalfWidth)
		merge!(solFileDict, getPeriodsStatsDict(periodStatsListLookup(priorityList)[1:numReps]))
		fileWriteDlmLine!(solFile, solFileHeader, solFileDict, priorityList)

		global priorityListsOutputFilename
		writePriorityListsFile(priorityListsOutputFilename, priorityListSols)
	end

	# print out all results from each finished local search
	doPrint && println()
	doPrint && println("Priority lists from completed local searches:")
	for i = 1:numSearches
		doPrint && println()
		doPrint && println("Search: ", i)
		priorityList = priorityListSols[i]
		doPrint && printPriorityList(priorityList)
		doPrint && println("objective value = ", mean(objValsLookup(priorityList)))
	end

	close(solFile)
	close(logFile)
end

# Local search of priority list, starting at priorityList.
# Move items in priority list, making changes that optimise the
# objective function (objFn) value for the given sense (:min or :max).
# Continue search until no improvement can be made.
# Mutates: sim, logFile
function localSearch!(sim::Simulation, priorityList::PriorityList, priorityListNumReps::Int, maxNumReps::Int, logFile::IOStream)::Tuple{PriorityList,Int}

	global minNumReps, maxNumRepsIncrease, logFileHeader, logFileDict, sense, doPrint

	# shorthand
	numAmbs = sim.numAmbs
	numStations = sim.numStations

	# track iterations
	iter = 1 # count iterations performed
	startTime = time()
	getSearchDuration() = round(time() - startTime, digits = 2)

	# track priority lists tried
	priorityListsTried = Set{PriorityList}()

	# keep track of how many reps were simulated for each priority list
	priorityListsRepsDone = Dict{PriorityList, Int}()
	getPriorityListRepsDone!(priorityList::PriorityList) = get!(priorityListsRepsDone, priorityList, 0)
	updatePriorityListsRepsDone!(priorityList::PriorityList, n::Int) =
		(priorityListsRepsDone[priorityList] = max(n, getPriorityListRepsDone!(priorityList)))

	# calculate objective value for starting point
	numReps = priorityListNumReps # do not just use minNumReps, as priorityList may have already been simulated a number of reps
	objVals, periods, numLookups = simObjVals!(sim, priorityList, 1:numReps)
	(objVal, objValHalfWidth) = calcMeanAndHalfWidth(objVals)
	updatePriorityListsRepsDone!(priorityList, numReps)
	push!(priorityListsTried, priorityList)
	doPrint && println("starting objective value: ", objVal)

	# starting solution is current best
	bestPriorityList = copy(priorityList)
	bestObjVals = objVals
	bestObjVal = objVal

	# write starting point
	for (key, value) in logFileDict
		if !in(key, ("search", "maxNumRepsListIndex", "maxNumReps"))
			logFileDict[key] = "" # reset value
		end
	end
	merge!(logFileDict, getPeriodsStatsDict(periods))
	merge!(logFileDict, Dict("iter" => iter, "numReps" => numReps, "numLookups" => numLookups, "searchDurationSeconds" => getSearchDuration(),
		"objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal, "objValDiff" => 0, "objValDiffHalfWidth" => 0))
	fileWriteDlmLine!(logFile, logFileHeader, logFileDict, priorityList)

	# local search
	doPrint && println("starting local search")
	endPoint = nothing
	iter += 1
	while true
		for i = (numAmbs+numStations):-1:1, j = 1:numAmbs, move = ["insert", "swap"]
			# starting with i > numAmbs will try insert/swap new items into priority list first (instead of only moving existing items within list); should help with speed

			# change priority list, using re-insertion / swapping
			# if i <= numAmbs, move item in priorityList[i] to priorityList[j]
			# if i > numAmbs, try move slot from station i-numAmbs (if num spare slots > 0)
			priorityList = copy(bestPriorityList)
			if move == "insert"
				insert!(priorityList, i, j)
			elseif move == "swap"
				swap!(priorityList, i, j)
			end

			newBestFound = false
			if !in(priorityList, priorityListsTried)
				doPrint && print(move, ": i = ", i, ", j = ", j)

				# simulate until there is a significant difference in objective value between priorityList and bestPriorityList,
				# or until running out of replications to simulate
				numReps = minNumReps
				meanDiff = halfWidth = 0.0 # init
				objVals = ObjVal[] # for priorityList
				periods = SimPeriodStats[] # for priorityList
				numLookups = 0
				while numReps <= maxNumReps
					newObjVals, newPeriods, newNumLookups = simObjVals!(sim, priorityList, (length(objVals)+1):numReps)
					push!(priorityListsTried, priorityList)
					push!(objVals, newObjVals...)
					push!(periods, newPeriods...)
					numLookups += newNumLookups
					updatePriorityListsRepsDone!(priorityList, numReps)
					if length(bestObjVals) < length(objVals)
						push!(bestObjVals, simObjVals!(sim, bestPriorityList, (length(bestObjVals)+1):numReps)[1]...)
						updatePriorityListsRepsDone!(bestPriorityList, numReps)
					end
					meanDiff, halfWidth = calcMeanAndHalfWidth(objVals[1:numReps] - bestObjVals[1:numReps])
					if numReps == maxNumReps
						break # ran out of reps, will have to guess which priority list is better based on mean
					elseif (sense == :max && meanDiff < -halfWidth) || (sense == :min && meanDiff > halfWidth)
						break # bestPriorityList is better, stop
					elseif (sense == :max && meanDiff > halfWidth) || (sense == :min && meanDiff < -halfWidth)
						# priorityList may be better
						# need to simulate priorityList up to much as bestPriorityList, to avoid possibility of infinite loop
						newNumReps = getPriorityListRepsDone!(bestPriorityList)
						if numReps < newNumReps
							numReps = min(maxNumReps, newNumReps)
						elseif numReps == newNumReps
							break # done, no need to simulate priorityList further
						else
							error("numReps > newNumReps, this should not happen.")
						end
					else # not sure which priority list is better, simulate more reps
						newNumReps = numReps * (halfWidth / meanDiff) ^ 2 # since standard error is roughly proportional to 1/sqrt(n) for n samples
						if isnan(newNumReps) newNumReps = maxNumReps end # meanDiff was zero, just use all reps
						numReps = ceil(Int, max(numReps + 1, min(newNumReps, numReps + maxNumRepsIncrease, maxNumReps))) # bound newNumReps
					end
				end

				objVal, objValHalfWidth = calcMeanAndHalfWidth(objVals)
				bestObjVal = mean(bestObjVals)
				# note that length(objVals) <= length(bestObjVals), as the priority lists may have been simulated different amounts (reps)
				doPrint && println("; objective value: ", objVal, " (best: ", bestObjVal, ")")

				usedMove = false
				if (sense == :max && meanDiff > 0) || (sense == :min && meanDiff < 0) # do not compare objVals and bestObjVals (or their means) as lengths may differ
					# keep change
					bestPriorityList = priorityList
					bestObjVals = objVals
					bestObjVal = objVal
					usedMove = true
					newBestFound = true

					resetPriorityListsPeriodStatsList!() # remove periods to save on memory; new solution has almost completely new neighbourhood
					setPriorityListPeriodStatsList!(priorityList, 1:length(periods), periods)

					doPrint && println("made ", move, ": i = ", i, ", j = ", j, "; new objective value: ", bestObjVal)
				end

				merge!(logFileDict, getPeriodsStatsDict(periods))
				merge!(logFileDict, Dict("iter" => iter, "move" => move, "i" => i, "j" => j, "numReps" => numReps, "numLookups" => numLookups, "usedMove" => Int(usedMove),
					"searchDurationSeconds" => getSearchDuration(), "objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal,
					"objValDiff" => meanDiff, "objValDiffHalfWidth" => halfWidth))
				fileWriteDlmLine!(logFile, logFileHeader, logFileDict, priorityList)

				iter += 1
			end

			if newBestFound || endPoint == nothing
				endPoint = (i,j,move)
			elseif endPoint == (i,j,move)
				doPrint && println("Found local optimum.")
				return bestPriorityList, getPriorityListRepsDone!(bestPriorityList)
			end
		end
	end
end

# for the given priority list, count how many slots are still available for the station
function countStationSpareSlots(priorityList::PriorityList, stationIndex::Int)
	global stationCapacities
	stationSpareSlots = stationCapacities[stationIndex]
	for i in priorityList
		if i == stationIndex
			stationSpareSlots -= 1
		end
	end
	return stationSpareSlots
end

# Insert item at index i into index j in priority list.
# If i > numAmbs, insertion is done with station index i-numAmbs.
# Return true if insertion completed, false otherwise.
function insert!(priorityList::PriorityList, i::Int, j::Int)
	if i == j return false end

	# shorthand
	numAmbs = length(priorityList)
	numStations = length(stationCapacities) # global stationCapacities

	if i <= numAmbs
		@assert(i >= 1)
		# re-insert item in priorityList[i] to priorityList[j]
		temp = priorityList[i]
		if i < j
			for k = i:j-1
				priorityList[k] = priorityList[k+1]
			end
		elseif i > j
			for k = i:-1:j+1
				priorityList[k] = priorityList[k-1]
			end
		else
			error()
		end
		priorityList[j] = temp
	else
		i -= numAmbs # so now i = station index
		@assert(1 <= i && i <= numStations)
		stationSpareSlots = countStationSpareSlots(priorityList, i)
		if stationSpareSlots == 0
			return false
		end
		# insert slot from station i to priorityList[j] (if stationSpareSlots > 0)
		temp = priorityList[numAmbs] # will be bumped from end of priorityList
		for k = numAmbs:-1:j+1
			priorityList[k] = priorityList[k-1]
		end
		priorityList[j] = i
	end

	return true
end

# Swap items at index i and index j in priority list.
# If i > numAmbs, swap is done with station index i-numAmbs.
# Return true if swap completed, false otherwise.
function swap!(priorityList::PriorityList, i::Int, j::Int)
	if i == j return false end

	# shorthand
	numAmbs = length(priorityList)
	numStations = length(stationCapacities) # global stationCapacities

	@assert(1 <= i && i <= numAmbs + numStations)
	@assert(1 <= j && j <= numAmbs)

	# perform swap if stations are different
	if i <= numAmbs
		if priorityList[i] == priorityList[j]
			return false
		end
		(priorityList[i], priorityList[j]) = (priorityList[j], priorityList[i])
	else # elseif i <= numAmbs + numStations
		i -= numAmbs # so now i = station index
		@assert(1 <= i && i <= numStations)
		stationSpareSlots = countStationSpareSlots(priorityList, i)
		if i == priorityList[j] || stationSpareSlots == 0
			return false
		end
		# swap i with priorityList[j]
		priorityList[j] = i
	end

	return true
end

# mutates: sim
function setSimStatsCapture!(sim::Simulation, periodDurationsIter::Stateful, warmUpDuration::Float)
	stats = sim.stats = deepcopy(sim.backup.stats)
	stats.doCapture = true
	(stats.periodDurationsIter, stats.warmUpDuration) = (periodDurationsIter, warmUpDuration)
	stats.nextCaptureTime = sim.startTime + (stats.warmUpDuration > 0 ? stats.warmUpDuration : first(stats.periodDurationsIter))
	sim.backup.stats = deepcopy(sim.stats)
end

# parameters
configFilename = "sim_config.xml"
outputFolder = "output"
numSearches = 1
priorityListRng = Random.MersenneTwister(0) # useful for reproducing results, if using random restarts
# replication parameters
warmUpDuration = 0.1
periodDuration = 1.0

# parameter checks
@assert(numSearches >= 1)
@assert(warmUpDuration >= 0)
@assert(periodDuration > 0)

# initialise sim
t = time()
@assert(isfile(configFilename))
println("Initialising simulation from config: ", configFilename)
sim = initSim(configFilename, createBackup = true)
@assert(sim.numReps >= 2, "Need at least two simulation replications (sim.reps) for statistics.")

# set sim stats capturing
setSimStatsCapture!(sim, Stateful(periodDuration), warmUpDuration)
periodEndTime = sim.startTime + warmUpDuration + periodDuration
for rep in sim.reps
	# statistics should stop collecting before last call arrives (if not earlier), to avoid cool-down period
	@assert(periodEndTime <= rep.calls[end].arrivalTime)
end

# generate priority lists
stationCapacities = [station.capacity for station in sim.stations]
priorityLists = [makeRandPriorityList(sim.numAmbs, sim.numStations; stationCapacities = stationCapacities, rng = priorityListRng) for i = 1:numSearches]

# set sim to use priority list for move up
setMoveUpModule!(sim, priorityListModule)

makeRepsRunnable!(sim) # copy changes to sim (stats, priority list) to sim.reps

# run
println("Starting local search(es).")
priorityListLocalSearch!(sim, priorityLists, outputFolder = outputFolder)
println("Total runtime: ", round(time()-t, digits = 2), " seconds")


text:
mutable struct ListNode{T} <: AbstractNode{T}
    prev::ListNode{T}
    next::ListNode{T}
    data::T
    ListNode{T}() where T =(x=new(); x.prev=x; x.next=x; x)
    ListNode{T}(p, n, d) where T =new(p, n, d)
end
ListNode(p, n, d::T) where T =ListNode{T}(p, n, d)

# Doubly linked list.
mutable struct LinkedList{T} <: AbstractList{T}
    node::ListNode{T}
    LinkedList{T}() where T =new(ListNode{T}())
end

mutable struct SListNode{T} <: AbstractNode{T}
    next::SListNode{T}
    data::T
    SListNode{T}() where T =(x=new(); x.next=x; x)
    SListNode{T}(n::SListNode{T}, d::T) where T =new(n, d)
end
SListNode(n::SListNode{T}, d::T) where T =SListNode{T}(n, d)

# Singly-linked list
mutable struct SLinkedList{T} <: AbstractList{T}
    # node is always the last element. Points to the first element.
    node::SListNode{T}
    SLinkedList{T}() where T =new(SListNode{T}())
end

function show(io::IO, l::AbstractList{T}) where T
    print(io, "$(typeof(l)){", string(T), "}(")
    middle=false
    for item in l
        if middle
            print(io, ", ")
        end
        show(io, item)
        middle=true
    end
    print(io, ")")
end

# The section titles work through sections of the manual.
#### Iteration

isdone(l::AbstractList, n::AbstractNode)=(n==l.node)
iterate(l::AbstractList, n::AbstractNode=l.node.next)= n==l.node ? nothing : (n.data, n.next)

struct ListIndexIterator{L}
    l::L
end

# Returns an iterator over indices.
# Use getindex, setindex! to find the item at this index.
keys(l::AbstractList)=ListIndexIterator(l)
length(l::ListIndexIterator) = length(l.l)
iterate(liter::ListIndexIterator, n::AbstractNode=liter.l.node.next)= n==liter.l.node ? nothing : (n, n.next)

#### General Collections
isempty(l::AbstractList)=(l.node.next==l.node)

function empty!(l::AbstractList)
    l.node.next=l.node
end

function length(l::AbstractList)
    cnt=0
    for n in l
        cnt+=1
    end
    cnt
end

# This is supposed to be an integer index, but
# we return the node as an index.
function lastindex(l::AbstractList)
    node=l.node.next
    while node.next!=l.node
        node=node.next
    end
    node
end
function lastindex(l::LinkedList)
    l.node.prev
end

#### Iterable Collections

function in(item, l::AbstractList)
    miss=false
    for e in l
        if e == item
            return true
        elseif e === missing
            miss=true
        end
    end
    miss ? missing : false
end

eltype(l::AbstractList{T}) where T = T

# Highest index in list for each value in a that is
# a member of the list.
function indexin(a, l::AbstractList)
    highest::Vector{Union{Nothing,T} where T<:AbstractNode}  = fill(nothing, length(a))
    for (node, data) in pairs(l)
        for (xidx, x) in enumerate(a)
            if isequal(x, data)
                highest[xidx]=node
                break
            end
        end
    end
    return highest
end

first(l::AbstractList)=l.node.next.data

"""
    last(l::AbstractList)

Traverse the list and return the value stored in the last node.
It is `O(n)` because of the traversal.
Fallback used for single-linked lists (`SLinkedList`).
"""
function last(l::AbstractList)
    lastd=l.node.data
    for d in l
        lastd=d
    end
    lastd
end

"""
    last(l::LinkedList)

Returns the value stored in the last node in `O(1)` time.
Does not make use of `lastindex`.
"""
function last(l::LinkedList)
    l.node.prev.data
end

#### Indexable Collections

# Treat the node as an index. It is also what
# is used for the state in iterators.
getindex(l::AbstractList, n::AbstractNode)=n.data
function setindex!(l::LinkedList, n::AbstractNode, d)
    n.data=d
end

#### Dequeues

# Breaking interface expectation to push multiple items
# so that we can return an index of the pushed item.
# Use append! for multiple items.
function push!(l::AbstractList, item)
    lnode=lastindex(l)
    lnode.next=SListNode(lnode.next, item)
    lnode.next
end
function push!(l::LinkedList, item)
    toadd=ListNode(l.node.prev, l.node, item)
    l.node.prev.next=toadd
    l.node.prev=toadd
    toadd
end

function pop!(l::AbstractList)
    node=l.node
    while node.next.next!=l.node
        node=node.next
    end
    d=node.next.data
    node.next=node.next.next
    d
end
function pop!(l::LinkedList)
    d=l.node.prev.data
    l.node.prev.prev.next=l.node
    l.node.prev=l.node.prev.prev
    d
end

# Breaking interface expectation because:
# Returns an index to the item instead of the collection.
# Takes only one value at a time. Use prepend! for multiple.
function pushfirst!(l::AbstractList, d)
    l.node.next=SListNode(l.node.next, d)
    l.node.next
end
function pushfirst!(l::LinkedList, d)
    toadd=ListNode(l.node, l.node.next, d)
    l.node.next.prev=toadd
    l.node.next=toadd
    toadd
end

function popfirst!(l::AbstractList)
    x=l.node.next.data
    l.node.next.next.prev=l.node
    l.node.next=l.node.next.next
    x
end
function popfirst!(l::SLinkedList)
    d=l.node.next.data
    l.node.next=l.node.next.next
    d
end

# Insert-after.
function insert!(l::AbstractList, n::AbstractNode, d)
    n.next=typeof(n)(n.next, d)
end
function insert!(l::LinkedList, n::ListNode, d)
    toadd=ListNode(n.prev, n, d)
    n.prev.next=toadd
    n.prev=toadd
    toadd
end

# Linear in the number of elements.
# Second argument is the state from an iterator.
function deleteat!(l::AbstractList, n::AbstractNode)
    prev=l.node
    while prev.next!=n
        prev=prev.next
    end
    prev.next=n.next
    l
end
function deleteat!(l::LinkedList, n::ListNode)
    n.prev.next=n.next
    n.next.prev=n.prev
    l
end

# Removal of a node, returning the value at that node.
function splice!(l::AbstractList, n::AbstractNode)
    prev=l.node
    while prev.next!=n
        prev=prev.next
    end
    prev.next=n.next
    n.data
end

function splice!(l::LinkedList, n::ListNode)
    n.prev.next=n.next
    n.next.prev=n.prev
    n.data
end

# Replacement of a node.
function splice!(l::AbstractList, n::AbstractNode, d)
    (d, n.data)=(n.data, d)
    d
end

function append!(l::AbstractList, items)
    lnode=lastindex(l)
    for i in items
        lnode.next=SListNode(lnode.next, i)
        lnode=lnode.next
    end
end
function append!(l::LinkedList, items)
    for i in items
        push!(l, i)
    end
end

function prepend!(l::AbstractList, items)
    for i in reverse(items)
        pushfirst!(l, i)
    end
    l
end
function prepend!(l::LinkedList, items)
    node=l.node # Invariant: Add after the node "node."
    for i in items
        toadd=ListNode(node, node.next, i)
        node.next.prev=toadd
        node.next=toadd
        node=toadd
    end
    l
end

#find the index (node) of the first element of l for which predicate returns true
function findfirst(predicate, l::AbstractList)
    for n in ListIndexIterator
        if predicate(n.data) return n end
    end
    return(nothing)
end

# returns the position of a node in a list
function indextoposition(n::AbstractNode, l::AbstractList)
    for (i, node) in enumerate(ListIndexIterator(l))
        if node === n return i end
    end
    return nothing
end
indextoposition(a::Vector, l::AbstractList) = map(x -> indextoposition(x,l) , a)
indextoposition(::Nothing, _) = nothing
function positiontoindex(i::Int, l::AbstractList)
    if i <= length(l)
        ii = 0
        for j in keys(l)
            ii += 1
            if ii === i
                return j
            end
        end
    else
        error("list is shorter than $i")
    end
end
positiontoindex(v::Vector, l::AbstractList) = map(x -> positiontoindex(x, l), v)
positiontoindex(::Nothing, _) = nothing

# getindex (positiontoindex will error if idx is invalid)
Base.getindex(lst::AbstractList, idx::Int) = lst[positiontoindex(idx, lst)]


text:
module Units

import Unitful
using Unitful: @unit
day = Unitful.d
week = Unitful.wk
year = Unitful.yr
@unit month "month" Month 2.628e6 * Unitful.s false

const days = day
const weeks = week
const months = month
const years = year
const January = 0month
const February = 1month
const March = 2months
const April = 3months
const May = 4months
const June = 5months
const July = 6months
const August = 7months
const September = 8months
const October = 9months
const November = 10months
const December = 11months

const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(Units)
end

export day, days, week, weeks, month, months, year, years, Rates,
January, February, March, April, May, June, July, August,
September, October, November, December

end


text:
#module API
using CEnum
using Libdl

# TODO: This is quite simplistic, better create a build of VEDA through Yggdrasil
const libveda = Libdl.find_library("libveda.so.0", ["/usr/local/ve/veda/lib64/"])
if isempty(libveda)
    @warn "Could not find libveda, VectorEngine.jl will not function"
end

# TODO: figure out why Clang.jl missed these
const __VEDAcontext = Cvoid
const __VEDAmodule = Cvoid
const veoargs = Cvoid
const VEDAhost_function = Ptr{Cvoid}

include(joinpath(@__DIR__, "..", "..", "gen", "libveda_common.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "libveda.jl"))
#end


text:
using Test
using Logging
using NetworkOptions
using NetworkOptions: CA_ROOTS_VARS, bundled_ca_roots, bundled_known_hosts

const pkg_dir = dirname(@__DIR__)

const TEST_URLS = [
    "" # not a valid host name
    "com"
    ".com" # not a valid host name
    "example.com"
    "user@example.com"
    "https://example.com"
    "xyz://example.com" # protocol doesn't matter
    "\x02H2ݼtOٲ\0RƆW9\e]B>#ǲR" # not a valid host name
]

const TRANSPORTS = [nothing, "ssl", "tls", "ssh", "xyz"]

host_variants(host::AbstractString) = Dict(
    "$host"              => true,
    ".$host"             => false,
    "$host."             => false,
    "/$host"             => false,
    "$host/"             => false,
    "user@$host"         => true,
    "user@$host:path"    => true,
    "user@$host:/path"   => true,
    "user@$host/"        => false,
    "user@$host/:"       => false,
    "user@$host/path"    => false,
    "user@$host/:path"   => false,
    "user@$host."        => false,
    "user@.$host"        => false,
    "@$host"             => false,
    "$host@"             => false,
    "://$host"           => false,
    "https://$host"      => true,
    "https://$host/"     => true,
    "https://$host/path" => true,
    "https://$host."     => false,
    "https://.$host"     => false,
    "https:/$host"       => false,
    "https:///$host"     => false,
    "xyz://$host"        => true,
    "xyz://$host."       => false,
    "xyz://.$host"       => false,
)

const VARIABLES = [
    "JULIA_SSL_CA_ROOTS_PATH"
    "SSL_CERT_DIR"
    "SSL_CERT_FILE"
    "JULIA_NO_VERIFY_HOSTS"
    "JULIA_SSH_NO_VERIFY_HOSTS"
    "JULIA_SSL_NO_VERIFY_HOSTS"
    "JULIA_ALWAYS_VERIFY_HOSTS"
    "SSH_DIR"
    "SSH_KEY_NAME"
    "SSH_KEY_PASS"
    "SSH_KEY_PATH"
    "SSH_PUB_KEY_PATH"
    "SSH_KNOWN_HOSTS_FILES"
]

const SAVED_VARS = Dict{String,Union{String,Nothing}}(
    var => nothing for var in VARIABLES
)

function save_env()
    for var in VARIABLES
        SAVED_VARS[var] = get(ENV, var, nothing)
    end
end

function reset_env()
    for var in VARIABLES
        val = get(SAVED_VARS, var, nothing)
        if val !== nothing
            ENV[var] = val
        else
            delete!(ENV, var)
        end
    end
end

function clear_env()
    for var in VARIABLES
        delete!(ENV, var)
    end
end

save_env()
clear_env()


text:
using Test
using LinearAlgebra: dot
using Random
using Distributions
using POMDPs
using FiniteHorizonPOMDPs
using POMDPModelTools
using POMDPModels
using POMDPPolicies
using FiVI
using PointBasedValueIteration

import Base.convert
import FiVI: backup, prob_o_given_b_a, b_o_a, init_belief_space, upper_bound_update, UB, corner_belief, z_k_t_a_o
import FiniteHorizonPOMDPs: InStageDistribution, FixedHorizonPOMDPWrapper

function convert_pbvi(::Type{Array{Float64, 1}}, b::FiniteHorizonPOMDPs.InStageDistribution{Array{Float64, 1}}, pomdp)
    belief = zeros(length(states(pomdp)))
    no_states = length(stage_states(pomdp, 1))
    belief[1 + (b.stage - 1) * no_states : b.stage * no_states] .= b.d
    return belief
end

#test elementary part of backup
@testset "FiVI" begin
    hor = 5
    pomdp = fixhorizon(MiniHallway(), hor)

    @test corner_belief(3, 1).b == [1, 0, 0]

    a = 1
    obs_wrong = (1, 1)
    obs_correct = (1, 4)
    st = 4
    α = AlphaVec([1, 2, 3, 10, 20, 30, 1, 2, 3, 10, 20, 30, 100], 1)

    # observation from wrong stage test
    @test z_k_t_a_o(pomdp, a, obs_wrong, st, α) == zeros(length(stage_states(pomdp, st)))

    # correct observation stage test
    @test z_k_t_a_o(pomdp, a, obs_correct, st, α) == [1., 0., 0., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0.]
end

# test whole backup of tigerPOMDP
@testset "FiVI.jl" begin
    hor = 5
    pomdp = fixhorizon(TigerPOMDP(), hor)
    b_init = convert(Array{Float64, 1}, initialstate(pomdp), pomdp)
    b1 = [1., 0.]; b2 = [0., 1.]
    r = LazyCachedSAR(pomdp)

    # test backups of terminal states
    @test backup(pomdp, b_init, 5, [], r) == AlphaVec([-0.5, -0.5], 0)
    @test backup(pomdp, b1, 5, [], r) == AlphaVec([10., -0.], 1)
    @test backup(pomdp, b2, 5, [], r) == AlphaVec([0., 10.], 2)

    Γ = [AlphaVec([-0.5, -0.5], 0), AlphaVec([10., -0.], 1), AlphaVec([0., 10.], 2)]

    # test backup of stage t <= horizon(t)
    @test backup(pomdp, b_init, 4, Γ, r) == AlphaVec([8., 8.], 0)

    o = (true, 5)
    p = prob_o_given_b_a(pomdp, 5, b_init, 0, o)

    @test p == 0.5

    a = 0
    new_b = b_o_a(pomdp, 5, b_init, a, o, p)

    @test isapprox(new_b, [0.15, 0.85])

    # create states with already backed up values (simpler version)
    Bs, _ = init_belief_space(pomdp, horizon(pomdp)+1)
    for b in Bs
        b.v = -Inf
        for a in ordered_actions(pomdp)
            v_temp = dot([r(s, a) for s in stage_states(pomdp, horizon(pomdp))], b.b)
            b.v = max(b.v, v_temp)
        end
    end

    ub = UB(new_b, Bs)
    @test isapprox(ub, 6.7)

    Bs2, _ = init_belief_space(pomdp, horizon(pomdp))

    v = upper_bound_update(pomdp, Bs2[3], Bs, horizon(pomdp), r)
    @test isapprox(v, 5.7)
end

@testset "Comparison with PBVI" begin
    pomdps = [fixhorizon(MiniHallway(), 5)]

    for pomdp in pomdps
        pomdp = fixhorizon(MiniHallway(), 5)
        solver = FiVISolver(1, 50.)
        policy, vu = solve(solver, pomdp)

        pbvi_solver = PBVISolver(15, 1., false)
        pbvi_policy = solve(pbvi_solver, pomdp)

        @testset "$(typeof(pomdp)) Value function comparison" begin
            B = []
            if typeof(pomdp.m) == MiniHallway
                reachable_beliefs = [([0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0], 2), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0], 3),
                    ([0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 4), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0], 5),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0], 2), ([0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 3),
                    ([0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 4), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0], 5),
                    ([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 4), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0], 4),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 5), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0], 4),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0], 5), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0], 5),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 4), ([0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 2),
                    ([0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0], 3), ([0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 3),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 3), ([0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 5),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0], 4), ([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 3),
                    ([0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 2), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0], 4),
                    ([0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.08333333333333333, 0.0], 1),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0], 3), ([0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0], 4),
                    ([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 5), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0], 5),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 4), ([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 2),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0], 4), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0], 2),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0], 3), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0], 2),
                    ([0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0], 5), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0], 4),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 5), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0], 5),
                    ([0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 5), ([0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.05555555555555555, 0.33333333333333337], 6),
                    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0], 3), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0], 4),
                    ([0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 5), ([0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 3),
                    ([0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 2), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0], 5),
                    ([0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 4), ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0], 5),
                    ([0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 5)]
            else
                reachable_beliefs = [(rand(length(stage_states(pomdp, 1))), rand(1:horizon(pomdp))) for _ in 1:100]
            end

            reachable_beliefs_dist = [FiniteHorizonPOMDPs.InStageDistribution(d[1], d[2]) for d in reachable_beliefs]

            fivi_vs = [value(policy, b.d, b.stage) for b in reachable_beliefs_dist]
            pbvi_vs = [value(pbvi_policy,
                        convert_pbvi(Array{Float64, 1},
                        b, pomdp)) for b in reachable_beliefs_dist]

            @test isapprox(fivi_vs, pbvi_vs, rtol=.1)
        end

        @testset "$(typeof(pomdp)) Simulation results comparison" begin
            no_simulations = typeof(pomdp.m) == MiniHallway ? 1 : 10_000
            for s in states(pomdp)
                # Commented parts do not work yet as the POMDPSimulators does not easily support FH
                # println(s)
                # @show value(policy, Deterministic(s))
                # @show value(sarsop_policy, Deterministic(s))
                #
                # @show action(policy, Deterministic(s))
                # @show action(sarsop_policy, Deterministic(s))
                #
                # @show mean([simulate(RolloutSimulator(max_steps = 100), pomdp, policy, updater(policy), Deterministic(s)) for i in 1:no_simulations])
                # @show mean([simulate(RolloutSimulator(max_steps = 100), pomdp, sarsop_policy, updater(sarsop_policy), Deterministic(s)) for i in 1:no_simulations])

                @test isapprox(value(policy, Deterministic(s[1]), s[2]), value(pbvi_policy, Deterministic(s)), rtol=0.1) || value(policy, Deterministic(s[1]), s[2]) >  value(pbvi_policy, Deterministic(s))
                # @test isapprox( mean([simulate(RolloutSimulator(max_steps = 100), pomdp, policy, updater(policy), Deterministic(s)) for i in 1:no_simulations]),
                #                 mean([simulate(RolloutSimulator(max_steps = 100), pomdp, pbvi_policy, updater(pbvi_policy), Deterministic(s)) for i in 1:no_simulations]),
                #                 rtol=0.1)
            end
        end
    end
end


text:
using Test
using MPI
using CLIMA.Mesh.Topologies

function main()
  MPI.Init()

  comm = MPI.COMM_WORLD
  crank = MPI.Comm_rank(comm)
  csize = MPI.Comm_size(comm)

  @assert csize == 2

  topology = StackedBrickTopology(comm, (1:4,5:8,9:12), periodicity=(false,true,false),
                                  boundary=((1,2),(3,4),(5,6)))

  elems       = topology.elems
  realelems   = topology.realelems
  ghostelems  = topology.ghostelems
  sendelems   = topology.sendelems
  elemtocoord = topology.elemtocoord
  elemtoelem  = topology.elemtoelem
  elemtoface  = topology.elemtoface
  elemtoordr  = topology.elemtoordr
  elemtobndy  = topology.elemtobndy
  nabrtorank  = topology.nabrtorank
  nabrtorecv  = topology.nabrtorecv
  nabrtosend  = topology.nabrtosend

  globalelemtoface = [1  1  1  2  2  2  2  2  2  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2
                      1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2
                      4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4
                      3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3
                      5  6  6  5  6  6  5  6  6  5  6  6  5  6  6  5  6  6  5  6  6  5  6  6  5  6  6
                      5  5  6  5  5  6  5  5  6  5  5  6  5  5  6  5  5  6  5  5  6  5  5  6  5  5  6]

  globalelemtoordr = ones(Int, size(globalelemtoface))

  globalelemtocoord = Array{Int}(undef, 3, 8, 27)

  globalelemtocoord[:,:, 1] = [1 2 1 2 1 2 1 2; 5 5 6 6 5 5 6 6;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:, 2] = [1 2 1 2 1 2 1 2; 5 5 6 6 5 5 6 6; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:, 3] = [1 2 1 2 1 2 1 2; 5 5 6 6 5 5 6 6; 11 11 11 11 12 12 12 12]
  globalelemtocoord[:,:, 4] = [2 3 2 3 2 3 2 3; 5 5 6 6 5 5 6 6;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:, 5] = [2 3 2 3 2 3 2 3; 5 5 6 6 5 5 6 6; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:, 6] = [2 3 2 3 2 3 2 3; 5 5 6 6 5 5 6 6; 11 11 11 11 12 12 12 12]
  globalelemtocoord[:,:, 7] = [2 3 2 3 2 3 2 3; 6 6 7 7 6 6 7 7;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:, 8] = [2 3 2 3 2 3 2 3; 6 6 7 7 6 6 7 7; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:, 9] = [2 3 2 3 2 3 2 3; 6 6 7 7 6 6 7 7; 11 11 11 11 12 12 12 12]
  globalelemtocoord[:,:,10] = [1 2 1 2 1 2 1 2; 6 6 7 7 6 6 7 7;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:,11] = [1 2 1 2 1 2 1 2; 6 6 7 7 6 6 7 7; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:,12] = [1 2 1 2 1 2 1 2; 6 6 7 7 6 6 7 7; 11 11 11 11 12 12 12 12]
  globalelemtocoord[:,:,13] = [1 2 1 2 1 2 1 2; 7 7 8 8 7 7 8 8;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:,14] = [1 2 1 2 1 2 1 2; 7 7 8 8 7 7 8 8; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:,15] = [1 2 1 2 1 2 1 2; 7 7 8 8 7 7 8 8; 11 11 11 11 12 12 12 12]
  globalelemtocoord[:,:,16] = [2 3 2 3 2 3 2 3; 7 7 8 8 7 7 8 8;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:,17] = [2 3 2 3 2 3 2 3; 7 7 8 8 7 7 8 8; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:,18] = [2 3 2 3 2 3 2 3; 7 7 8 8 7 7 8 8; 11 11 11 11 12 12 12 12]
  globalelemtocoord[:,:,19] = [3 4 3 4 3 4 3 4; 7 7 8 8 7 7 8 8;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:,20] = [3 4 3 4 3 4 3 4; 7 7 8 8 7 7 8 8; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:,21] = [3 4 3 4 3 4 3 4; 7 7 8 8 7 7 8 8; 11 11 11 11 12 12 12 12]
  globalelemtocoord[:,:,22] = [3 4 3 4 3 4 3 4; 6 6 7 7 6 6 7 7;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:,23] = [3 4 3 4 3 4 3 4; 6 6 7 7 6 6 7 7; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:,24] = [3 4 3 4 3 4 3 4; 6 6 7 7 6 6 7 7; 11 11 11 11 12 12 12 12]
  globalelemtocoord[:,:,25] = [3 4 3 4 3 4 3 4; 5 5 6 6 5 5 6 6;  9  9  9  9 10 10 10 10]
  globalelemtocoord[:,:,26] = [3 4 3 4 3 4 3 4; 5 5 6 6 5 5 6 6; 10 10 10 10 11 11 11 11]
  globalelemtocoord[:,:,27] = [3 4 3 4 3 4 3 4; 5 5 6 6 5 5 6 6; 11 11 11 11 12 12 12 12]

  globalelemtobndy = [1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
                      0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 2 2 2 2
                      0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                      0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                      5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0
                      0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6]

  if crank == 0
    nrealelem = 12
    globalelems = [ 1, 2, 3, # 1
                    4, 5, 6, # 2
                    7, 8, 9, # 3
                   10,11,12, # 4
                   13,14,15, # 5
                   16,17,18, # 6
                   22,23,24, # 8
                   25,26,27] # 9

    elemtoelem_expect = [  1   2   3   1   2   3  10  11  12  10  11  12 13 14 15 16 17 18 19 20 21 22 23 24
                           4   5   6  22  23  24  19  20  21   7   8   9 13 14 15 16 17 18 19 20 21 22 23 24
                          13  14  15  16  17  18   4   5   6   1   2   3 13 14 15 16 17 18 19 20 21 22 23 24
                          10  11  12   7   8   9  16  17  18  13  14  15 13 14 15 16 17 18 19 20 21 22 23 24
                           1   1   2   4   4   5   7   7   8  10  10  11 13 14 15 16 17 18 19 20 21 22 23 24
                           2   3   3   5   6   6   8   9   9  11  12  12 13 14 15 16 17 18 19 20 21 22 23 24]

  nabrtorank_expect = [1]
  nabrtorecv_expect = UnitRange{Int}[1:12]
  nabrtosend_expect = UnitRange{Int}[1:12]
  elseif crank == 1
    nrealelem = 15

    globalelems = [13,14,15, # 5
                   16,17,18, # 6
                   19,20,21, # 7
                   22,23,24, # 8
                   25,26,27, # 9
                    1, 2, 3, # 1
                    4, 5, 6, # 2
                    7, 8, 9, # 3
                   10,11,12] # 4

    elemtoelem_expect = [ 1   2   3   1   2   3   4   5   6  22  23  24  19  20  21 16 17 18 19 20 21 22 23 24 25 26 27
                          4   5   6   7   8   9   7   8   9  10  11  12  13  14  15 16 17 18 19 20 21 22 23 24 25 26 27
                         25  26  27  22  23  24  10  11  12  13  14  15   7   8   9 16 17 18 19 20 21 22 23 24 25 26 27
                         16  17  18  19  20  21  13  14  15   7   8   9  10  11  12 16 17 18 19 20 21 22 23 24 25 26 27
                          1   1   2   4   4   5   7   7   8  10  10  11  13  13  14 16 17 18 19 20 21 22 23 24 25 26 27
                          2   3   3   5   6   6   8   9   9  11  12  12  14  15  15 16 17 18 19 20 21 22 23 24 25 26 27]

  nabrtorank_expect = [0]
  nabrtorecv_expect = UnitRange{Int}[1:12]
  nabrtosend_expect = UnitRange{Int}[1:12]
  end

  @test elems == 1:length(globalelems)
  @test realelems == 1:nrealelem
  @test ghostelems == nrealelem+1:length(globalelems)

  @test elemtocoord == globalelemtocoord[:,:,globalelems]
  @test elemtoface[:,realelems] == globalelemtoface[:,globalelems[realelems]]
  @test elemtoelem == elemtoelem_expect
  @test elemtobndy == globalelemtobndy[:,globalelems]
  @test elemtoordr == ones(eltype(elemtoordr), size(elemtoordr))
  @test nabrtorank == nabrtorank_expect
  @test nabrtorecv == nabrtorecv_expect
  @test nabrtosend == nabrtosend_expect
end

main()

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text: error("This purposefully dies")