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

# write sim object and output files

function writeAmbsFile(filename::String, ambulances::Vector{Ambulance}; writeOutputFields::Bool = false)
	header = ["index", "stationIndex", "class", "attributes"]
	row1(a::Ambulance) = [a.index, a.stationIndex, Int(a.class), json(a.attributes)]

	if writeOutputFields
		counts = (:numCallsTreated, :numCallsTransported, :numDispatches, :numDispatchesFromStation, :numDispatchesWhileMobilising, :numDispatchesOnRoad, :numDispatchesOnFree, :numRedispatches, :numMoveUps, :numMoveUpsFromStation, :numMoveUpsOnRoad, :numMoveUpsOnFree, :numMoveUpsReturnToPrevStation)
		countHeaders = [string(c) for c in counts]

		statuses = (setdiff(instances(AmbStatus), (ambNullStatus,))..., instances(AmbStatusSet)...)
		travelStatuses = (ambStatusSets[ambTravelling]..., instances(AmbStatusSet)...)
		statusDurationHeaders = [string("statusDurations_", string(s)) for s in statuses]
		statusDistanceHeaders = [string("statusDistances_", string(s)) for s in travelStatuses]

		header = vcat(header, countHeaders, statusDurationHeaders, statusDistanceHeaders)
		row2(a::Ambulance) = vcat([getfield(a,c) for c in counts], [a.statusDurations[s] for s in statuses], [a.statusDistances[s] for s in travelStatuses])

		# skipped: statusTransitionCounts (array)
	end

	row(a::Ambulance) = writeOutputFields ? vcat(row1(a), row2(a)) : row1(a)
	table = Table("ambulances", header; rows = [row(a) for a in ambulances])
	writeTablesToFile(filename, table)
end

function writeArcsFile(filename::String, arcs::Vector{Arc}, travelTimes::Array{Float,2}, arcForm::String)
	@assert(arcForm == "directed" || arcForm == "undirected")
	numModes = size(travelTimes,1)
	miscTable = Table("miscData", ["arcForm", "numModes"]; rows = [[arcForm, numModes]])
	mainHeaders = ["index", "fromNode", "toNode", "distance", ["mode_$i" for i = 1:numModes]...]
	fieldNames = setdiff(collect(keys(arcs[1].fields)), mainHeaders) # assume fields are same for all arcs
	arcsTable = Table("arcs", [mainHeaders..., fieldNames...];
		rows = [vcat(a.index, a.fromNodeIndex, a.toNodeIndex, a.distance, travelTimes[:,a.index]..., [a.fields[f] for f in fieldNames]...) for a in arcs])
	writeTablesToFile(filename, [miscTable, arcsTable])
end

function writeCallsFile(filename::String, startTime::Float, calls::Vector{Call}; writeOutputFields::Bool = false)
	@assert(length(calls) >= 1)
	miscTable = Table("miscData", ["startTime"]; rows = [[startTime]])

	header = ["index", "priority", "x", "y", "arrivalTime", "dispatchDelay", "onSceneDuration", "transport", "hospitalIndex", "handoverDuration", "attributes"]
	row1(c::Call) = [c.index, Int(c.priority), c.location.x, c.location.y, c.arrivalTime, c.dispatchDelay, c.onSceneDuration, Int(c.transport), c.hospitalIndex, c.handoverDuration, json(c.attributes)]

	if writeOutputFields
		header = vcat(header, ["dispatchTime", "ambArrivalTime", "hospitalArrivalTime", "numBumps", "wasQueued", "ambDispatchLoc.x", "ambDispatchLoc.y", "ambStatusBeforeDispatch", "chosenHospitalIndex", "queuedDuration", "bumpedDuration", "waitingForAmbDuration", "responseDuration", "ambGoingToCallDuration", "transportDuration", "serviceDuration"])
		row2(c::Call) = [c.dispatchTime, c.ambArrivalTime, c.hospitalArrivalTime, c.numBumps, Int(c.wasQueued), c.ambDispatchLoc.x, c.ambDispatchLoc.y, string(c.ambStatusBeforeDispatch), c.chosenHospitalIndex, c.queuedDuration, c.bumpedDuration, c.waitingForAmbDuration, c.responseDuration, c.ambGoingToCallDuration, c.transportDuration, c.serviceDuration]
	end

	row(c::Call) = writeOutputFields ? vcat(row1(c), row2(c)) : row1(c)
	callsTable = Table("calls", header; rows = [row(c) for c in calls])
	writeTablesToFile(filename, [miscTable, callsTable])
end

function writeDemandFile(filename::String, demand::Demand)
	demandRastersTable = Table("demandRasters", ["rasterIndex", "rasterFilename"];
		rows = [[i, demand.rasterFilenames[i]] for i = 1:demand.numRasters])

	demandModesTable = Table("demandModes", ["modeIndex", "rasterIndex", "priority", "arrivalRate"];
		rows = [[m.index, m.rasterIndex, string(m.priority), m.arrivalRate] for m in demand.modes])
	@assert(all(i -> demand.modes[i].index == i, 1:demand.numModes))

	# demand sets table
	dml = demand.modeLookup # shorthand
	@assert(size(dml) == (demand.numSets, numPriorities)) # should have value for each combination of demand mode and priority
	demandSetsTable = Table("demandSets", ["setIndex", "modeIndices"];
		rows = [[i, hcat(dml[i,:]...)] for i = 1:demand.numSets])

	# demand sets timing table
	startTimes = demand.setsStartTimes # shorthand
	setIndices = demand.setsTimeOrder # shorthand
	demandSetsTimingTable = Table("demandSetsTiming", ["startTime", "setIndex"];
		rows = [[startTimes[i], setIndices[i]] for i = 1:length(startTimes)])

	writeTablesToFile(filename, [demandRastersTable, demandModesTable, demandSetsTable, demandSetsTimingTable])
end

function writeDemandCoverageFile(filename::String, demandCoverage::DemandCoverage)
	dc = demandCoverage # shorthand
	coverTimesTable = Table("coverTimes", ["demandPriority", "coverTime"];
		rows = [[string(priority), coverTime] for (priority, coverTime) in dc.coverTimes])

	demandRasterCellNumPointsTable = Table("demandRasterCellNumPoints", ["rows", "cols"];
		rows = [[dc.rasterCellNumRows, dc.rasterCellNumCols]])

	writeTablesToFile(filename, [coverTimesTable, demandRasterCellNumPointsTable])
end

function writeHospitalsFile(filename::String, hospitals::Vector{Hospital}; writeOutputFields::Bool = false)
	header = ["index", "x", "y", "attributes"]
	row1(h::Hospital) = [h.index, h.location.x, h.location.y, json(h.attributes)]

	if writeOutputFields
		header = vcat(header, ["numCalls"])
		row2(h::Hospital) = [h.numCalls]
	end

	row(h::Hospital) = writeOutputFields ? vcat(row1(h), row2(h)) : row1(h)
	table = Table("hospitals", header, rows = [row(h) for h in hospitals])
	writeTablesToFile(filename, table)
end

function writeMapFile(filename::String, map::Map)
	table = Table("map", ["xMin", "xMax", "yMin", "yMax", "xScale", "yScale"];
		rows = [[map.xMin, map.xMax, map.yMin, map.yMax, map.xScale, map.yScale]])
	writeTablesToFile(filename, table)
end

function writeMultiCompTableFile(filename::String, multiCompTable::MultiCompTable)
	mct = multiCompTable # shorthand
	checkMultiCompTable(mct)
	numAmbs = length(mct)
	numStations = length(mct[1][1])
	miscTable = Table("miscData", ["numAmbs", "numStations"], rows = [[numAmbs, numStations]])
	table = Table("multiCompTable", vcat("numAmbs", ["station_$j" for j = 1:numStations]...);
		rows = vcat([[[i, row...] for row in mct[i]] for i=1:length(mct)]...))
	writeTablesToFile(filename, [miscTable, table])
end

function writeNodesFile(filename::String, nodes::Vector{Node})
	mainHeaders = ["index", "x", "y", "offRoadAccess"]
	fieldNames = setdiff(collect(keys(nodes[1].fields)), mainHeaders) # assume fields are same for all nodes
	table = Table("nodes", [mainHeaders..., fieldNames...];
		rows = [[n.index, n.location.x, n.location.y, Int(n.offRoadAccess), [n.fields[f] for f in fieldNames]...] for n in nodes])
	writeTablesToFile(filename, table)
end

function writeRedispatchFile(filename::String, redispatch::Redispatch)
	miscTable = Table("miscData", ["allowRedispatch"], rows = [[Int(redispatch.allow)]])
	conditionsTable = Table("redispatchConditions", ["fromCallPriority", "toCallPriority", "allowRedispatch"],
		rows = [[[string(p1), string(p2), Int(redispatch.conditions[Int(p1),Int(p2)])] for p1 in priorities, p2 in priorities]...])
	writeTablesToFile(filename, [miscTable, conditionsTable])
end

function writeRNetTravelsFile(filename::String, rNetTravels::Vector{NetTravel})
	n = length(rNetTravels)
	@assert(all(i -> rNetTravels[i].isReduced, 1:n))
	@assert(all(i -> rNetTravels[i].modeIndex == i, 1:n))
	# save only some field values to file
	rNetTravelsSave = [NetTravel(true) for i = 1:n]
	for i = 1:n, fname in (:modeIndex, :arcTimes, :arcDists, :fadjList, :spFadjIndex, :spNodePairArcIndex, :spFadjArcList)
		setfield!(rNetTravelsSave[i], fname, getfield(rNetTravels[i], fname))
	end
	serializeToFile(filename, rNetTravelsSave)
end

function writePrioritiesFile(filename::String, targetResponseDurations::Vector{Float})
	table = Table("priorities", ["priority", "name", "targetResponseDuration"];
		rows = [[i, string(Priority(i)), targetResponseDurations[i]] for i = 1:length(targetResponseDurations)])
	writeTablesToFile(filename, table)
end

function writePriorityListFile(filename::String, priorityList::PriorityList)
	numAmbs = length(priorityList)
	table = Table("priorityList", ["item", "stationIndex"];
		cols = [collect(1:numAmbs), priorityList])
	writeTablesToFile(filename, table)
end

function writePriorityListsFile(filename::String, priorityLists::Vector{PriorityList})
	n = length(priorityLists)
	numAmbs = length(priorityLists[1])
	@assert(all(priorityList -> length(priorityList) == numAmbs, priorityLists))
	table = Table("priorityLists", ["item", ["priorityList_$i stationIndex" for i = 1:n]...];
		cols = [collect(1:numAmbs), priorityLists...])
	writeTablesToFile(filename, table)
end

# note that this writes the raster to file in xyz format (plain text)
function writeRasterFile(filename::String, raster::Raster)
	open(filename, "w") do file
		for iy = 1:raster.ny, ix = 1:raster.nx
			x, y, z = raster.x[ix], raster.y[iy], raster.z[ix,iy]
			write(file, "$x $y $z", newline)
		end
	end
end

function writeStationsFile(filename::String, stations::Vector{Station})
	table = Table("stations", ["index", "x", "y", "capacity", "attributes"];
		rows = [[s.index, s.location.x, s.location.y, s.capacity, json(s.attributes)] for s in stations])
	writeTablesToFile(filename, table)
end

function writeTravelFile(filename::String, travel::Travel)
	# travel modes table
	travelModesTable = Table("travelModes", ["travelModeIndex", "offRoadSpeed"];
		rows = [[tm.index, tm.offRoadSpeed] for tm in travel.modes])

	# travel sets table
	tml = travel.modeLookup # shorthand
	@assert(size(tml) == (length(travel.modes), numPriorities)) # should have value for each combination of travel mode and priority
	travelSetsTable = Table("travelSets", ["travelSetIndex", "priority", "travelModeIndex"];
		rows = [[[i, string(Priority(j)), tml[i,j]] for i = 1:size(tml,1), j = 1:size(tml,2)]...])

	# travel sets timing table
	startTimes = travel.setsStartTimes # shorthand
	setIndices = travel.setsTimeOrder # shorthand
	travelSetsTimingTable = Table("travelSetsTiming", ["startTime", "travelSetIndex"];
		rows = [[startTimes[i], setIndices[i]] for i = 1:length(startTimes)])

	writeTablesToFile(filename, [travelModesTable, travelSetsTable, travelSetsTimingTable])
end

function writeZhangIpParamsFile(filename::String, zhangIpData::ZhangIpData)
	# shorthand
	zid = zhangIpData
	@unpack marginalBenefits, stationCapacities = zid
	numStations = length(marginalBenefits)
	@assert(numStations == length(stationCapacities))

	miscTable = Table("miscParams", ["travelTimeCost", "onRoadMoveUpDiscountFactor", "regretTravelTimeThreshold", "expectedHospitalHandoverDuration"];
		rows = [[zid.travelTimeCost, zid.onRoadMoveUpDiscountFactor, zid.regretTravelTimeThreshold, zid.expectedHospitalHandoverDuration]])

	stationCapacitiesTable = Table("stationCapacities", ["stationIndex", "capacity"];
		rows = [[i, stationCapacities[i]] for i = 1:numStations])

	maxNumSlots = maximum(stationCapacities)
	mb = Array{Any,2}(undef, numStations, maxNumSlots)
	mb[:] .= ""
	for i = 1:numStations, j = 1:length(marginalBenefits[i])
		mb[i,j] = marginalBenefits[i][j]
	end
	marginalBenefitsTable = Table("stationMarginalBenefits", ["stationIndex", ["slot_$i" for i = 1:maxNumSlots]...];
		rows = [[i, mb[i,:]...] for i = 1:numStations])

	writeTablesToFile(filename, [miscTable, stationCapacitiesTable, marginalBenefitsTable])
end

# opens output files for writing during simulation
# note: should have field sim.resim.use set/fixed before calling this function
function openOutputFiles!(sim::Simulation)
	if !sim.writeOutput; return; end

	outputFilePath(name::String) = sim.outputFiles[name].path

	# create output path if it does not already exist
	if !isdir(sim.outputPath)
		mkdir(sim.outputPath)
	end

	# open output files for writing
	# currently, only need to open events file
	if !sim.resim.use # otherwise, existing events file is used for resimulation
		sim.outputFiles["events"].iostream = open(outputFilePath("events"), "w")
		sim.eventsFile.io = sim.outputFiles["events"].iostream # shorthand

		# save checksum of input files
		inputFiles = sort([name for (name, file) in sim.inputFiles])
		fileChecksumStrings = [string("'", sim.inputFiles[name].checksum, "'") for name in inputFiles]
		writeTablesToFile!(sim.eventsFile.io, Table("inputFiles", ["name", "checksum"]; cols = [inputFiles, fileChecksumStrings]))

		# save events with a key, to reduce file size
		eventForms = instances(EventForm)
		eventKeys = [Int(eventForm) for eventForm in eventForms]
		eventNames = [string(eventForm) for eventForm in eventForms]
		eventFilterValues = [Int(sim.eventsFile.eventFilter[eventForm]) for eventForm in eventForms]
		writeTablesToFile!(sim.eventsFile.io, Table("eventDict", ["key", "name", "filter"]; cols = [eventKeys, eventNames, eventFilterValues]))

		# write events table name and header
		writeDlmLine!(sim.eventsFile.io, "events")
		writeDlmLine!(sim.eventsFile.io, "index", "parentIndex", "time", "eventKey", "ambIndex", "callIndex", "stationIndex")
	end
end

function writeOutputFiles(sim::Simulation)
	writeMiscOutputFiles(sim)
	if !isempty(sim.reps)
		writeStatsFiles(sim, sim.reps)
	else
		writeStatsFiles(sim)
	end
end

function writeMiscOutputFiles(sim::Simulation)
	outputFileKeys = keys(sim.outputFiles)
	outputFilePath(name::String) = sim.outputFiles[name].path
	if in("ambulances", outputFileKeys) writeAmbsFile(outputFilePath("ambulances"), sim.ambulances; writeOutputFields = true) end
	if in("calls", outputFileKeys) writeCallsFile(outputFilePath("calls"), sim.startTime, sim.calls; writeOutputFields = true) end
	if in("hospitals", outputFileKeys) writeHospitalsFile(outputFilePath("hospitals"), sim.hospitals; writeOutputFields = true) end
end

function closeOutputFiles!(sim::Simulation)
	if !sim.writeOutput; return; end

	if !sim.resim.use
		writeDlmLine!(sim.eventsFile.io, "end")
		close(sim.eventsFile.io)
	end
end

function writeEventToFile!(sim::Simulation, event::Event)
	if !sim.writeOutput || sim.resim.use; return; end
	if sim.eventsFile.eventFilter[event.form]
		writeDlmLine!(sim.eventsFile.io, event.index, event.parentIndex, @sprintf("%0.5f", event.time), Int(event.form), event.ambIndex, event.callIndex, event.stationIndex)
		# flush(sim.eventsFile.io)
	end
end

# write deployments to file
function writeDeploymentsFile(filename::String, deployments::Vector{Deployment}, numStations::Int)
	numAmbs = length(deployments[1])
	@assert(numStations >= maximum([maximum(d) for d in deployments]))
	numDeployments = length(deployments)

	miscTable = Table("miscData", ["numStations", "numDeployments"]; rows = [[numStations, numDeployments]])
	deploymentsTable = Table("deployments",
		["ambIndex", ["deployment_$i stationIndex" for i = 1:numDeployments]...];
		cols = [collect(1:numAmbs), deployments...])
	writeTablesToFile(filename, [miscTable, deploymentsTable])
end

# save batch mean response durations to file
function writeBatchMeanResponseDurationsFile(filename::String, batchMeanResponseDurations::Array{Float,2};
	batchTime = nullTime, startTime = nullTime, endTime = nullTime, responseDurationUnits::String = "minutes")
	@assert(batchTime != nullTime && startTime != nullTime && endTime != nullTime)
	x = batchMeanResponseDurations # shorthand
	(numRows, numCols) = size(x) # numRows = numSims, numCols = numBatches
	miscTable = Table("misc_data",
		["numSims", "numBatches", "batchTime", "startTime", "endTime", "response_duration_units"];
		rows=[[numRows, numCols, batchTime, startTime, endTime, responseDurationUnits]])
	avgBatchMeansTable = Table("avg_batch_mean_response_durations",
		["sim_index", "avg_batch_mean_response_duration", "standard_error"];
		rows = [[i, mean(x[i,:]), sem(x[i,:])] for i = 1:numRows])
	batchMeansTable = Table("batch_mean_response_durations",
		["batch_index", ["sim_$i" for i = 1:numRows]...];
		rows = [[i, x[:,i]...] for i = 1:numCols])
	writeTablesToFile(filename, [miscTable, avgBatchMeansTable, batchMeansTable])
end

# for batches in a single simulation replication
function writeStatsFiles(sim::Simulation)
	outputFileKeys = keys(sim.outputFiles)
	outputFilePath(name::String) = sim.outputFiles[name].path
	stats = sim.stats # shorthand
	if in("ambulancesStats", outputFileKeys) writeAmbsStatsFile(outputFilePath("ambulancesStats"), stats) end
	if in("callsStats", outputFileKeys) writeCallsStatsFile(outputFilePath("callsStats"), stats) end
	if in("hospitalsStats", outputFileKeys) writeHospitalsStatsFile(outputFilePath("hospitalsStats"), stats) end
	if in("stationsStats", outputFileKeys) writeStationsStatsFile(outputFilePath("stationsStats"), stats) end
	if in("statsDict", outputFileKeys) writeStatsDictFile(outputFilePath("statsDict"), getPeriodStatsList(sim)) end
end

# for independent simulation replications
function writeStatsFiles(sim::Simulation, reps::Vector{Simulation}; periodIndex::Int = nullIndex)
	outputFileKeys = keys(sim.outputFiles)
	outputFilePath(name::String) = sim.outputFiles[name].path
	periods = periodIndex == nullIndex ? getRepsPeriodStatsList(reps) : getRepsPeriodStatsList(reps, periodIndex = periodIndex)
	@assert(!isempty(periods))
	if in("ambulancesStats", outputFileKeys) writeAmbsStatsFile(outputFilePath("ambulancesStats"), periods) end
	if in("callsStats", outputFileKeys) writeCallsStatsFile(outputFilePath("callsStats"), periods) end
	if in("hospitalsStats", outputFileKeys) writeHospitalsStatsFile(outputFilePath("hospitalsStats"), periods) end
	if in("stationsStats", outputFileKeys) writeStationsStatsFile(outputFilePath("stationsStats"), periods) end
	if in("statsDict", outputFileKeys) writeStatsDictFile(outputFilePath("statsDict"), periods) end
end

function simStatsTimestampsTable(stats::SimStats)::Table
	return Table("timestamps", ["simStartTime", "warmUpEndTime", "lastCallArrivalTime", "simEndTime"];
		rows = [[stats.simStartTime, stats.warmUpEndTime, stats.lastCallArrivalTime, stats.simEndTime]])
end

function simStatsPeriodTable(period::SimPeriodStats)::Table
	return Table("period", ["startTime", "endTime", "duration"];
		rows = [[period.startTime, period.endTime, period.duration]])
end
function simStatsPeriodTable(periods::Vector{SimPeriodStats})::Table
	@assert(length(periods) >= 1)
	@unpack startTime, endTime, duration = periods[1]
	@assert(all(p -> isapprox(p.startTime, startTime), periods))
	@assert(all(p -> isapprox(p.endTime, endTime), periods))
	@assert(all(p -> isapprox(p.duration, duration), periods))
	return simStatsPeriodTable(periods[1])
end

function simStatsPeriodsTable(periods::Vector{SimPeriodStats})::Table
	periodsTable = Table("periods", ["periodIndex", "startTime", "endTime", "duration"];
		rows = [[i, p.startTime, p.endTime, p.duration] for (i,p) in enumerate(periods)])
	return periodsTable
end

function makeAmbsStatsTables(periods::Vector{SimPeriodStats})::Vector{Table}
	numAmbs = length(periods[1].ambulances) # shorthand
	@assert(all(p -> length(p.ambulances) == numAmbs, periods))

	ambulanceTables = Table[]
	counts = (:numCallsTreated, :numCallsTransported, :numDispatches, :numDispatchesFromStation, :numDispatchesWhileMobilising, :numDispatchesOnRoad, :numDispatchesOnFree, :numRedispatches, :numMoveUps, :numMoveUpsFromStation, :numMoveUpsOnRoad, :numMoveUpsOnFree, :numMoveUpsReturnToPrevStation)
	statuses = (setdiff(instances(AmbStatus), (ambNullStatus,))..., instances(AmbStatusSet)...)
	travelStatuses = (sort(collect(ambStatusSets[ambTravelling]))..., instances(AmbStatusSet)...)
	countHeaders = [string(c) for c in counts]
	statusDurationHeaders = [string("duration_", string(s)) for s in statuses]
	statusDistanceHeaders = [string("distance_", string(s)) for s in travelStatuses]
	getAmb(period::SimPeriodStats, ambIndex::Int) = ambIndex == 0 ? period.ambulance : period.ambulances[ambIndex]
	header = vcat("periodIndex", countHeaders, statusDurationHeaders, statusDistanceHeaders)
	row(a::AmbulanceStats) = vcat([getfield(a,c) for c in counts], [a.statusDurations[s] for s in statuses], [a.statusDistances[s] for s in travelStatuses])
	# skipped: statusTransitionCounts
	for i = 0:numAmbs
		name = i == 0 ? "ambulance" : "ambulances[$i]"
		ambulanceTable = Table(name, header;
			rows = [vcat(j, row(getAmb(p,i))) for (j,p) in enumerate(periods)])
		push!(ambulanceTables, ambulanceTable)
	end

	return ambulanceTables
end

# for batches in a single simulation replication
function writeAmbsStatsFile(filename::String, stats::SimStats)
	# shorthand
	periods = stats.periods
	numAmbs = length(periods[1].ambulances)

	miscTable = Table("miscData", ["mode", "numAmbs"]; rows = [["batch", numAmbs]])
	timestampsTable = simStatsTimestampsTable(stats)
	periodsTable = simStatsPeriodsTable(periods)
	ambulanceTables = makeAmbsStatsTables(periods)
	writeTablesToFile(filename, [miscTable, timestampsTable, periodsTable, ambulanceTables...])
end

# for independent simulation replications
function writeAmbsStatsFile(filename::String, periods::Vector{SimPeriodStats})
	numAmbs = length(periods[1].ambulances) # shorthand
	@assert(all(p -> length(p.ambulances) == numAmbs, periods))

	miscTable = Table("miscData", ["mode", "numAmbs"]; rows = [["replication", numAmbs]])
	periodTable = simStatsPeriodTable(periods)
	ambulanceTables = makeAmbsStatsTables(periods)
	writeTablesToFile(filename, [miscTable, periodTable, ambulanceTables...])
end

function makeCallsStatsTables(periods::Vector{SimPeriodStats})::Vector{Table}
	tables = Table[]
	fnames = setdiff(fieldnames(CallStats), (:callIndex, :responseDurationHist))
	getCallStats(period::SimPeriodStats, priority::Priority) = priority == nullPriority ? period.call : period.callPriorities[priority]

	for priority in (nullPriority, priorities...)
		name = priority == nullPriority ? "call" : "callPriorities[$priority]"
		push!(tables, Table(name, vcat("periodIndex", collect(string.(fnames)));
			rows = [vcat(i, [getfield(getCallStats(p, priority), fname) for fname in fnames]) for (i,p) in enumerate(periods)]))
	end

	# add response duration histogram tables
	if all(p -> p.call.responseDurationHist != nullHist, periods)
		# periods may not have same responseDurationHist size if using replications
		(n,i) = findmax([length(p.call.responseDurationHist.weights) for p in periods])
		edges = periods[i].call.responseDurationHist.edges[1]
		@assert(edges[1] == 0) # otherwise binWidth calc is wrong
		binWidth = edges[2]
		push!(tables, Table("responseDurationHist_params", ["binWidth", "binEdges"]; rows = [[string(binWidth), string(edges)]]))

		for priority in (nullPriority, priorities...)
			weights(h::Histogram) = vcat(h.weights, zeros(Int, n-length(h.weights)))
			name = priority == nullPriority ? "call.responseDurationHist" : "callPriorities[$priority].responseDurationHist"
			push!(tables, Table(name, vcat("periodIndex", edges[1:n]);
				rows = [vcat(i, weights(getCallStats(p, priority).responseDurationHist)) for (i,p) in enumerate(periods)]))
		end
	end

	return tables
end

# for batches in a single simulation replication
function writeCallsStatsFile(filename::String, stats::SimStats)
	# shorthand
	periods = stats.periods
	numCalls = stats.captures[end].call.numCalls

	miscTable = Table("miscData", ["mode", "numCalls"]; rows = [["batch", numCalls]])
	timestampsTable = simStatsTimestampsTable(stats)
	periodsTable = simStatsPeriodsTable(periods)
	callsTables = makeCallsStatsTables(periods)
	writeTablesToFile(filename, [miscTable, timestampsTable, periodsTable, callsTables...])
end

# for independent simulation replications
function writeCallsStatsFile(filename::String, periods::Vector{SimPeriodStats})
	miscTable = Table("miscData", ["mode"]; rows = [["replication"]]) # note that number of calls may vary between replications
	periodTable = simStatsPeriodTable(periods)
	callsTables = makeCallsStatsTables(periods)
	writeTablesToFile(filename, [miscTable, periodTable, callsTables...])
end

function makeHospitalsStatsTables(periods::Vector{SimPeriodStats})::Vector{Table}
	numHospitals = length(periods[1].hospitals) # shorthand
	@assert(all(p -> length(p.hospitals) == numHospitals, periods))

	hospitalsTables = Table[]
	fnames = setdiff(fieldnames(HospitalStats), (:hospitalIndex,))
	getHospital(period::SimPeriodStats, hospitalIndex::Int) = hospitalIndex == 0 ? period.hospital : period.hospitals[hospitalIndex]
	for i = 0:numHospitals
		name = i == 0 ? "hospital" : "hospitals[$i]"
		hospitalTable = Table(name, vcat("periodIndex", collect(string.(fnames)));
			rows = [vcat(j, [getfield(getHospital(p,i), fname) for fname in fnames]) for (j,p) in enumerate(periods)])
		push!(hospitalsTables, hospitalTable)
	end

	return hospitalsTables
end

# for batches in a single simulation replication
function writeHospitalsStatsFile(filename::String, stats::SimStats)
	# shorthand
	periods = stats.periods
	numHospitals = length(periods[1].hospitals)

	miscTable = Table("miscData", ["mode", "numHospitals"]; rows = [["batch", numHospitals]])
	timestampsTable = simStatsTimestampsTable(stats)
	periodsTable = simStatsPeriodsTable(periods)
	hospitalsTables = makeHospitalsStatsTables(periods)
	writeTablesToFile(filename, [miscTable, timestampsTable, periodsTable, hospitalsTables...])
end

# for independent simulation replications
function writeHospitalsStatsFile(filename::String, periods::Vector{SimPeriodStats})
	numHospitals = length(periods[1].hospitals) # shorthand
	@assert(all(p -> length(p.hospitals) == numHospitals, periods))

	miscTable = Table("miscData", ["mode", "numHospitals"]; rows = [["replication", numHospitals]])
	periodTable = simStatsPeriodTable(periods)
	hospitalsTables = makeHospitalsStatsTables(periods)
	writeTablesToFile(filename, [miscTable, periodTable, hospitalsTables...])
end

function makeStationsStatsTables(periods::Vector{SimPeriodStats})::Vector{Table}
	numStations = length(periods[1].stations) # shorthand
	@assert(all(p -> length(p.stations) == numStations, periods))

	getStation(period::SimPeriodStats, stationIndex::Int) = stationIndex == 0 ? period.station : period.stations[stationIndex]
	stationsNumIdleAmbsTotalDurationTable = Table("stations_numIdleAmbsTotalDuration", vcat("periodIndex", "station", ["stations[$i]" for i = 1:numStations]);
		rows = [vcat(j, [string(getStation(p,i).numIdleAmbsTotalDuration) for i = 0:numStations]) for (j,p) in enumerate(periods)])
	return [stationsNumIdleAmbsTotalDurationTable]
end

# for batches in a single simulation replication
function writeStationsStatsFile(filename::String, stats::SimStats)
	# shorthand
	periods = stats.periods
	numStations = length(periods[1].stations)

	miscTable = Table("miscData", ["mode", "numStations"]; rows = [["batch", numStations]])
	timestampsTable = simStatsTimestampsTable(stats)
	periodsTable = simStatsPeriodsTable(periods)
	stationsTables = makeStationsStatsTables(periods)
	writeTablesToFile(filename, [miscTable, timestampsTable, periodsTable, stationsTables...])
end

# for independent simulation replications
function writeStationsStatsFile(filename::String, periods::Vector{SimPeriodStats})
	numStations = length(periods[1].stations) # shorthand
	@assert(all(p -> length(p.stations) == numStations, periods))

	miscTable = Table("miscData", ["mode", "numStations"]; rows = [["replication", numStations]])
	periodTable = simStatsPeriodTable(periods)
	stationsTables = makeStationsStatsTables(periods)
	writeTablesToFile(filename, [miscTable, periodTable, stationsTables...])
end

# stats dict values into table
function makeStatsDictTable(statsDict::Dict{String,Any})
	statsDictFlat = flatten(statsDict)
	@assert(all(v -> isa(v, MeanAndHalfWidth), values(statsDictFlat))) # check that statsDictFlat really is flat
	statsDictTable = Table("statsDict", ["stat", "mean", "halfWidth"];
		rows = [[k, statsDictFlat[k].mean, statsDictFlat[k].halfWidth] for k in sort(collect(keys(statsDictFlat)))])
	return statsDictTable
end

function writeStatsDictFile(filename::String, periods::Vector{SimPeriodStats}; conf::Float = 0.95)
	periodDuration = 0.0
	if !isempty(periods)
		periodDuration = periods[1].duration
		@assert(all(p -> isapprox(p.duration, periodDuration), periods))
	end
	statsDictFlat = statsDictFromPeriodStatsList(periods, conf = conf) |> flatten
	statsDictTable = makeStatsDictTable(statsDictFlat)
	miscTable = Table("misc", ["conf", "numPeriods", "periodDuration"], rows = [[conf, length(periods), periodDuration]])
	writeTablesToFile(filename, [miscTable, statsDictTable])
end


text:
function Rhyolite.command_buffers(fstate::FrameState, app::Application, render_set::RenderSet)
    command_buffer, _... = unwrap(allocate_command_buffers(device(fstate), _CommandBufferAllocateInfo(CommandPool(app.gr.command_pools), COMMAND_BUFFER_LEVEL_PRIMARY, 1)))

    @record command_buffer begin
        cmd_begin_render_pass(
            RenderPassBeginInfo(
                fstate.render_pass,
                fstate.current_frame[].framebuffer,
                Rect2D(Offset2D(0, 0), Extent2D(extent(main_window(app.wm))...)),
                [ClearValue(ClearColorValue((0.05f0, 0.01f0, 0.1f0, 0.1f0)))],
            ),
            SUBPASS_CONTENTS_INLINE,
        )

        render_set = app.gr.render_set
        prepare!(render_set)
        bind_state = BindState(nothing, nothing, [], nothing, nothing)
        for render_info in render_set.render_infos
            bind_state = bind(command_buffer, BindRequirements(render_info), bind_state)
            draw_args = render_info.draw_args
            if !isnothing(bind_state.index_buffer)
                cmd_draw_indexed(draw_args[1:3]..., 0, draw_args[4])
            else
                cmd_draw(draw_args...)
            end
        end

        cmd_end_render_pass()
    end

    [command_buffer]
end

# type piracy.
@nospecialize
Vulkan.PrimitiveTopology(::Type{<:IndexList{Line}}) = PRIMITIVE_TOPOLOGY_LINE_LIST
Vulkan.PrimitiveTopology(::Type{<:Strip{Line}}) = PRIMITIVE_TOPOLOGY_LINE_STRIP
Vulkan.PrimitiveTopology(::Type{<:IndexList{Triangle}}) = PRIMITIVE_TOPOLOGY_TRIANGLE_LIST
Vulkan.PrimitiveTopology(::Type{<:Strip{Triangle}}) = PRIMITIVE_TOPOLOGY_TRIANGLE_STRIP
Vulkan.PrimitiveTopology(::Type{<:Fan{Triangle}}) = PRIMITIVE_TOPOLOGY_TRIANGLE_FAN
@specialize

function Vulkan.GraphicsPipelineCreateInfo(
        device::Created{Device,DeviceCreateInfo},
        render_pass::RenderPass,
        extent,
        shaders::Vector{Shader},
        mesh::Type{MeshVertexEncoding{I,T}},
        pipeline_layout::PipelineLayout,
    ) where {I,T}

    require_feature(device, :sampler_anisotropy)

    shader_stages = PipelineShaderStageCreateInfo.([shaders.vertex, shaders.fragment], shaders.specialization_constants)
    vertex_input_state = PipelineVertexInputStateCreateInfo([VertexInputBindingDescription(T, 0)], vertex_input_attribute_descriptions(T, 0))
    input_assembly_state = PipelineInputAssemblyStateCreateInfo(PrimitiveTopology(I), false)
    viewport_state = PipelineViewportStateCreateInfo(viewports = [Viewport(0, 0, extent..., 0, 1)], scissors = [Rect2D(Offset2D(0, 0), Extent2D(extent...))])
    rasterizer = PipelineRasterizationStateCreateInfo(false, false, POLYGON_MODE_FILL, FRONT_FACE_CLOCKWISE, false, 0.0, 0.0, 0.0, 1.0, cull_mode = CULL_MODE_BACK_BIT)
    multisample_state = PipelineMultisampleStateCreateInfo(SAMPLE_COUNT_1_BIT, false, 1.0, false, false)
    color_blend_attachment = PipelineColorBlendAttachmentState(
        true,
        BLEND_FACTOR_SRC_ALPHA,
        BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
        BLEND_OP_ADD,
        BLEND_FACTOR_SRC_ALPHA,
        BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
        BLEND_OP_ADD;
        color_write_mask = COLOR_COMPONENT_R_BIT | COLOR_COMPONENT_G_BIT | COLOR_COMPONENT_B_BIT,
    )
    color_blend_state = PipelineColorBlendStateCreateInfo(false, LOGIC_OP_CLEAR, [color_blend_attachment], Float32.((0.0, 0.0, 0.0, 0.0)))
    GraphicsPipelineCreateInfo(
        shader_stages,
        rasterizer,
        pipeline_layout,
        render_pass,
        0,
        0;
        vertex_input_state,
        multisample_state,
        color_blend_state,
        input_assembly_state,
        viewport_state,
    )
end

function setup_swapchain(device::Device, surface::SurfaceKHR)
    surface_formats = unwrap(get_physical_device_surface_formats_khr(device.physical_device, surface))
    format = first(surface_formats)
    attachment = AttachmentDescription(
        format.format,
        SAMPLE_COUNT_1_BIT,
        ATTACHMENT_LOAD_OP_CLEAR,
        ATTACHMENT_STORE_OP_STORE,
        ATTACHMENT_LOAD_OP_DONT_CARE,
        ATTACHMENT_STORE_OP_DONT_CARE,
        IMAGE_LAYOUT_UNDEFINED,
        IMAGE_LAYOUT_PRESENT_SRC_KHR,
    )
    capabilities = unwrap(get_physical_device_surface_capabilities_khr(device.physical_device, surface))
    swapchain_ci = SwapchainCreateInfoKHR(
        surface,
        3,
        format.format,
        format.color_space,
        Extent2D(capabilities.current_extent.vks.width, capabilities.current_extent.vks.height),
        1,
        IMAGE_USAGE_COLOR_ATTACHMENT_BIT,
        SHARING_MODE_EXCLUSIVE,
        [],
        capabilities.current_transform,
        COMPOSITE_ALPHA_OPAQUE_BIT_KHR,
        PRESENT_MODE_IMMEDIATE_KHR,
        false,
    )
    swapchain = Created(unwrap(create_swapchain_khr(device, swapchain_ci)), swapchain_ci)

    render_pass = RenderPass(
        device,
        [attachment],
        [
            SubpassDescription(
                PIPELINE_BIND_POINT_GRAPHICS,
                [],
                [AttachmentReference(0, IMAGE_LAYOUT_ATTACHMENT_OPTIMAL_KHR)],
                [],
            ),
        ],
        [
            SubpassDependency(
                SUBPASS_EXTERNAL,
                0;
                src_stage_mask = PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT,
                dst_stage_mask = PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT,
                dst_access_mask = ACCESS_COLOR_ATTACHMENT_WRITE_BIT,
            ),
        ],
    )

    swapchain, render_pass
end

function Vulkan.GraphicsPipelineCreateInfo(app::Application, widget)
    spec = gr
    layout_info = PipelineLayoutCreateInfo()
    GraphicsPipelineCreateInfo(
        gr.device,
        get(gr.pipeline_layout_cache),
    )
end

ext: jl


text:
run(`sudo apt-get install default-jre`)
Pkg.add("Taro")

Pkg.add("ZipFile")
Pkg.add("Requests")
Pkg.add("DataFrames")
Pkg.add("Iterators")


text:
# str_load.jl
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.6.0
# 16.12.2016
# Last Edit: 07.07.2021

"""# str_load

Load data from Aerodye STR files generated by TDLWintel

---

### Examples

`time,D = str_load(\"K:\\\\Data\\\\140113_135354.str\")`\n
`time,D = str_load([\"K:\\\\Data\\\\140113_135354.str\",\"K:\\\\Data\\\\140114_000000.str\"])`\n
`time,D = str_load(\"K:\\\\Data\\\\\")` # Load all STR files in the given directory and subdirectories\n
`time,D = str_load(\"K:\\\\Data\\\\\",verbose=true,cols=[\"Traw\",\"X1\"])` # Load only Traw and X1 from all files in the given directory and display information as it is processed\n
`time,D = str_load(\"K:\\\\Data\\\\\",DateTime(2014,6,27,11,58,0))` # Load files starting at the given timestamp\n
`time,D = str_load(\"K:\\\\Data\\\\\",DateTime(2014,6,27,11,58,0),DateTime(2014,7,4,12,0,0))` # Load files between the given timestamps\n

---

`time, D = str_load(F)`\n
* **time**::Array{DateTime,1} = Parsed time values from the STR file
* **D**::DataFrame = Data from the STR file
* **F**::Array{String,1} = Array of STR files to load

`time, D = str_load(F;verbose=false,cols=[])`\n
* **timeStr**::Array{DateTime,1} =
* **F**::String = File name and location
* **verbose**::Bool (optional) = Display what is happening, FALSE is default
* **cols**::Array{String,1} or Array{Symbol,1} (optional) = List of columns to return, [] is default

`time, D = str_load(F,mindate;verbose=false,cols=[])`\n
`time, D = str_load(F,mindate,maxdate;verbose=false,cols=[])`\n
* **mindate**::DateTime = Load files including and after this date and time
* **maxdate**::DateTime = Load files up to this date and time
"""
function str_load(Dr::String,mindate::DateTime;verbose::Bool=false,cols::Array{String,1}=String[],recur::Int=9999)
	## Load STR files in the given directory (Dr) including and beyond the given date ##
	return str_load(Dr,mindate,DateTime(9999,1,1,0,0,0);verbose=verbose,cols=cols,recur=recur)
end # End of str_load(Dr::String,mindate::DateTime;verbose::Bool=false,cols::Array{String,1}=String[],recur::Int=9999)

function str_load(Dr::String,mindate::DateTime,maxdate::DateTime;verbose::Bool=false,cols::Array{String,1}=String[],recur::Int=9999)
	## Load STR files in the given directory (Dr) between the given dates ##

	# Check for Directory
	if isdir(Dr) == false
		error("First input should be a directory")
	end

	# List Files
	verbose ? println("Listing Files") : nothing
	(Fstr,folders) = dirlist(Dr,regex=r"\d{6}_\d{6}\.str$",recur=recur) # List STR files

	# Parse Times
	(Tstr,Fstr) = aerodyne_parsetime(Fstr)

	# Sort Info
	f1 = findlast(Tstr .<= mindate)
	if isnothing(f1)
		f1 = 1
	end
	f2 = findlast(Tstr .< maxdate)
	f = collect(f1:f2)
	Tstr = Tstr[f]
	Fstr = Fstr[f]

	# Remove Files Out of Range
	begin
		# STR Files
		f = findall(mindate .<= Tstr .< maxdate)
		Tstr = Tstr[f]
		Fstr = Fstr[f]
	end

	# Load Data
	if verbose
		println("Loading:")
	end
	(Tstr,Dstr) = str_load(Fstr,verbose=verbose,cols=cols)

	# Remove Time Values Out of Range
	begin
		f = findall(mindate .<= Tstr .< maxdate)
		Tstr = Tstr[f]
		Dstr = Dstr[f,:]
	end

	return Tstr,Dstr
end # End of str_load(Dr::String,mindate::DateTime,maxdate::DateTime;verbose::Bool=false,cols::Array{String,1}=String[],recur::Int=9999)

function str_load(F::Array{T,1};verbose::Bool=false,cols::Array{String,1}=String[]) where T <: String
	## Load a list of file names ##

	# Sort files by dates
	(Ftime,F) = aerodyne_parsetime(F)

	# Load and concatinate data
	if !isempty(F)
		(t,D) = str_load(F[1],verbose=verbose,cols=cols) # Initial load
		for i=2:1:length(F)
			(tempT,tempD) = str_load(F[i],verbose=verbose,cols=cols)
			t = [t;tempT]
			D = [D;tempD]
		end
	else
		# No Files Found
		t = Vector{DateTime}(undef,0)
		D = DataFrame()
	end

	return t,D
end # End str_load(F::Array{T,1};verbose::Bool=false,cols::Array{String,1}=String[]) where T <: String

function str_load(F::String;verbose::Bool=false,cols::Array{String,1}=String[])
	ext = splitext(F)[2]

	verbose ? println("  " * F) : nothing

	## Check for a proper file
	if isdir(F) == true
		return str_load(F,DateTime(0,1,1,0,0,0),DateTime(9999,1,1,0,0,0),verbose=verbose,cols=cols)
	elseif isempty(ext) == true
		error("No file extension")
	elseif ext != ".str"
		error("Extension is not STR. Returning nothing...")
		return
	end

	#################
	##  Load Data  ##
	#################
	## Load Header Information
	fid = open(F,"r")
	h1 = readline(fid)
	h1 = split(h1,"SPEC:")[2] # Remove everything including and before the SPEC:
	h2 = readdlm(IOBuffer("\"" * replace(h1,"," => "\",\"") * "\""),',')
	h = Array{String}(undef,length(h2) + 2)
	h[1] = "time"
	h[end] = "Empty"
	for i=1:1:length(h2)
		h[i+1] = strip(String(h2[i])) # Remove leading and trailing whitespace
		h[i+1] = replace(h[i+1]," " => "_") # Replace remaining whitespace with _
	end
	close(fid)

	## Check for duplicate column names and adjust them
	unames = unique(h) # unique column names
	for i in unames
		check = 0
		for j = 1:length(h)
			if i == h[j]
				check += 1

				if check > 1
					h[j] = h[j] * "-" * string(check)
				end
			end
		end
	end

	## List column types
	coltypes = Any[Float64 for i=1:length(h)]
	for i=1:1:length(h)
		if h[i] == "SPEFile"
			coltypes[i] = Union{Missing,String}
		end

		if h[i] == "StatusW"
			coltypes[i] = Int
		end
	end

	## Todo: Replace column names with reasonable names

	## Load data
	D = DataFrame()
	try
		D = CSV.read(F,DataFrame;delim=" ",header=h,datarow=2)[:,1:end-1]
		#D = DataFrame!(CSV.File(F;delim=" ",header=h,datarow=2))[:,1:end-1]
	catch
		println("Cannot load " * F)
		error("ERROR loading files")
	end

	#########################
	##  Parse Time Format  ##
	#########################
	time = Array{DateTime}(undef,length(D.time)) # Preallocate time column
	secs = Dates.Second # Initialize so it doesn't have to do it every time in the loop
	millisecs = Dates.Millisecond # Initialize so it doesn't have to do it every time in the loop
	for i=1:1:length(D.time)
		secs = Dates.Second(floor(D.time[i]))
		millisecs = Dates.Millisecond(floor(1000(D.time[i] - floor(D.time[i]))))
		time[i] = DateTime(1904,1,1,0,0,0) + secs + millisecs
	end

	##################################
	##  Keep only selected columns  ##
	##################################
	if !isempty(cols)
		# Check cols' type
		if typeof(cols) != Array{Symbol,1} && typeof(cols) != Symbol
			temp = Array{Symbol}(undef,length(cols))
			for i=1:1:length(cols)
				temp[i] = Symbol(cols[i]) # Convert all the values to symbols
			end
			cols = temp
		end

		# Make sure each entry exists
		fields = Symbol.(names(D))
		cols_bool = fill!(Array{Bool}(undef,length(cols)),false) # Preallocate false
		for i=1:1:length(cols)
			for j=1:1:length(fields)
				if fields[j] == cols[i]
					cols_bool[i] = true
				end
			end
		end
		cols = cols[cols_bool]

		# Remove Unwanted column
		if isempty(cols)
			D = DataFrame()
			time = Array{DateTime}(0)
		else
			D = D[:,cols]
		end
	end

	return time,D
end # End of str_load(F::String;verbose::Bool=false,cols::Array{String,1}=String[])


text:
#########################################
####        JLMS and BC index        ####
#########################################

#? --------------- Truncated Normal --------------

function jlmsbc(::Type{Trun}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
                Y::Matrix, X::Matrix, Z::Matrix, q, W::Matrix, V::Matrix, # dum1,
                dum2)

    x_pre = X*coef[pos.begx : pos.endx]
    z_pre = Z*coef[pos.begz : pos.endz]  # mu
  # q_pre = Q*coef[pos.begq : pos.endq]
    w_pre = W*coef[pos.begw : pos.endw] # log_σᵤ²
    v_pre = V*coef[pos.begv : pos.endv] # log_σᵥ²

    μ   = z_pre
    σᵤ² = @. exp(w_pre)
    σᵤ  = @. exp(0.5 * w_pre)
    σᵥ² = @. exp(v_pre)
    σ²  = σᵤ² + σᵥ²
    ϵ   = PorC*(Y - x_pre)
    μₛ  = @. (σᵥ² * μ - σᵤ² * ϵ) / σ²
    σₛ  = @. sqrt((σᵥ² * σᵤ²) / σ²)

   jlms = @. (σₛ * normpdf(μₛ / σₛ )) / normcdf(μₛ / σₛ) + μₛ
     bc = @. exp( -μₛ + 0.5 * (σₛ)^2  ) * ( normcdf( (μₛ / σₛ) - σₛ ) / normcdf(μₛ / σₛ)  )

   return jlms, bc
end

#? ------------------ Half normal model -------------

function jlmsbc(::Type{Half}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
                Y::Matrix, X::Matrix, z, q, W::Matrix, V::Matrix, # dum1,
                dum2)

    x_pre = X*coef[pos.begx : pos.endx]
  # z_pre = Z*coef[pos.begz : pos.endz]  # mu
  # q_pre = Q*coef[pos.begq : pos.endq]
    w_pre = W*coef[pos.begw : pos.endw] # log_σᵤ²
    v_pre = V*coef[pos.begv : pos.endv] # log_σᵥ²

    μ   = 0.0
    σᵤ² = @. exp(w_pre)
    σᵥ² = @. exp(v_pre)
    σ²  = σᵤ² + σᵥ²
    ϵ   = PorC*(Y - x_pre)
    μₛ  = @. (σᵥ² * μ - σᵤ² * ϵ) / σ²
    σₛ  = @. sqrt((σᵥ² * σᵤ²) / σ²)

   jlms = @. (σₛ * normpdf(μₛ / σₛ )) / normcdf(μₛ / σₛ) + μₛ
     bc = @. exp( -μₛ + 0.5 * (σₛ)^2  ) * ( normcdf( (μₛ / σₛ) - σₛ ) / normcdf(μₛ / σₛ)  )

   return jlms, bc
end

#? ------------------ Exponential model -------------

function jlmsbc(::Type{Expo}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
                Y::Matrix, X::Matrix, z, q, W::Matrix, V::Matrix, # dum1,
                dum2)

    x_pre = X*coef[pos.begx : pos.endx]
  # z_pre = Z*coef[pos.begz : pos.endz]  # mu
  # q_pre = Q*coef[pos.begq : pos.endq]
    w_pre = W*coef[pos.begw : pos.endw] # log_σᵤ²
    v_pre = V*coef[pos.begv : pos.endv] # log_σᵥ²

    λ   = @. exp(0.5*w_pre)
    σᵥ² = @. exp(v_pre)
    σᵥ  = @. exp(0.5*v_pre)

    ϵ   = PorC*(Y - x_pre)
    μₛ  = (- ϵ) -(σᵥ² ./ λ) # don't know why this line cannot use @.

    jlms = @. (σᵥ * normpdf(μₛ / σᵥ )) / normcdf(μₛ / σᵥ) + μₛ
      bc = @. exp( -μₛ + 0.5 * σᵥ²  ) * ( normcdf( (μₛ / σᵥ) - σᵥ ) / normcdf(μₛ / σᵥ)  )

    return jlms, bc
end

#? -------------- Scaling Property Model ------------

function jlmsbc(::Type{Trun_Scale}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
                Y::Matrix, X::Matrix, Z::Matrix, Q::Matrix, W::Matrix, V::Matrix, # dum1,
                dum2)

    x_pre = X*coef[pos.begx : pos.endx]
    z_pre = Z*coef[pos.begz : pos.endz]  # mu
    q_pre = Q*coef[pos.begq : pos.endq]
    w_pre = W*coef[pos.begw : pos.endw] # log_σᵤ²
    v_pre = V*coef[pos.begv : pos.endv] # log_σᵥ²

    μ   = @. (z_pre) * exp(q_pre)             # Q is _cons; μ here is the after-mutiplied-by-scaling function
    σᵤ² = @. exp(w_pre + 2*q_pre)           # The notation of σᵤ² is different from the paper. W is _cons
    σᵥ² = @. exp(v_pre)
    σ²  = σᵤ² + σᵥ²
    ϵ   = PorC*(Y - x_pre)
    μₛ  = @. (σᵥ² * μ - σᵤ² * ϵ) / σ²
    σₛ  = @. sqrt((σᵥ² * σᵤ²) / σ²)

   jlms = @. (σₛ * normpdf(μₛ / σₛ )) / normcdf(μₛ / σₛ) + μₛ
     bc = @. exp( -μₛ + 0.5 * (σₛ)^2  ) * ( normcdf( (μₛ / σₛ) - σₛ ) / normcdf(μₛ / σₛ)  )

   return jlms, bc
end

#? ---------- panel FE Wang and Ho (2010 JoE), Half normal -----------

function jlmsbc(::Type{PFEWHH}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
                Y::Vector, X::Matrix, z, Q::Matrix, W::Matrix, V::Matrix, idt::Matrix{Any})

    x_pre = X*coef[pos.begx : pos.endx]
  # z_pre = Z*coef[pos.begz : pos.endz]  # mu
    q_pre = Q*coef[pos.begq : pos.endq]
    w_pre =   coef[pos.begw] # log_σᵤ²
    v_pre =   coef[pos.begv] # log_σᵥ²

    ϵ̃   = PorC*(Y - x_pre)
    h   = exp.(q_pre)
    σᵤ² = exp(w_pre)
    σᵥ² = exp(v_pre)

    μ = 0.0

    nofobs = length(Y)
    N = size(idt)[1]

    jlms = zeros(nofobs, 1)
      bc = zeros(nofobs, 1)

    for i = 1:N
      # @views Π⁻  = INVtrM[i]*(1/σᵥ²)
        @views ind = idt[i,1]
      # @views trM = Matrix(I, idt[i,2], idt[i,2]) .- (1/idt[i,2])
        @views h̃   = sf_demean(h[ind]) # trM*h[ind]
              σₛₛ² = σₛₛ² = 1.0/(h̃'*h̃*(1/σᵥ²) + 1/σᵤ²)  # h̃ is already de-meaned, so skip Π⁻ but mind 1/σᵥ² # 1.0/(h̃'*Π⁻*h̃ + 1/σᵤ²)
              σₛₛ  = sqrt(σₛₛ²)
        @views μₛₛ = (μ/σᵤ² - ϵ̃[ind]'*h̃*(1/σᵥ²)) * σₛₛ² # (μ/σᵤ² - ϵ̃[ind]'*Π⁻*h̃) * σₛₛ²

        jlms[ind] = @. h[ind] * (μₛₛ + normpdf(μₛₛ/σₛₛ)*σₛₛ/normcdf(μₛₛ/σₛₛ))
        bc[ind]   = @. ((normcdf(μₛₛ/σₛₛ - h[ind]*σₛₛ))/normcdf(μₛₛ/σₛₛ))*exp(-h[ind]*μₛₛ + 0.5*(h[ind]^2)*σₛₛ²)
    end # for_i

    return jlms, bc
end

#? ---------- panel FE Wang and Ho (2010 JoE), truncated normal -----------

function jlmsbc(::Type{PFEWHT}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
    Y::Vector, X::Matrix, Z::Matrix, Q::Matrix, W::Matrix, V::Matrix, idt::Matrix{Any})

    x_pre = X*coef[pos.begx : pos.endx]
    z_pre =   coef[pos.begz]  # mu, a scalar
    q_pre = Q*coef[pos.begq : pos.endq]
    w_pre =   coef[pos.begw] # log_σᵤ²
    v_pre =   coef[pos.begv] # log_σᵥ²

    ϵ̃   = PorC*(Y - x_pre)
    μ   = z_pre
    h   = exp.(q_pre)
    σᵤ² = exp(w_pre)
    σᵥ² = exp(v_pre)

    nofobs = length(Y)
         N = size(idt)[1]

    jlms = zeros(nofobs, 1)
      bc = zeros(nofobs, 1)

    for i = 1:N
      # @views Π⁻  = INVtrM[i]*(1/σᵥ²)
        @views ind = idt[i,1]
      # @views trM = Matrix(I, idt[i,2], idt[i,2]) .- (1/idt[i,2])
        @views h̃   = sf_demean(h[ind]) # trM*h[ind]
              σₛₛ² = 1.0/(h̃'*h̃*(1/σᵥ²) + 1/σᵤ²)  # 1.0/(h̃'*Π⁻*h̃ + 1/σᵤ²)
              σₛₛ  = sqrt(σₛₛ²)
        @views μₛₛ = (μ/σᵤ² - ϵ̃[ind]'*h̃*(1/σᵥ²)) * σₛₛ² # (μ/σᵤ² - ϵ̃[ind]'*Π⁻*h̃) * σₛₛ²

        jlms[ind] = @. h[ind] * (μₛₛ + normpdf(μₛₛ/σₛₛ)*σₛₛ/normcdf(μₛₛ/σₛₛ))
        bc[ind]   = @. ((normcdf(μₛₛ/σₛₛ - h[ind]*σₛₛ))/normcdf(μₛₛ/σₛₛ))*exp(-h[ind]*μₛₛ + 0.5*(h[ind]^2)*σₛₛ²)
    end # for_i

    return jlms, bc
end

#? ---------- panel time decay model of Battese and Coelli (1992) -----------

function jlmsbc(::Type{PanDecay}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
  Y::Matrix, X::Matrix, Z::Matrix, Q::Matrix, W::Matrix, V::Matrix, # dum1,
  idt::Matrix{Any})

  x_pre = X*coef[pos.begx : pos.endx]
  z_pre = Z*coef[pos.begz : pos.endz]
  q_pre = Q*coef[pos.begq : pos.endq]
  w_pre =   coef[pos.begw] # log_σᵤ²
  v_pre =   coef[pos.begv] # log_σᵥ²

  nofobs = length(Y)
  N = size(idt)[1]

  σᵤ² = exp(w_pre)    # should be W*log_σᵤ² where W is a _cons =1; make a short cut here
  σᵥ² = exp(v_pre)
  σₛ² = σᵤ² + σᵥ²

  jlms_ui  = zeros(nofobs, 1)
  jlms_uit = zeros(nofobs, 1)
  bc_uit   = zeros(nofobs, 1)

  @inbounds for i = 1:N
      @views ind = idt[i,1]
      @views Tᵢ  = idt[i,2]
      @views ε   = PorC*(Y[ind] - x_pre[ind])
      @views μ   = (z_pre[ind])[1]  # ok becuase it is time-invariant
      @views Gₜ  = exp.(q_pre[ind])
      Gε  = Gₜ .* ε

      ΣGε  = sum_kbn(Gε)
      ΣGₜ² = sum_kbn(Gₜ.^2)

      μₛ   = (μ * σᵥ² - ΣGε * σᵤ²) / (σᵥ² + σᵤ² * ΣGₜ²)
      σₛ²  =  σᵥ² * σᵤ² / (σᵥ² + σᵤ² * ΣGₜ²)
      σₛ   = sqrt(σₛ²)

             jlms_ui[ind] .= μₛ + σₛ*(normpdf(μₛ/σₛ) / normcdf(μₛ/σₛ))  # identical within the panel
      @views jlms_uit[ind] = Gₜ .* jlms_ui[ind]

      bc_uit[ind] = (normcdf.(μₛ/σₛ .- Gₜ*σₛ  ) ./ normcdf( μₛ/σₛ )) .* exp.(-Gₜ.*μₛ + 0.5*((Gₜ).^2)*σₛ²)

  end

  return jlms_uit, bc_uit
end

#? ---------- panel Kumbhakar 1990 -----------

function jlmsbc(::Type{PanKumb90}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
  Y::Matrix, X::Matrix, Z::Matrix, Q::Matrix, W::Matrix, V::Matrix, # dum1,
  idt::Matrix{Any})

  x_pre = X*coef[pos.begx : pos.endx]
  z_pre = Z*coef[pos.begz : pos.endz]
  q_pre = Q*coef[pos.begq : pos.endq]
  w_pre =   coef[pos.begw] # log_σᵤ²
  v_pre =   coef[pos.begv] # log_σᵥ²

  nofobs = length(Y)
  N = size(idt)[1]

  σᵤ² = exp(w_pre)    # should be W*log_σᵤ² where W is a _cons =1; make a short cut here
  σᵥ² = exp(v_pre)
  σₛ² = σᵤ² + σᵥ²

  jlms_ui  = zeros(nofobs, 1)
  jlms_uit = zeros(nofobs, 1)
  bc_uit   = zeros(nofobs, 1)

  @inbounds for i = 1:N
      @views ind = idt[i,1]
      @views Tᵢ  = idt[i,2]
      @views ε   = PorC*(Y[ind] - x_pre[ind])
      @views μ   = (z_pre[ind])[1]  # ok becuase it is time-invariant
      @views Gₜ  = 2 ./ (1 .+ exp.(q_pre[ind])  ) # Gₜ  = exp.(q_pre[ind])
      Gε  = Gₜ .* ε

      ΣGε  = sum_kbn(Gε)
      ΣGₜ² = sum_kbn(Gₜ.^2)

      μₛ   = (μ * σᵥ² - ΣGε * σᵤ²) / (σᵥ² + σᵤ² * ΣGₜ²)
      σₛ²  =  σᵥ² * σᵤ² / (σᵥ² + σᵤ² * ΣGₜ²)
      σₛ   = sqrt(σₛ²)

             jlms_ui[ind] .= μₛ + σₛ*(normpdf(μₛ/σₛ) / normcdf(μₛ/σₛ))  # identical within the panel
      @views jlms_uit[ind] = Gₜ .* jlms_ui[ind]

      bc_uit[ind] = (normcdf.(μₛ/σₛ .- Gₜ*σₛ  ) ./ normcdf( μₛ/σₛ )) .* exp.(-Gₜ.*μₛ + 0.5*((Gₜ).^2)*σₛ²)

  end

  return jlms_uit, bc_uit
end

#? ---------- panel FE CSW (JoE 2014), Half normal -----------

function jlmsbc(::Type{PFECSWH}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
  Y::Vector, X::Matrix, z, q, w, v, idt::Matrix{Any})

   nofobs = length(Y)
   N = size(idt,1)
  #=
    for i in 1:N
        ind = idt[i,1]
        Y[ind] = Y[ind] .+ idt[i,3][1]  # add back the mean for the panel
        X[ind] = X[ind] .+ idt[i,3][2:end]
    end
  =#

    x_pre = X*coef[pos.begx : pos.endx]
    # z_pre = Z*coef[pos.begz : pos.endz]  # mu
    # q_pre = Q*coef[pos.begq : pos.endq]
    w_pre =   coef[pos.begw] # log_σᵤ²
    v_pre =   coef[pos.begv] # log_σᵥ²

    αϵ   = (Y - x_pre)  # αᵢ + ϵᵢ
    σᵤ² = exp(w_pre)
    σᵥ² = exp(v_pre)
    σᵤ = exp(0.5*w_pre)

    μ = 0.0

    jlms = zeros(nofobs, 1)
    bc = zeros(nofobs, 1)

    for i = 1:N

        @views ind = idt[i,1]
        @views αᵢᴹ = mean(αϵ[ind]) + PorC*( sqrt(2/π)*σᵤ )
        @views ϵ   = αϵ[ind] .- αᵢᴹ
               μₛ  = -(σᵤ²*ϵ) ./ (σᵤ² + σᵥ²)
               σₛ² = (σᵤ² * σᵥ²)/(σᵤ² + σᵥ²)
               r   =  μₛ ./ sqrt(σₛ²)
         jlms[ind] =  (sqrt.(σₛ²) * normpdf.(r) ./ normcdf.(r)) .+μₛ
         bc[ind]   =  exp.( -μₛ .+ 0.5 .* σₛ²  ) .* ( normcdf.( r .- sqrt.(σₛ²) ) ./ normcdf.(r)  )

    end # for_i

    return jlms, bc
end

#? ----- Panel TRE, Half Normal --------------

function jlmsbc(::Type{PTREH}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
  Y::Matrix, X::Matrix, z, q, w, v, idt::Matrix{Any})

   #! only BC, no JLMS

    nofobs = length(Y)

    β   = coef[1:pos.endx]

    σₐ² = exp(coef[pos.begq])
    σᵤ² = exp(coef[pos.begw])
    σᵥ² = exp(coef[pos.begv])

    σᵤ  = exp(0.5*coef[pos.begw])
    σᵥ  = exp(0.5*coef[pos.begv])

    ε   = PorC*(Y - X * β)

  # μ   = (z[ind, :] * δ) #! truncated normal
    μ   = zeros(nofobs, 1)

    bc  = zeros(nofobs, 1)

    nofid = size(idt, 1)

    for i in 1:nofid
        @views ind  = idt[i, 1]
        @views T    = idt[i, 2]

              Σ     = ones(T,T)*σₐ² + Matrix(σᵥ²*I,T,T)
              Σᵤ    = Matrix(σᵤ²*I,T,T)
              invΣ  = (1/σᵥ²) * (Matrix(I,T,T) - σₐ²/(σᵥ² + T*σₐ²)*ones(T,T))
              invΣᵤ = (1/σᵤ²) * Matrix(I,T,T)
        @views K    = ε[ind]'*invΣ - μ[ind]'*invΣᵤ

      # analytic form of G = inv(inv(Σ)+inv(Σᵤ))
        γ²   = 1/(1/σᵤ²+1/σᵥ²)
        g    = -T*1/(1+σₐ²/σᵥ²*T)*σₐ²/σᵥ²^2*1/(1/σᵤ²+1/σᵥ²)
        ρ²   = 1/(1+g) * 1/(1/σᵤ²+1/σᵥ²)^2 * 1/(1 + σₐ²/σᵥ²*T) *(σₐ²/σᵥ²^2)
        G    = γ² * Matrix(I,T,T) + ρ² * ones(T,T)

        B    = (-K*G)'
        γ    = sqrt(γ²)
        ρ    = sqrt(ρ²)

        MCDF1 = quadgk(u -> normpdf(u)*prod(normcdf.((B.-ρ*u)/γ)), -Inf, Inf, rtol=1e-8)

        h1 = MCDF1[1]*exp(0.5*(B'*(invΣ + invΣᵤ)*B))

        bc_i = zeros(T)
        for j in 1:T
           l_T    = zeros(T)
           l_T[j] = 1.0
           K2 = K + l_T'
           B2 = (-K2*G)'
           MCDF2 = quadgk(u -> normpdf(u)*prod(normcdf.((B2.-ρ*u)/γ)), -Inf, Inf, rtol=1e-8)
           h2 = MCDF2[1]*exp(0.5*(B2'*(invΣ + invΣᵤ)*B2))
           bc_i[j] = h2/h1
        end
           bc[ind] = bc_i

      end
           jlms = 0

    return jlms, bc
end

#? ----- Panel TRE, truncated normal ---------------

function jlmsbc(::Type{PTRET}, PorC::Int64, pos::NamedTuple, coef::Array{Float64, 1},
  Y::Matrix, X::Matrix, Z::Matrix, q, w, v, idt::Matrix{Any})

    #! only bc, no jlms

    nofobs = length(Y)

    β   = coef[1 : pos.endx]
    δ   = coef[pos.begz : pos.endz]

    σₐ² = exp(coef[pos.begq])
    σᵤ² = exp(coef[pos.begw])
    σᵥ² = exp(coef[pos.begv])

    ε   = PorC*(Y - X * β)

    μ   = (Z * δ)

    bc  = zeros(nofobs, 1)

    nofid = size(idt, 1)

    for i in 1:nofid
        @views ind  = idt[i, 1]
        @views T    = idt[i, 2]
              invΣ  = (1/σᵥ²) * (Matrix(I,T,T) - σₐ²/(σᵥ² + T*σₐ²)*ones(T,T))
              invΣᵤ = (1/σᵤ²) * Matrix(I,T,T)
        @views K     = ε[ind]'*invΣ - μ[ind]'*invΣᵤ

      # the analytic form of G = inv(inv(Σ)+inv(Σᵤ))
        γ²   = 1/(1/σᵤ²+1/σᵥ²)
        g    = -T*1/(1+σₐ²/σᵥ²*T)*σₐ²/σᵥ²^2*1/(1/σᵤ²+1/σᵥ²)
        ρ²   = 1/(1+g) * 1/(1/σᵤ²+1/σᵥ²)^2 * 1/(1 + σₐ²/σᵥ²*T) *(σₐ²/σᵥ²^2)
        G    = γ² * Matrix(I,T,T) + ρ² * ones(T,T)

        B    = (-K*G)'
        γ    = sqrt(γ²)
        ρ    = sqrt(ρ²)

        MCDF1 = quadgk(u -> normpdf(u)*prod(normcdf.((B.-ρ*u)/γ)), -Inf, Inf, rtol=1e-8)

        h1 = MCDF1[1]*exp(0.5*(B'*(invΣ + invΣᵤ)*B))

        bc_i = zeros(T)
        for j in 1:T
           l_T    = zeros(T)
           l_T[j] = 1.0
           K2 = K + l_T'
           B2 = (-K2*G)'
           MCDF2 = quadgk(u -> normpdf(u)*prod(normcdf.((B2.-ρ*u)/γ)), -Inf, Inf, rtol=1e-8)
           h2 = MCDF2[1]*exp(0.5*(B2'*(invΣ + invΣᵤ)*B2))
           bc_i[j] = h2/h1
        end
           bc[ind] = bc_i

      end
           jlms = 0

    return jlms, bc
end


text:
function gen_dr( r, center )
    Npoints = size(r)[2]
    dr = Array{Float64}(undef,Npoints)
    #
    for ip=1:Npoints
        dx2 = ( r[1,ip] - center[1] )^2
        dy2 = ( r[2,ip] - center[2] )^2
        dz2 = ( r[3,ip] - center[3] )^2
        dr[ip] = sqrt( dx2 + dy2 + dz2 )
    end
    return dr
end


text:
module QMDP

using POMDPs
using POMDPModelTools
using POMDPPolicies
using DiscreteValueIteration
using POMDPLinter: @POMDP_require, @subreq

import POMDPs: Solver
import POMDPs: solve

export
    QMDPSolver,
    solve

include("vanilla.jl")

end # module


text:
type Foo
    a
    b
    c

    version
    txn
end
Foo(a, b, c) = Foo(a, b, c, nothing, nothing)

FooTable = DynamoTable(Foo, string_key("a"), real_key("b"),
                       reads=10, writes=10,
                       version_column="version", txn_column="txn")
FooLSI = DynamoLocalIndex(FooTable, real_key("c"), reads=10, writes=10)
FooGSI = DynamoGlobalIndex(FooTable, real_key("c"), reads=10, writes=10)
create_table(FooTable, FooLSI, FooGSI)

load_foo(args...) = get_item(FooTable, args...)
load_foo_lsi(args...) = get_item(FooLSI, args...)
load_foos_by_c(args...) = get_items(FooGSI, args...)
persist(f :: Foo, args...) = put_item(FooTable, f, args...)

persist(Foo("one", 2, 3))

load_foo("one", 2)
load_foo_lsi("one", 3)
persist(Foo("one", 2, 3)) # boom ... missing version

persist(Foo("two", 3, 4))

item = load_foo("one", 2)
item.b = 20
persist(item)

dynamo_transaction(retries=3) do |t|
    one = load_foo("one", 2)
    two = load_foo("two", 3)

    # we want an atomic write on both
    one.b += 1
    two.b -= 1

    persist(t, one)
    persist(t, two)

    rollback(t) # undooooo!
    error("blah") # same as rollback

    commit(t) # writes it
end

f = load_foo("one", 2)
f.c += 1
persist(f, (attr("c") >= 10) && size(attr("a")) <= 4)


text:
"""
    struct MultiFieldFEFunction <: CellField
      # private fields
    end
"""
struct MultiFieldFEFunction{T<:MultiFieldCellField} <: CellField
  single_fe_functions::Vector{<:SingleFieldFEFunction}
  free_values::AbstractArray
  cell_dof_values::AbstractArray
  fe_space::MultiFieldFESpace
  multi_cell_field::T

  function MultiFieldFEFunction(
    free_values::AbstractVector,
    space::MultiFieldFESpace,
    single_fe_functions::Vector{<:SingleFieldFEFunction})

    multi_cell_field = MultiFieldCellField(map(i->i.cell_field,single_fe_functions))
    T = typeof(multi_cell_field)

    f(i) = get_cell_axes(get_fe_space(i))
    cell_axes = create_array_of_blocked_axes(map(f,single_fe_functions)...)
    blocks = Tuple(map(get_cell_values,single_fe_functions))
    blockids = [(i,) for i in 1:length(single_fe_functions)]
    cell_dof_values = VectorOfBlockArrayCoo(blocks,blockids,cell_axes)

    new{T}(
      single_fe_functions,
      free_values,
      cell_dof_values,
      space,
      multi_cell_field)
  end
end

#function MultiFieldFEFunction(
#  space::MultiFieldFESpace,
#  blocks::Vector{<:SingleFieldFEFunction})
#  fv = zero_free_values(space)
#  xh0 = MultiFieldFEFunction(fv,space,blocks)
#end

Arrays.get_array(a::MultiFieldFEFunction) = get_array(a.multi_cell_field)

CellData.get_memo(a::MultiFieldFEFunction) = get_memo(a.multi_cell_field)

CellData.get_cell_map(a::MultiFieldFEFunction) = get_cell_map(a.multi_cell_field)

CellData.get_cell_axes(a::MultiFieldFEFunction) = get_cell_axes(a.multi_cell_field)

CellData.RefStyle(::Type{<:MultiFieldFEFunction{T}}) where T = RefStyle(T)

CellData.MetaSizeStyle(::Type{<:MultiFieldFEFunction{T}}) where T = MetaSizeStyle(T)

FESpaces.FEFunctionStyle(::Type{<:MultiFieldFEFunction}) = Val{true}()

FESpaces.get_free_values(f::MultiFieldFEFunction) = f.free_values

FESpaces.get_fe_space(f::MultiFieldFEFunction) = f.fe_space

Base.length(f::MultiFieldFEFunction) = length(f.multi_cell_field)

"""
    num_fields(m::MultiFieldFEFunction)
"""
num_fields(m::MultiFieldFEFunction) = length(m.single_fe_functions)

Base.iterate(m::MultiFieldFEFunction) = iterate(m.single_fe_functions)

Base.iterate(m::MultiFieldFEFunction,state) = iterate(m.single_fe_functions,state)

Base.getindex(m::MultiFieldFEFunction,field_id::Integer) = m.single_fe_functions[field_id]

function Geometry.restrict(a::MultiFieldFEFunction,trian::Triangulation)
  f = (ai) -> restrict(ai,trian)
  blocks = map(f,a.single_fe_functions)
  blocks
end

function FESpaces.get_cell_values(f::MultiFieldFEFunction)
  f.cell_dof_values
end


text:
# Test utililties.jl

# Tests could be more comprehensive (of course).  In some cases all
# that is done is to execute the function in some way just to make
# sure no exceptions are thrown as a result of something breaking with
# a new version of Julia.

println("testing utilites")

# bbspline
P = [0 1 2 3
     0 1 2 3]
k = 3
N = 10
S = bbspline(P, k, N)
dx = gradient(vec(S[1,:]))
dy = gradient(vec(S[2,:]))
slope = dy./dx
@test maximum(slope - 1) < 1e-3

# pbspline
P = [0 1 2 3 0
     0 1 2 3 0]
k = 3
N = 100
S = pbspline(P, k, N)

# gaussfilt1d
a = zeros(20,1)
a[10] = 1
sigma = 2
sa = gaussfilt1d(a, sigma)
@test abs(sa[10] - 0.282) < 1e-3
@test abs(sum(sa) - 1) < 1e-3

# interp1
x = [0, 2, 7]
y = 2*x
xi = [1.1, 0.5, 6]
yi = interp1(x,y,xi)
@test maximum(abs(yi - 2*xi) .< eps())

# normalise
a = rand(5,3)
an = normalise(a)
@test isapprox(minimum(an), 0)
@test isapprox(maximum(an), 1)

reqmean = 5
reqvar = 7
an = normalise(a, reqmean, reqvar)
@test isapprox(mean(an), reqmean)
@test isapprox(var(an), reqvar)

# histtruncate
rows = 100;
cols = 200;
img = rand(rows, cols)
lHistCut = 2
uHistCut = 4
htimg = histtruncate(img, lHistCut, uHistCut)
minval = minimum(htimg)
maxval = maximum(htimg)
# Check that the number of saturated pixels at each extreme is
# approximately lHistCut and uHistCut. Note some integer rounding has
# to occur.
v = find(abs(htimg-minval).<eps())
@test abs(length(v)/(rows*cols) * 100 - lHistCut) < 1

v = find(abs(htimg-maxval).<eps())
@test abs(length(v)/(rows*cols) * 100 - uHistCut) < 1

# sineramp
sr = sineramp()
sr = sineramp((256, 512), 12.5, 8, 2)

# circlesineramp
(img, alpha) = circlesineramp()
(img, alpha) = circlesineramp(512, pi/10, 8, 2, false)

# meshgrid
xrange = -3:2
yrange = 2:.5:6
(x,y) = meshgrid(xrange, yrange) # two arguments
(rows,cols) = size(x)
for r = 1:rows
    @test all(x[r:r,:] .== xrange')
end
for c = 1:cols
    @test all(y[:,c:c] .== yrange)
end

arange = 2:.1:4
(x,y) = meshgrid(arange)  # one argument
(rows,cols) = size(x)
for r = 1:rows
    @test all(x[r:r,:] .== arange')
end
for c = 1:cols
    @test all(y[:,c:c] .== arange)
end


text:
#  Copyright 2016, 2017 Eric S. Tellez <eric.tellez@infotec.mx>
#
#    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.

# abstract Sequential

import Base: push!

export Sequential, search, push!, fit, optimize!

"""
    Sequential{T}

A simple exhaustive search index
"""
struct Sequential{T} <: Index
    db::Vector{T}
end

function fit(::Type{Sequential}, db::Vector{T}) where T
    Sequential(db)
end

function search(index::Sequential{T}, dist::Function, q::T, res::KnnResult) where T
    i::Int32 = 1
    d::Float64 = 0.0
    for obj in index.db
        d = dist(q, obj)
        push!(res, i, convert(Float32, d))
        i += 1
    end

    res
end

function push!(index::Sequential{T}, dist::Function, item::T) where T
    push!(index.db, item)
    length(index.db)
end

"""
    optimize!(index::Sequential{T}, dist::Function, recall::Float64)

Optimizes the index for the required quality
"""
function optimize!(index::Sequential{T}, dist::Function, recall::Float64) where T
    # do nothing for sequential
end


text:
export girth, girth_cycle

"""
`girth_cycle(G)` returns a shorest cycle of `G` as an array listing
the vertices on that cycle, or an empty array if `G` is acyclic.

*Warning*: This implementation is quite inefficient.
"""
function girth_cycle(G::SimpleGraph{T}) where {T}
  if cache_check(G,:girth_cycle)
    return cache_recall(G,:girth_cycle)
  end
  best_path = T[]
  if is_acyclic(G)
    cache_save(G,:girth_cycle,best_path)
    return best_path
  end
  best = NV(G)+1

  GG = deepcopy(G)

  for e in elist(GG)
    u,v = e
    SimpleGraphs.delete!(GG,u,v)
    P = find_path(GG,u,v)
    add!(GG,u,v)
    nP = length(P)
    if  0 < nP < best
      best = nP
      best_path = P
    end
  end
  cache_save(G,:girth_cycle,best_path)
  return best_path
end

"""
`girth(G)` computes the length of a shortest cycle in `G` or returns `0`
if `G` is acyclic.

**Warning**: This implementation is quite inefficient.
"""
function girth(G::SimpleGraph)
  if cache_check(G,:girth)
    return cache_recall(G,:girth)
  end
  g = length(girth_cycle(G))
  cache_save(G,:girth,g)
  return g
end


text:
function math_expr(op, op1, op2)
    expr = Expr(:call, op, op1, op2)
    return expr
end

ex = math_expr(:+, 1, Expr(:call, :*, 4, 5))

eval(ex)

function make_expr2(op, opr1, opr2)
    opr1f, opr2f = map(x -> isa(x, Number) ? 2*x : x, (opr1, opr2))
    retexpr = Expr(:call, op, opr1f, opr2f)
    return retexpr
end

make_expr2(:+, 1, 2)

ex = make_expr2(:+, 1, Expr(:call, :*, 5, 8))

eval(ex)

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP


text:
"""
SamplerState for the multidimensional stochastic gradient Nosé-Hoover thermostat algorithm following

Ding, N., Fang, Y., Babbush, R., Chen, C., Skeel, R. D., & Neven, H. (2014). Bayesian Sampling Using Stochastic Gradient Thermostats. In Advances in Neural Information Processing Systems (pp. 3203–3211).

and

Leimkuhler, B., & Shang, X. (2016). Adaptive Thermostats for Noisy Gradient Systems. SIAM Journal on Scientific Computing, 38(2), A712–A736.

In contrast to SGNHT this sampler has a thermostat variable per dimension.

Fields:

    - x: the current state
    - p: current momentum
    - xi: thermostat variable (initialised to be D*ones)
    - t: iterations counter
    - stepsize: s.stepsize(s.t) gives the stepsize to be used in iteration s.t
    - mu: Constant in the thermostat equation
    - D: friction parameters
Constructors:
    - construct using `mSGNHTState(x,p,stepsize::Function;mu=1.,D=1.0)` or `mSGNHTState(x,p,stepsize::Float64;mu=1.,D=1.0)` for a fixed stepsize if p is not specified then          p = randn(length(x))

"""
type mSGNHTState <: SamplerState
    x::Array{Float64,1}
    p::Array{Float64,1}
    xi::Array{Float64,1}

    t::Int
    stepsize::Function
    mu::Float64
    D::Float64
    function mSGNHTState(x,p,stepsize::Function;mu=1.,D=1.0)
        @assert size(x) = size(p)
        xi = ones(length(x))*D
        new(x,p,xi,0,stepsize,mu,D)
    end
    function mSGNHTState(x,p,stepsize::Float64;mu=1.,D=1.0)
        @assert size(x) = size(p)
        xi = ones(length(x))*D
        s(niters) = stepsize
        new(x,p,xi,0,s,mu,D)
    end
    function mSGNHTState(x,stepsize::Function;mu=1.,D=1.0)
        p = randn(length(x))
        xi = ones(length(x))*D
        new(x,p,xi,0,stepsize,mu,D)
    end
    function mSGNHTState(x,stepsize::Float64;mu=1.,D=1.0)
        p = randn(length(x))
        xi = ones(length(x))*D
        s(niters) = stepsize
        new(x,p,xi,0,s,mu,D)
    end

end

"""
sample! performs one mSGNHTState update and returns the updated mSGNHTState.
"""
function sample!(s::mSGNHTState,grad::Function)
    # increment iterations counter
    s.t += 1

    # improve readability
    nparams = length(s.x)
    stepsize = s.stepsize(s.t)
    x = s.x
    p = s.p
    xi = s.xi
    D = s.D
    mu = s.mu

    # second order symmetric splitting scheme from Leimkuhler et al.
    add_inplace!(.5*stepsize,grad(x),p)
    add_inplace!(0.5*stepsize,p,x)
    add_inplace!(.5*stepsize/mu, (p.*p - 1),xi)
    ids = abs(stepsize*xi) .> 1e-6
    p[ids] = exp(-stepsize*xi[ids]).*p[ids] + sqrt(.5.*(1-exp(-2.0*stepsize*xi[ids]))./xi[ids] * D) .* randn(sum(ids))
    p[!ids] += sqrt(stepsize*D)*randn(sum(!ids))
    add_inplace!(.5*stepsize/mu, (p.*p - 1),xi)
    add_inplace!(.5*stepsize,p,x)
    add_inplace!(.5*stepsize,grad(x),p)
    s
end


text:
"""
    PoincareInvariants

A Julia library for the computation of Poincaré integral invariants.
"""
module PoincareInvariants

export AbstractPoincareInvariant, compute!
export FirstPoincareInvariant, CanonicalFirstPI, FirstPI
export SecondPoincareInvariant, CanonicalSecondPI, SecondPI

export getpoints, getpointspec, getpointnum
export getdim, getform, getplan

export FirstFinDiffPlan, FirstFourierPlan
export SecondChebyshevPlan, SecondFinDiffPlan

export canonical_one_form, CanonicalSymplecticMatrix, canonical_two_form

include("utils.jl")

## AbstractPoincareInvariant ##

"""
    AbstractPoincareInvariant{T, D}

represents a Poincare integral invariant in a phase space of dimension D using numeric type
T for calculations.
"""
abstract type AbstractPoincareInvariant{T, D} end

# these types will be accepted as a point cloud representing an area or line
# Any other type will get fed to a PoincareInvariantIter constructor
# The rows of an AbstractMatrix are points and
# the elements of the AbstractVector{AbstractVector} are points
const PointCollection{T} = Union{AbstractMatrix{T}, AbstractVector{AbstractVector{T}}}

checkpoints(T::Type) = @argcheck T <: PointCollection "Points must be a matrix or vector of vectors"

writepoints!(pinv::AbstractPoincareInvariant, points::AbstractMatrix) = (pinv.points .= points)
function writepoints!(pinv::AbstractPoincareInvariant, points::AbstractVector{AbstractVector})
    for i in 1:length(points)
        pinv.points[:, i] .= points[i]
    end
end

"""
    compute!(pinv::AbstractPoincareInvariant, args...)

computes a Poincaré invariant.

implementations should define a method `compute!(pinv, t, p)`, which acts on the internal
points storage.
"""
function compute!(
    pinv::AbstractPoincareInvariant,
    points::PointCollection,
    t::Real, p
)
    writepoints!(pinv, points)
    compute!(pinv, t, p)
end

function compute!(pinv::AbstractPoincareInvariant, data, times::AbstractVector, p)
    _compute_iter(pinv, data, times, p, Base.IteratorEltype(data))
end

function _compute_iter(pinv, data, times, p, ::Base.HasEltype)
    checkpoints(eltype(data))
    map(enumerate(data)) do (i, points)
        writepoints!(pinv, points)
        compute!(pinv, times[i], p)
    end
end

function _compute_iter(pinv, data, times, p, ::Base.EltypeUnknown)
    map(enumerate(data)) do (i, points)
        checkpoints(typeof(points))
        writepoints!(pinv, points)
        compute!(pinv, times[i], p)
    end
end

#=
function compute!(picache::ThreadCache{<:AbstractPoincareInvariant{T}}, data, ts, p) where T
    n = datalength(data)
    lk = ReentrantLock()

    out = Vector{T}(undef, n)
    @threads for i in 1:n
        # each thread gets one pinv
        pinv = picache[Threads.threadid()]

        # lock in case data has mutable state
        lock(_ -> writepoints!(pinv, data, i), lk)

        out[i] = compute!(pinv, ts[i], p)
    end

    out
end

function compute!(pinv::AbstractPoincareInvariant, data, t, p)
    iter = PoincareInvariantIter(pinv, pointsiter, args...)
    map(iter) do (points, t, p)
        compute!(pinv, points, t, p)
    end
end
=#

"""
    getplan(pinv::AbstractPoincareInvariant)

returns plan that will be used to `compute!` `pinv`.
"""
function getplan end

"""
    getpointspec(pinv::AbstractPoincareInvariant)

get point specification, which may, for example, be a tuple specifying a grid or
a number giving the number of points used to sample in phase space.
"""
function getpointspec end

"""
    getpoints(pinv::AbstractPoincareInvariant)

returns points on which to evaluate the phase space line or surface parameterisation
so as to `compute!` `pinv`.
"""
function getpoints(f::Function, pinv::AbstractPoincareInvariant{T, D}) where {T, D}
    getpoints(f, T, getpointspec(pinv), typeof(getplan(pinv)))
end

"""
    getpointnum(pinv::AbstractPoincareInvariant)

returns number of points to sample in phase space to `compute!` `pinv`.
"""
getpointnum(pinv::AbstractPoincareInvariant) = getpointnum(getpointspec(pinv))
getpointnum(N::Integer) = N

"""
    getdim(pinv::AbstractPoincareInvariant)

returns dimension of phase space to `compute!` `pinv` in.
"""
function getdim end

"""
    getform(pinv::AbstractPoincareInvariant)

get invariant one- or two-form.
"""
function getform end

## Canonical Forms ##

include("CanonicalSymplecticForms.jl")

using .CanonicalSymplecticForms: canonical_one_form, CanonicalSymplecticMatrix,
    canonical_two_form

## FirstPoincareInvariant ##

struct FirstPoincareInvariant{T, D, θT, P} <: AbstractPoincareInvariant{T, D}
    θ::θT
    N::Int
    plan::P
    points::Matrix{T}
end

function FirstPoincareInvariant{T, D}(θ::θT, N::Integer, plan::P) where {T, D, θT, P}
    n = getpointspec(N, P) |> getpointnum
    points = Matrix{T}(undef, n, D)
    FirstPoincareInvariant{T, D, θT, P}(θ, n, plan, points)
end

function FirstPoincareInvariant{T, D}(
    θ::θT, N::Integer, P::Type=DEFAULT_FIRST_PLAN
) where {T, D, θT}
    plan = P{T, D}(θ, getpointspec(N, P))
    FirstPoincareInvariant{T, D}(θ, N, plan)
end

function FirstPoincareInvariant{T, D, typeof(canonical_one_form)}(
    N, P=DEFAULT_FIRST_PLAN
) where {T, D}
    FirstPoincareInvariant{T, D}(canonical_one_form, N, P)
end

const FirstPI = FirstPoincareInvariant
const CanonicalFirstPI{T, D} = FirstPoincareInvariant{T, D, typeof(canonical_one_form)}

# Interface

getpointspec(pinv::FirstPoincareInvariant) = pinv.N
getpoints(pinv::FirstPoincareInvariant) = getpoints(identity, pinv)
getform(pinv::FirstPoincareInvariant) = pinv.θ
getplan(pinv::FirstPoincareInvariant) = pinv.plan
getdim(::FirstPoincareInvariant{<:Any, D}) where D = D

# Implementations

include("FirstFinDiffPlans.jl")
include("FirstFourierPlans.jl")

using .FirstFinDiffPlans: FirstFinDiffPlan
using .FirstFourierPlans: FirstFourierPlan

const DEFAULT_FIRST_PLAN = FirstFourierPlan

## SecondPoincareInvariant ##

struct SecondPoincareInvariant{
    T,  # phase space and return type
    D,  # dimension of phase space
    ωT,
    PS,
    P
} <: AbstractPoincareInvariant{T, D}
    ω::ωT  # symplectic matrix or function returning one
    pointspec::PS  # specifies how many points
    plan::P  # plan for chebyshev transform, differentiation, etc...
    points::Matrix{T}
end

function SecondPoincareInvariant{T, D}(ω::ωT, N, plan::P) where {T, D, ωT, P}
    ps = getpointspec(N, P)
    points = Matrix{T}(undef, getpointnum(ps), D)
    SecondPoincareInvariant{T, D, ωT, typeof(ps), P}(ω, ps, plan, points)
end

function SecondPoincareInvariant{T, D}(
    ω, N, P::Type=DEFAULT_SECOND_PLAN
) where {T, D}
    plan = P{T, D}(ω, getpointspec(N, P))
    SecondPoincareInvariant{T, D}(ω, N, plan)
end

function SecondPoincareInvariant{T, D, CanonicalSymplecticMatrix{T}}(
    N, P::Type=DEFAULT_SECOND_PLAN
) where {T, D}
    ω = CanonicalSymplecticMatrix{T}(D)
    SecondPoincareInvariant{T, D}(ω, N, P)
end

const SecondPI = SecondPoincareInvariant
const CanonicalSecondPI{T, D} = SecondPoincareInvariant{T, D, CanonicalSymplecticMatrix{T}}

# Interface

getdim(::SecondPoincareInvariant{<:Any, D, <:Any, <:Any, <:Any}) where D = D
getform(pinv::SecondPoincareInvariant) = pinv.ω
getplan(pinv::SecondPoincareInvariant) = pinv.plan

getpoints(pinv::SecondPoincareInvariant) = getpoints((x, y) -> (x, y), pinv)

getpointspec(pinv::SecondPoincareInvariant) = pinv.pointspec

# Implementations

include("SecondChebyshevPlans.jl")
include("SecondFinDiffPlans.jl")

using .SecondChebyshevPlans: SecondChebyshevPlan
using .SecondFinDiffPlans: SecondFinDiffPlan

const DEFAULT_SECOND_PLAN = SecondChebyshevPlan

end  # module


text:
"""
    Object

Defines what is drawn in a defined frame range.

# Fields
- `frames::Frames`: A range of frames for which the `Object` is called
- `func::Function`: The drawing function which draws something on the canvas.
    It gets called with the arguments `video, object, frame`
- `start_pos::Union{Object, Point}` defines the origin of the object. It gets translated to this point
- `actions::Vector{AbstractAction}` a list of actions applied to this object
- `current_setting`:: The current state of the object see [`ObjectSetting`](@ref)
- `opts::Any` can hold any options defined by the user
- `change_keywords::Dict{Symbol,Any}` the modified keywords changed by `change`
- `result::Vector` the result of the object (if something gets returned)
"""
struct Object <: AbstractObject
    frames::Frames
    func::Function
    start_pos::Union{Object,Point}
    actions::Vector{AbstractAction}
    current_setting::ObjectSetting
    opts::Dict{Symbol,Any}
    change_keywords::Dict{Symbol,Any}
    result::Vector
end

"""
    CURRENT_OBJECT

holds the current object in an array to be declared as a constant
The current object can be accessed using CURRENT_OBJECT[1]
"""
const CURRENT_OBJECT = Array{Object,1}()

Object(func::Function, args...; kwargs...) = Object(:same, func, args...; kwargs...)

Object(frames, func::Function; kwargs...) = Object(frames, func, O; kwargs...)

"""
    Object([frames], func::Function, [start_pos]; kwargs...)

# Arguments
- frames can be a `Symbol`, a `UnitRange` or a relative way to define frames see [`RFrames`](@ref)
    - **Default:** If not defined it will be the same as the previous [`Object`](@ref).
    - **Important:** The first `Object` needs the frames specified as a `UnitRange`.
    - It defines for which frames the object is active
- func is a `Function` and the only required argument
    - This defines the actual object that gets drawn.
    - The function takes the following three arguments:
        - video, object, frame
    - If you don't need them you can write `(args...)->your_function(arg1, arg2)`

# Example
```julia
function ground(args...)
    background("black")
    sethue("white")
end

video = Video(500, 500)
Background(1:100, ground)
Object((args...)->circle(O, 50, :fill))
render(video; pathname="test.gif")
```

Here the [`Background`](@ref) uses the named way of defining the function whereas
the circle object is defined in the anonymous function `(args...)->circle(O, 50, :fill)`.
It basically depends whether you want to have a simple Luxor object or something more complex.
"""
function Object(frames, func::Function, start_pos::Union{Object,Point}; kwargs...)
    if isempty(CURRENT_VIDEO)
        throw(ErrorException("A `Video` must be defined before an `Object`"))
    end
    CURRENT_VIDEO[1].defs[:last_frames] = frames
    opts = Dict(kwargs...)

    if get(opts, :in_global_layer, false) && frames isa UnitRange
        CURRENT_VIDEO[1].background_frames =
            union(CURRENT_VIDEO[1].background_frames, frames)
    end
    object = Object(
        frames,
        func,
        start_pos,
        AbstractAction[],
        ObjectSetting(),
        opts,
        Dict{Symbol,Any}(),
        Any[nothing],
    )
    push!(CURRENT_VIDEO[1].objects, object)
    return object
end

"""
    act!

Adds an [`Action`] or a list of actions to an [`Object`](@ref) or a list of objects.

# Example
```julia
Background(1:100, ground)
obj = Object((args...) -> rect(O, 50, 50, :fill), Point(100, 0))
act!(obj, Action(1:50, anim_scale(1.5)))
```

Here the scaling is applied to the rectangle for the first fifty frames.

# Options
1. A single object and action:

   `act!(object::AbstractObject, action::AbstractAction)`
   - `object::AbstractObject` - the object the action is applied to
   - `action::AbstractAction` - the action applied to the object

2. A single object and a list of actions:

   `act!(object::AbstractObject, action::Vector{<:AbstractAction})`
   - `object::AbstractObject` - the object actions are applied to
   - `action::Vector{<:AbstractAction}` - the actions applied to an object

3. A list of objects and a list of actions:

   `act!(object::Vector{<:AbstractObject}, action::Vector{<:AbstractAction})`
   - `object::Vector{<:AbstractObject}` - the objects actions are applied to
   - `action::Vector{<:AbstractAction}` - the actions applied to the objects
"""
function act!(object::AbstractObject, action::AbstractAction)
    push!(object.actions, copy(action))
end

function act!(object::AbstractObject, actions::Vector{<:AbstractAction})
    for action in actions
        act!(object, action)
    end
end

function act!(objects::Vector{<:AbstractObject}, action)
    for object in objects
        act!(object, action)
    end
end

"""
    Background(frames, func)

The `Background` is internally just an [`Object`](@ref) and can be defined the same way.
In contrast to an object this a `Background` will change the global canvas and not just
a layer. Normally it's used to define defaults and the `background` color. See Luxor.background

# Example
```julia
function ground(args...)
    background("black")
    sethue("white")
end

video = Video(500, 500)
Background(1:100, ground)
Object((args...)->circle(O, 50, :fill))
render(video; pathname="test.gif")
```

This draws a white circle on a black background as `sethue` is defined for the global frame.
"""
function Background(frames, func::Function, args...; kwargs...)
    Object(frames, func, args...; in_global_layer = true, kwargs...)
end


text:
function _call end

module NoJuliaInterpreter

module JuliaInterpreter
macro bp() end
end

import ..Recalls: _call
include("juliainterpreter.jl")

end


text:
module FFEstimator

using Unitful
using DataStructures

include("./io.jl")

export FFreqEstimator

# Autocorellation function of lag
function r(samples::AbstractVector{T}, lag::Int, windowSize::Int) where {T <: Number}
    return sum(@view(samples[1:windowSize]) .* @view(samples[lag+1:lag+windowSize]))
end

# Struct for optimizing `d` function execution time (iterative sum)
mutable struct DCache{T <: Number}
    prevSuffixSum::Vector{T} # terms that are in both previous and current sum
    # lag range from last iteration for removal of stale values
    lastStart::Int
    lastEnd::Int

    DCache{T}() where {T <: Number} = new{T}([], 0, 0)
end

# Stride cache one step function
function strideCache(d_cache::DCache{T}) where {T <: Number}
    d_cache.prevSuffixSum[1:d_cache.lastStart-1] .= T(-1)
    resize!(d_cache.prevSuffixSum, d_cache.lastEnd)
    d_cache.lastStart = d_cache.lastEnd = 0
end

# Difference function of lag (optimized)
function d(d_cache::DCache{T}, samples::AbstractVector{T}, lag::Int, windowSize::Int, windowStride::Int) where {T <: Number}
    if d_cache.lastStart == 0
        d_cache.lastStart = lag
    end
    d_cache.lastEnd = lag

    func = (b, e) -> sum((@view(samples[b:e]) .- @view(samples[lag+b:lag+e])) .^ 2)
    if windowStride * 3 <= windowSize
        # window stride is small compared to window size, can optimize
        prefix = func(1, windowStride)
        if length(d_cache.prevSuffixSum) < lag
            push!(d_cache.prevSuffixSum, -1)
        end
        if d_cache.prevSuffixSum[lag] < 0
            # previous value is not evaluated => evaluate everything now
            res = prefix + func(windowStride+1, windowSize)
        else
            # use previous value
            suffix = func(windowSize - windowStride + 1, windowSize)
            res = d_cache.prevSuffixSum[lag] + suffix
        end
        # save for future generations
        d_cache.prevSuffixSum[lag] = res - prefix
        return res
    else
        # optimization is not worth it
        return func(1, windowSize)
    end
end

# Cumulative mean normalized difference function
function d_prime(d_cache::DCache{T}, samples::AbstractVector{T}, prefixSum::T, lag::Int, windowSize::Int, windowStride::Int) where {T <: Number}
    if lag == 0
        return (T(1), T(0))
    else
        d_lag = d(d_cache, samples, lag, windowSize, windowStride)
        return (d_lag / (1 / lag * (prefixSum + d_lag)), prefixSum + d_lag)
    end
end

# Make argmin, argmax work for generators
Base.keys(g::Base.Generator) = g.iter

# Autocorellation function peak
function ACFMethod(samples::AbstractVector{T}, startLag::Int, endLag::Int, windowSize::Int) where {T <: Number}
    return argmax(r(samples, lag, windowSize) for lag in startLag:endLag)
end

# Difference function dip
function DFMethod(d_cache::DCache{T}, samples::AbstractVector{T}, startLag::Int, endLag::Int, windowSize::Int, windowStride::Int) where {T <: Number}
    return argmin(d(d_cache, samples, lag, windowSize, windowStride) for lag in startLag:endLag)
end

# Cumulative mean normalized difference function dip
function CMNDFMethod(d_cache::DCache{T}, samples::AbstractVector{T}, startLag::Int, endLag::Int, windowSize::Int, windowStride::Int) where {T <: Number}
    prefixSum = sum(d(d_cache, samples, lag, windowSize, windowStride) for lag in 1:startLag-1)
    foundThreshold = false
    foundThresholdDip = false
    lastRes, minRes = typemax(T), typemax(T)
    toret = argmin(
        let
            if foundThresholdDip
                # evade too high errors and drop out
                typemax(T)
            else
                res, prefixSum = d_prime(d_cache, samples, prefixSum, lag, windowSize, windowStride)
                if res < minRes
                    minRes = res
                end
                if res < 0.1
                    # found acceptable value already, wait for minimum
                    foundThreshold = true
                end
                if foundThreshold && !foundThresholdDip && res > lastRes
                    # already found minimum after acceptable value, drop out
                    foundThresholdDip = true
                    res = typemax(T)
                end
                lastRes = res
                res
            end
        end
        for lag in startLag:endLag)
    return (toret, minRes)
end

# Fundamental frequency estimator struct (with data related to current audio stream)
struct FFreqEstimator{InT <: Number, EvalT <: Number}
    inFile
    inStreamFormat::streamFormat{InT}

    # fundamental period search range
    startLag::Int
    endLag::Int
    # window data
    windowSize::Int
    windowStride::Int
    # radius for best local estimate search
    shoppingRadius::Int

    # samples buffer
    buf::CircularBuffer{EvalT}

    function FFreqEstimator{EvalT}(inFile, inStreamFormat::streamFormat{InT},
                                   startFreq::Unitful.Frequency, endFreq::Unitful.Frequency,
                                   windowSize::Int, windowStride::Int,
                                   filterDuration::Unitful.Time) where {InT <: Number, EvalT <: Number}
        startLag, endLag, shoppingRadius = let
            noUnitsToInt = x -> trunc(Int, convert(Float64, x))

            min(noUnitsToInt(inStreamFormat.samplerate / endFreq), windowSize),
            min(noUnitsToInt(inStreamFormat.samplerate / startFreq), windowSize),
            noUnitsToInt((filterDuration / 2 * inStreamFormat.samplerate) / windowSize)
        end

        new{InT, EvalT}(inFile, inStreamFormat, startLag, endLag,
                        windowSize, windowStride, shoppingRadius,
                        CircularBuffer{EvalT}(2 * windowSize + shoppingRadius * windowStride))
    end
end

# Estimate fundamental period with best method function
function estimate(data::FFreqEstimator{InT, EvalT}, d_cache::DCache{EvalT}, offset::Int, startLag::Int = data.startLag, endLag::Int = data.endLag) where {InT <: Number, EvalT <: Number}
    res = CMNDFMethod(d_cache, @view(data.buf[offset+1:offset+2*data.windowSize]), startLag, endLag, data.windowSize, data.windowStride)
    strideCache(d_cache)
    return res
end

# Interface for fundamental frequency estimation (begin iteration)
function Base.iterate(data::FFreqEstimator{InT, EvalT}) where {InT <: Number, EvalT <: Number}
    shoppedEsts = Deque{Tuple{Int, EvalT}}()
    Cache1, Cache2 = DCache{EvalT}(), DCache{EvalT}()
    try
        x = normalizedChannelRead(data.inFile, data.inStreamFormat.endianBOM, data.inStreamFormat.nchannels,
                                  InT, EvalT, 2 * data.windowSize + data.shoppingRadius * data.windowStride)
        append!(data.buf, x)

        for _ in 1:data.shoppingRadius
            push!(shoppedEsts, (0, typemax(InT)))
        end
        for i in 0:data.shoppingRadius
            push!(shoppedEsts, estimate(data, Cache1, i*data.windowStride))
        end
    catch e
        if !isa(e, EOFError)
            throw(e)
        end
        return nothing
    end
    return Base.iterate(data, (0, 0, shoppedEsts, Cache1, Cache2, 400))
end

# Interface for fundamental frequency estimation (continue iteration)
function Base.iterate(data::FFreqEstimator{InT}, (samplePos, timeStamp, shoppedEsts, Cache1, Cache2, currentPeriod)) where {InT <: Number}
    # get predicted period from around here
    (predictedPeriod, predictedVal) = argmin(val_ for (period_, val_) in shoppedEsts)

    # Method from paper does not work really well with music, so just select best period from local range
    # (period, val) = estimate(data, Cache2, 0,
    #                          trunc(Int, max(predictedPeriod * 0.8, data.startLag)),
    #                          trunc(Int, min(predictedPeriod * 1.2, data.endLag)))
    currentPeriod = predictedPeriod
    freq = ustrip((data.inStreamFormat.samplerate / currentPeriod) |> u"Hz")

    # move window
    for _ in 1:data.windowStride popfirst!(data.buf) end
    x = nothing
    try
        x = channelRead(data.inFile, data.inStreamFormat.endianBOM, data.inStreamFormat.nchannels, InT, data.windowStride)
    catch e
        if !isa(e, EOFError)
            throw(e)
        end
        return nothing
    end
    append!(data.buf, x)

    popfirst!(shoppedEsts)
    push!(shoppedEsts, estimate(data, Cache1, data.shoppingRadius*data.windowStride))

    certainty = let
        from, to = 0.75, 1.25
        clamp((to - predictedVal) / (to - from), 0, 1)
    end

    return ((currentPeriod, freq, samplePos, timeStamp, certainty),
            (
                samplePos += data.windowStride,
                timeStamp += ustrip((data.windowStride / data.inStreamFormat.samplerate) |> u"ms"),
                shoppedEsts, Cache1, Cache2, currentPeriod
            ))
end

end # module


text:
using ConvDiffMIPDECO
using jInv.Mesh
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.LinearSolvers
using MUMPSjInv
using MAT
# using jInvVis
using LinearAlgebra

# read results from l-curve
dataset    = "Peaks2D"
reg        = "wTVReg"
m          = [256 128]
noiseLevel = 0.1

filename= joinpath(dirname(pathof(ConvDiffMIPDECO)),"..","examples","Peaks2D.mat")
data = matread(filename)
domain = data["domain"]
v      = data["v"]
sig    = data["sig"]
dtrue  = data["dtrue"]
rec    = data["rec"];

M = getRegularMesh(domain,m)

resFile = joinpath(dirname(pathof(ConvDiffMIPDECO)),"..","examples",
                    "$(dataset)-$(reg)-noise-$(noiseLevel)-$(m[1])x$(m[2])-penalty.mat")
res = matread(resFile)
dobs = res["dobs"]

srcRelaxed   = res["SourcesRelaxed"]
alphaRelaxed = res["alphaRelaxed"]
Ms           = res["Ms"]

println("solve using alpha = $(alphaRelaxed)")
noiseLevel = norm(dobs-dtrue)/norm(dtrue)
println("estimated noise level: $noiseLevel")
if noiseLevel > 0.5
    error("noise level too large, double check that data was loaded correctly. ")
end

# build inverse problem
x1,x2 = getNodalAxes(M)
P = interpmat(x1,x2, rec);
pFor = getConvDiffFEMParam(M,v=v,sig=sig,P=P,Ainv=getMUMPSsolver());

## configure misfit
Wt         = ones(size(dobs))
sigback    = 0.0
pMis       = getMisfitParam(pFor,Wt,dobs,SSDFun)

## Configure regularization
mref       = zeros(M.nc)
reg        = (m,mr,M,I=1.0) -> wTVReg(m,mr,M,eps=1e-8)

## Configure optimization
maxIter    = 50
minUpdate  = 1e-3
HesPrec    = getSSORRegularizationPreconditioner(1.0,1e-15,50)
cgit       = 5
pcgTol     = 1e-1
modFun     = identityMod
boundsLow  = 0*ones(M.nc)
boundsHigh = 1*ones(M.nc)
maxStep	   = 0.1*maximum(boundsHigh)

## store the configuration
pInv       = getInverseParam(M,modFun,reg,alphaRelaxed,mref,
                             boundsLow,boundsHigh,maxStep=maxStep,
                            pcgMaxIter=cgit,pcgTol=pcgTol,minUpdate=minUpdate,maxIter=maxIter,
                            HesPrec=HesPrec);

neighborhoods = (s->trues(size(s)), s->dilation(s,pInv.MInv,1), s->dilation(s,pInv.MInv,1,0.1,0.9))
neighborhoodsNames = ["all","dilation1", "dilation2"];

nrow = size(Ms,2)
ncol = length(neighborhoods)

results = zeros(m[1],m[2],nrow,length(neighborhoods))
init = zeros(m[1],m[2],nrow,length(neighborhoods))
His = zeros(pInv.maxIter,3,nrow,length(neighborhoods))
times = zeros(nrow,length(neighborhoods))

println("solve using alpha = $(pInv.alpha)")

for k1=1:nrow
    for k2=1:length(neighborhoods)
        println("\n --- mipdecoHeuristic starting with Ms[:,$k1] using $(neighborhoodsNames[k2])--")
        (t0,nsolve0) = (pFor.Ainv.solveTime, pFor.Ainv.nSolve)

        src0 = naiveRounding(Ms[:,k1])
        init[:,:,k1,k2] = src0

        times[k1,k2] = @elapsed begin
            mcTR,DcTR,flagTR,his = mipdecoHeuristic(src0,pInv,pMis,getNeighborhood=neighborhoods[k2],out=0)
            results[:,:,k1,k2] = mcTR
            His[:,:,k1,k2] = his;
		end
        nIter = findlast(His[:,1,k1,k2].>0)

        # p1 = viewImage2D(Ms[:,k1],M)
        # title!("relaxed solution")
        # plot!(rec[:,1],rec[:,2],seriestype=:scatter,legend=false)
        # p2 = viewImage2D(1.0.*src0,M)
        # title!("starting guess")
        # p3 = viewImage2D(1.0.*mcTR,M)
        # title!("MIPDECO result")
        # p4 = viewImage2D(1.0.*(mcTR-src0),M)
        # title!("update")
        # pt = plot(p1,p2,p3,p4)
        # display(pt)

        println("\t\t initial obj. :\t\t$(His[1,1,k1,k2]) ")
        println("\t\t final obj. :\t\t$(His[nIter,1,k1,k2]) ")
        println("\t\t relative :\t\t\t$(His[nIter,1,k1,k2]/His[1,1,k1,k2]) ")
        println("\t\t no. PDE solves:\t$(2*nIter) ")
        println("\t\t runtime: \t\t\t$(times[k1,k2])")
    end
end

res["SourcesMIPDECOPenalty"]=results
res["InitMIPDECOPenalty"]=init
res["HisMIPDECOPenalty"] = His
res["timesMIPDECOPenalty"] = times
matwrite(resFile,res)

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO