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

# read simulation input files

function readAmbsFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["ambulances"]
	n = size(table.data,1) # number of ambulances
	@assert(n >= 1)

	# create ambulances from data in table
	ambulances = Vector{Ambulance}(undef, n)
	columns = table.columns # shorthand
	attributes = parseAttributesColumn(table)
	for i = 1:n
		ambulances[i] = Ambulance()
		ambulances[i].index = columns["index"][i]
		ambulances[i].stationIndex = columns["stationIndex"][i]
		ambulances[i].class = AmbClass(columns["class"][i])
		ambulances[i].attributes = attributes[i]
		@assert(ambulances[i].index == i)
	end

	return ambulances
end

function readArcsFile(filename::String;
	keepAllFields::Bool = false)
	tables = readTablesFromFile(filename)

	# get arcForm, numModes
	table = tables["miscData"]
	arcForm = table.columns["arcForm"][1] # "directed" or "undirected" arcs
	@assert(in(arcForm, ["directed", "undirected"]))
	numModes = table.columns["numModes"][1]
	@assert(numModes >= 1)

	# swap to arcs table
	table = tables["arcs"]

	# initialise arcs, travelTimes
	numArcs = size(table.data,1) # number of arcs, this may double if arcForm == "undirected"
	@assert(numArcs >= 1)
	arcs = []
	travelTimes = []
	if arcForm == "directed"
		arcs = Vector{Arc}(undef, numArcs)
		travelTimes = Array{Float,2}(undef, numModes, numArcs)
	elseif arcForm == "undirected"
		arcs = Vector{Arc}(undef, numArcs * 2)
		travelTimes = Array{Float,2}(undef, numModes, numArcs * 2)
	end

	# check that table header has all travel times and modes
	for j = 1:numModes
		@assert(in(string("mode_", j), table.header), "Missing travel mode $j")
	end

	# fill arc fields with data in table
	c = table.columns # shorthand
	(indexCol = c["index"]); (fromNodeCol = c["fromNode"]); (toNodeCol = c["toNode"]) # shorthand, to avoid repeated dict lookups
	(hasDistCol, distCol) = haskey(c, "distance") ? (true, c["distance"]) : (false, fill(NaN, numArcs)) # shorthand
	modeCols = [c[string("mode_", i)] for i = 1:numModes] # shorthand
	fieldNames = setdiff(table.header, keepAllFields ? [] : ["index", "fromNode", "toNode", "distance", [string("mode_", i) for i = 1:numModes]...])
	rowsFields = tableRowsFieldDicts(table, fieldNames)
	# get arc data from table
	for i = 1:numArcs
		arcs[i] = Arc()
		arcs[i].index = indexCol[i]
		arcs[i].fromNodeIndex = fromNodeCol[i]
		arcs[i].toNodeIndex = toNodeCol[i]
		arcs[i].distance = distCol[i]
		arcs[i].fields = rowsFields[i]
		# read travel times
		for j = 1:numModes
			travelTimes[j,i] = modeCols[j][i]
		end
		@assert(arcs[i].index == i)
		if hasDistCol @assert(arcs[i].distance >= 0) end
	end

	# if arcs are undirected, make them directed
	if arcForm == "undirected"
		# copy arcs, swap directions
		for i = 1:numArcs
			j = i + numArcs
			arcs[j] = Arc()
			arcs[j].index = j
			arcs[j].fromNodeIndex = arcs[i].toNodeIndex
			arcs[j].toNodeIndex = arcs[i].fromNodeIndex
			arcs[j].distance = arcs[i].distance
			arcs[j].fields = deepcopy(arcs[i].fields)
		end
		# copy travel times
		travelTimes[:, (1:numArcs) .+ numArcs] = travelTimes[:, 1:numArcs]
	end

	# travel times should be positive
	@assert(all(travelTimes .> 0))

	return arcs, travelTimes
end

function readCallsFile(filename::String)
	tables = readTablesFromFile(filename)

	# calls file has sim start time
	startTime = tables["miscData"].columns["startTime"][1]
	@assert(startTime > nullTime)

	# create calls from data in calls table
	table = tables["calls"]
	data = table.data # shorthand
	n = size(data,1) # number of calls
	@assert(n >= 1)
	c = table.columns # shorthand
	(indexCol = c["index"]); (priorityCol = c["priority"]); (xCol = c["x"]); (yCol = c["y"]); (arrivalTimeCol = c["arrivalTime"]); (dispatchDelayCol = c["dispatchDelay"]); (onSceneDurationCol = c["onSceneDuration"]); (hospitalIndexCol = c["hospitalIndex"]) # shorthand, to avoid repeated dict lookups
	attributes = parseAttributesColumn(table)
	if haskey(c, "transport") && haskey(c, "handoverDuration")
		(transportCol = c["transport"]); (handoverDurationCol = c["handoverDuration"]);
	else # compat
		@assert(haskey(c, "transfer") && haskey(c, "transferDuration"))
		@warn("`transfer` and `transferDuration` headers in calls file are deprecated, use `transport` and `handoverDuration` instead.")
		(transportCol = c["transfer"]); (handoverDurationCol = c["transferDuration"]);
	end

	calls = Vector{Call}(undef, n)
	for i = 1:n
		calls[i] = Call()
		calls[i].index = indexCol[i]
		calls[i].priority = Priority(Int(priorityCol[i]))
		calls[i].location.x = xCol[i]
		calls[i].location.y = yCol[i]
		calls[i].arrivalTime = arrivalTimeCol[i]
		calls[i].dispatchDelay = dispatchDelayCol[i]
		calls[i].onSceneDuration = onSceneDurationCol[i]
		calls[i].transport = transportCol[i]
		calls[i].handoverDuration = calls[i].transport ? handoverDurationCol[i] : 0.0
		calls[i].hospitalIndex = hospitalIndexCol[i]
		calls[i].attributes = attributes[i]

		@assert(calls[i].index == i)
		@assert(calls[i].priority != nullPriority)
		@assert(calls[i].arrivalTime >= 0)
		@assert(calls[i].dispatchDelay >= 0)
		@assert(calls[i].onSceneDuration >= 0)
		@assert(calls[i].handoverDuration >= 0)
	end

	# check that calls are ordered by arrival time
	for i = 1:n-1
		@assert(calls[i].arrivalTime <= calls[i+1].arrivalTime)
	end

	# check that first call is not before start time
	@assert(startTime <= calls[1].arrivalTime)

	# check that either: no calls have same arrival time, or dispatch delay values are all > 0
	# This is to avoid potential bug where two calls arrive at same time with dispatch delay 0,
	# leading to calls being responded to in order of occurence in calls file, rather than by priority.
	@assert(all(calls[i].arrivalTime < calls[i+1].arrivalTime for i = 1:n-1)
		|| all(calls[i].dispatchDelay > 0 for i = 1:n))

	return calls, startTime
end

function readCompTableFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["compTable"]
	numAmbs = size(table.data,1)
	numStations = size(table.data,2) - 1
	columns = table.columns # shorthand

	# check that numAmbs column in table contains all ambulance counts
	@assert(all(sort(columns["numAmbs"]) .== [1:numAmbs;]))

	# check that table header contains all stations
	for j = 1:numStations
		@assert(in(string("station_", j), table.header))
	end

	# create compliance table
	compTable = CompTable(undef,numAmbs, numStations)
	for i in columns["numAmbs"]
		for j = 1:numStations
			compTable[i,j] = columns[string("station_", j)][i]
		end
	end

	@assert(all(compTable .>= 0))

	# check that compliance table rows add up to row indices
	@assert(vec(sum(compTable, dims=2)) == [1:numAmbs;])

	return compTable
end

function readDemandFile(filename::String)
	tables = readTablesFromFile(filename)

	# create demand object, will populate gradually
	demand = Demand()

	# demand rasters
	table = tables["demandRasters"]
	n = numRasters = demand.numRasters = size(table.data,1)
	@assert(n >= 1)
	columns = table.columns # shorthand
	demand.rasters = Vector{Raster}(undef, n)
	demand.rasterFilenames = Vector{String}(undef, n)
	allunique(columns["rasterFilename"]) || @warn("There are duplicate raster filenames in the demand file, this is unnecessary.")
	demandFileDir = dirname(realpath(filename))
	for i = 1:n
		@assert(columns["rasterIndex"][i] == i)

		rasterFilename = String(columns["rasterFilename"][i])
		rasterFilename = joinPathIfNotAbs(demandFileDir, interpolateString(rasterFilename)) # raster filename should be absolute (after interpolating), or relative to the demand file directory
		@assert(isfile(rasterFilename))
		demand.rasterFilenames[i] = rasterFilename
		demand.rasters[i] = readRasterFile(rasterFilename)
	end

	# demand modes
	table = tables["demandModes"]
	n = numModes = demand.numModes = size(table.data,1)
	@assert(n >= 1)
	columns = table.columns # shorthand
	demand.modes = Vector{DemandMode}(undef, n)
	for i = 1:n
		@assert(columns["modeIndex"][i] == i)

		demand.modes[i] = DemandMode()
		demand.modes[i].index = i
		rasterIndex = demand.modes[i].rasterIndex = columns["rasterIndex"][i]
		@assert(1 <= rasterIndex <= numRasters)
		demand.modes[i].raster = demand.rasters[rasterIndex]
		demand.modes[i].priority = columns["priority"][i] |> Meta.parse |> eval
		demand.modes[i].arrivalRate = columns["arrivalRate"][i]
		@assert(demand.modes[i].arrivalRate >= 0)
		demand.modes[i].rasterMultiplier = demand.modes[i].arrivalRate / sum(demand.modes[i].raster.z)
	end

	# demand sets
	table = tables["demandSets"]
	n = size(table.data,1)
	@assert(n >= 1)
	columns = table.columns # shorthand
	numSets = demand.numSets = maximum(columns["setIndex"])
	@assert(sort(unique(columns["setIndex"])) == [1:numSets;]) # all sets from 1:numSets should be used
	demand.modeLookup = fill(nullIndex, numSets, numPriorities)
	for i = 1:n
		setIndex = columns["setIndex"][i]
		@assert(setIndex == i)
		modeIndices = columns["modeIndices"][i] |> Meta.parse |> eval
		for modeIndex in modeIndices
			@assert(1 <= modeIndex <= numModes)
			priorityIndex = Int(demand.modes[modeIndex].priority)
			@assert(demand.modeLookup[setIndex,priorityIndex] == nullIndex) # value should not yet be filled
			demand.modeLookup[setIndex,priorityIndex] = modeIndex
		end
	end
	@assert(all(demand.modeLookup .!= nullIndex)) # check that demand.modeLookup is filled correctly
	# have not checked that all modes from 1:numModes are used

	# demand sets timing
	table = tables["demandSetsTiming"]
	columns = table.columns # shorthand
	demand.setsStartTimes = columns["startTime"]
	@assert(demand.setsStartTimes[1] > nullTime)
	@assert(issorted(demand.setsStartTimes, lt=<=)) # times should be strictly increasing
	demand.setsTimeOrder = columns["setIndex"]
	@assert(sort(unique(demand.setsTimeOrder)) == [1:numSets;]) # each demand set should be used at least once

	# misc
	demand.recentSetsStartTimesIndex = 1

	return demand
end

function readDemandCoverageFile(filename::String)
	tables = readTablesFromFile(filename)

	dc = demandCoverage = DemandCoverage()

	# demand cover times
	table = tables["coverTimes"]
	n = size(table.data,1) # number of priorities
	@assert(n == numPriorities)
	columns = table.columns # shorthand
	for i = 1:n
		priority = eval(Meta.parse(columns["demandPriority"][i]))
		coverTime = dc.coverTimes[priority] = columns["coverTime"][i]
		@assert(coverTime >= 0)
	end
	@assert(all(p -> haskey(dc.coverTimes, p), priorities), "coverTimes not set for all priorities")

	# number of points to represent demand per raster cell
	table = tables["demandRasterCellNumPoints"]
	@assert(size(table.data,1) == 1) # should be single row
	columns = table.columns # shorthand
	rows = dc.rasterCellNumRows = columns["rows"][1]
	cols = dc.rasterCellNumCols = columns["cols"][1]
	@assert(rows >= 1)
	@assert(cols >= 1)

	return dc
end

function readEventsFile(filename::String)
	tables = readTablesFromFile(filename)

	# input files
	table = tables["inputFiles"]
	inputFiles = table.columns["name"]
	fileChecksums = [parse(UInt, checksumString[2:end-1]) for checksumString in table.columns["checksum"]] # remove quote chars around number, convert to UInt

	# event dictionary and filter
	table = tables["eventDict"]
	eventNames = table.columns["name"]
	eventKeys = table.columns["key"]
	n = length(eventKeys) # number of event types
	filter = haskey(table.columns, "filter") ? Bool.(table.columns["filter"]) : fill(true,n)
	@assert(length(eventNames) == n)
	@assert(allunique(eventNames))
	@assert(allunique(eventKeys))
	eventDict = Dict{Int,EventForm}()
	eventFilter = Dict([e => true for e in instances(EventForm)]) # eventFilter[eventForm] = true if events of eventForm are included in file, false otherwise
	for i = 1:n
		eventForm = eval(Meta.parse(eventNames[i]))
		eventDict[eventKeys[i]] = eventForm
		eventFilter[eventForm] = filter[i]
	end

	# create events from data in events table
	table = tables["events"]
	c = table.columns # shorthand
	(eventIndexCol = c["index"]); (parentIndexCol = c["parentIndex"]); (eventKeyCol = c["eventKey"]); (timeCol = c["time"]); (ambIndexCol = c["ambIndex"]); (callIndexCol = c["callIndex"]); (stationIndexCol = c["stationIndex"]) # shorthand, to avoid repeated dict lookups
	data = table.data # shorthand
	fileEnded = (data[end,1] == "end") # true if writing all events to file ended before file was closed
	n = size(data,1) - fileEnded # number of events in events table, can be less than actual number of events simulated
	events = Vector{Event}(undef, n)
	for i = 1:n
		event = events[i] = Event()
		event.index = eventIndexCol[i]
		event.parentIndex = parentIndexCol[i]
		event.form = eventDict[eventKeyCol[i]]
		event.time = timeCol[i]
		event.ambIndex = ambIndexCol[i]
		event.callIndex = callIndexCol[i]
		event.stationIndex = stationIndexCol[i]

		@assert(event.index > event.parentIndex)
		if i > 1 @assert(event.index > events[i-1].index) end
		@assert(event.form != nullEvent)
		@assert(eventFilter[event.form])
		@assert(event.time > nullTime)
		@assert(event.ambIndex != nullIndex || event.callIndex != nullIndex)
	end

	# set eventsChildren; eventsChildren[i] gives events that are children of the ith event that occurred
	maxEventIndex = isempty(events) ? 0 : events[end].index
	eventsChildren = [Vector{Event}() for i = 1:maxEventIndex]
	for event in events
		j = event.parentIndex
		if j != nullIndex
			push!(eventsChildren[j], event)
		end
	end

	return events, eventsChildren, eventFilter, fileEnded, inputFiles, fileChecksums
end

# read geographic data from file and apply f to data
function readGeoFile(f::Function, filename::String)
	return ArchGDAL.read(filename) do dataset
		f(dataset)
	end
end

function readHospitalsFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["hospitals"]
	n = size(table.data,1) # number of hospitals
	@assert(n >= 1)

	# create hospitals from table data
	columns = table.columns # shorthand
	attributes = parseAttributesColumn(table)
	hospitals = Vector{Hospital}(undef, n)
	for i = 1:n
		hospitals[i] = Hospital()
		hospitals[i].index = columns["index"][i]
		hospitals[i].location.x = columns["x"][i]
		hospitals[i].location.y = columns["y"][i]
		hospitals[i].attributes = attributes[i]
		@assert(hospitals[i].index == i)
	end

	return hospitals
end

function readMapFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["map"]
	columns = table.columns # shorthand

	# create map from data in table
	map = Map()
	map.xMin = columns["xMin"][1]
	map.xMax = columns["xMax"][1]
	map.yMin = columns["yMin"][1]
	map.yMax = columns["yMax"][1]
	map.xScale = columns["xScale"][1]
	map.yScale = columns["yScale"][1]

	@assert(map.xMin < map.xMax)
	@assert(map.yMin < map.yMax)
	@assert(0 < map.xScale)
	@assert(0 < map.yScale)

	map.xRange = map.xMax - map.xMin
	map.yRange = map.yMax - map.yMin

	return map
end

function readMobilisationDelayFile(filename::String)::MobilisationDelay
	tables = readTablesFromFile(filename)

	table = tables["mobilisationDelay"]
	use = table.columns["use"][1]
	distr = table.columns["distribution"][1] |> Meta.parse |> eval
	seed = table.columns["seed"][1]
	expectedDuration = table.columns["expectedDuration"][1]

	return MobilisationDelay(use, DistrRng(distr, seed = seed), expectedDuration)
end

function readMultiCompTableFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["miscData"]
	numAmbs = table.columns["numAmbs"][1]
	numStations = table.columns["numStations"][1]

	table = tables["multiCompTable"]
	columns = table.columns # shorthand
	n = size(table.data,1)

	# check that numAmbs column in table contains all ambulance counts at least once
	@assert(all(in(columns["numAmbs"]), 1:numAmbs))

	# check that table header contains all stations
	@assert(all(j -> in("station_$j", table.header), 1:numStations))

	# create multi-compliance table
	mct = [Vector{Vector{Int}}() for i = 1:numAmbs] # mct[i][j][k] is jth option for number of ambs to assign to station k with i free ambs
	@assert(isa(mct, MultiCompTable))
	for (i,numAmbs) in enumerate(columns["numAmbs"])
		state = zeros(Int, numStations)
		for j = 1:numStations
			state[j] = columns["station_$j"][i]
		end
		if !in(state, mct[numAmbs])
			push!(mct[numAmbs], state)
		end
	end

	checkMultiCompTable(mct, numAmbs = numAmbs, numStations = numStations)

	return mct
end

function readNodesFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["nodes"]
	n = size(table.data,1) # number of nodes
	@assert(n >= 2)
	c = table.columns # shorthand
	(indexCol = c["index"]); (xCol = c["x"]); (yCol = c["y"]) # shorthand, to avoid repeated dict lookups
	accessCol = (haskey(c, "offRoadAccess") ? convert(Vector{Bool}, c["offRoadAccess"]) : fill(true, n))
	fieldNames = setdiff(table.header, ["index", "x", "y", "offRoadAccess"])
	rowsFields = tableRowsFieldDicts(table, fieldNames)

	# create nodes from data in table
	nodes = Vector{Node}(undef, n)
	for i = 1:n
		nodes[i] = Node()
		nodes[i].index = indexCol[i]
		nodes[i].location.x = xCol[i]
		nodes[i].location.y = yCol[i]
		nodes[i].offRoadAccess = accessCol[i]
		nodes[i].fields = rowsFields[i]

		@assert(nodes[i].index == i)
	end

	return nodes
end

function readPrioritiesFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["priorities"]
	n = size(table.data,1) # number of priorities
	@assert(n >= 1)
	@assert(n <= numPriorities)

	# read data from table
	columns = table.columns # shorthand
	targetResponseDurationString = haskey(columns, "targetResponseDuration") ? "targetResponseDuration" : "targetResponseTime" # compat for "targetResponseTime"
	targetResponseDurations = Vector{Float}(undef, n)
	responseTravelPriorities = Dict([p => p for p in priorities]) # default is to have response travel priority equal to call priority
	for i = 1:n
		@assert(columns["priority"][i] == i)
		@assert(eval(Meta.parse(columns["name"][i])) == Priority(i))
		targetResponseDurations[i] = columns[targetResponseDurationString][i]

		if haskey(columns, "responseTravelPriority")
			responseTravelPriorities[Priority(i)] = eval(Meta.parse(columns["responseTravelPriority"][i]))
		end
	end

	# target response durations should be positive
	@assert(all(targetResponseDurations .> 0))

	return targetResponseDurations, responseTravelPriorities
end

function readPriorityListFile(filename::String)::PriorityList
	tables = readTablesFromFile(filename)

	table = tables["priorityList"]
	n = size(table.data,1) # number of ambulances
	@assert(n >= 1)

	# create priority list from data in table
	columns = table.columns # shorthand
	itemCol = haskey(columns, "item") ? columns["item"] : columns["numAmbs"] # allow old column name "numAmbs"
	priorityList = PriorityList(undef,n)
	for i = 1:n
		@assert(itemCol[i] == i)
		priorityList[i] = columns["stationIndex"][i]
	end

	# station indices should be positive
	@assert(all(priorityList .>= 1))

	return priorityList
end

function readPriorityListsFile(filename::String)::Vector{PriorityList}
	tables = readTablesFromFile(filename)

	table = tables["priorityLists"]
	(numAmbs, n) = size(table.data)
	@assert(numAmbs >= 1)
	numPriorityLists = n - 1
	@assert(numPriorityLists >= 1) # number of priority lists

	# create priority list from data in table
	columns = table.columns # shorthand
	itemCol = haskey(columns, "item") ? columns["item"] : columns["numAmbs"] # allow old column name "numAmbs"
	priorityLists = PriorityList[]
	for i = 1:numPriorityLists
		priorityList = PriorityList(undef,numAmbs)
		for j = 1:numAmbs
			@assert(itemCol[j] == j)
			priorityList[j] = columns["priorityList_$i stationIndex"][j]
		end

		# station indices should be positive
		@assert(all(priorityList .>= 1))

		push!(priorityLists, priorityList)
	end

	return priorityLists
end

# read raster file using ArchGDAL package, return as custom Raster type
function readRasterFile(rasterFilename::String)
	@assert(isfile(rasterFilename))

	# open raster, get data, close raster
	# only gets first raster band, may change this later
	(geoTransform, z) = readGeoFile(rasterFilename) do dataset
		rasterband = ArchGDAL.getband(dataset, 1)
		return (ArchGDAL.getgeotransform(dataset), ArchGDAL.read(rasterband))
	end

	# shorthand:
	(nx, ny) = size(z) # number of cells in x and y directions
	(x1, dxdi, dxdj, y1, dydi, dydj) = geoTransform # dxdi is change in x per change in index i, for z[i,j]; likewise for dxdj, dydi, dydj
	dx = dxdi # width of cells in x direction
	dy = dydj # height of cells in y direction (may be negative)

	# data checks
	@assert(dxdj == 0 && dydi == 0) # otherwise raster is sloping, so changing x index changes y value, and vice-versa
	@assert(dx > 0)
	@assert(dy != 0)

	# convert data for easier use
	# find x and y vectors to represent raster cell centres, make sure values are increasing
	xMin = x1 + 0.5*dx # we know dx > 0
	# (xMin, dx, z) = (dx > 0) ? (x1 + 0.5*dx, dx, z) : (x1 + (nx-0.5)*dx, -dx, reverse(z, dims = 1))
	(yMin, dy, z) = (dy > 0) ? (y1 + 0.5*dy, dy, z) : (y1 + (ny-0.5)*dy, -dy, reverse(z, dims = 2))
	x = collect(range(xMin, step=dx, length=nx))
	y = collect(range(yMin, step=dy, length=ny))

	return Raster(x, y, z)
end

function readRedispatchFile(filename::String)::Redispatch
	tables = readTablesFromFile(filename)

	redispatch = Redispatch()

	table = tables["miscData"]
	redispatch.allow = Bool(table.columns["allowRedispatch"][1])

	table = tables["redispatchConditions"]
	n = size(table.data,1)
	priorityPairUsed = fill(false, numPriorities, numPriorities)
	for i = 1:n
		p1 = eval(Meta.parse(table.columns["fromCallPriority"][i]))
		p2 = eval(Meta.parse(table.columns["toCallPriority"][i]))
		redispatch.conditions[Int(p1),Int(p2)] = Bool(table.columns["allowRedispatch"][i])

		# check that condition will not overwrite existing condition
		@assert(priorityPairUsed[Int(p1),Int(p2)] == false, "Redispatch file has multiple conditions for call priority pair [$p1, $p2].")
		priorityPairUsed[Int(p1),Int(p2)] = true

		p1 == p2 && @assert(!redispatch.conditions[Int(p1),Int(p2)], "Redispatch file has condition that allows redispatch for calls of same priority ($p1).")
	end

	# should not allow redispatch from priority: p1 -> p1, or p1 -> p2 -> p1, or p1 -> p2 -> p3 -> p1, etc, otherwise may redispatch indefinitely
	@assert(all(x -> x == 0, redispatch.conditions ^ numPriorities), "Redispatch file has conditions with a cycle that may cause unlimited redispatching.")

	return redispatch
end

function readRNetTravelsFile(filename::String)
	rNetTravels = deserializeFile(filename)
	for (i, rNetTravel) in enumerate(rNetTravels)
		@assert(isa(rNetTravel, NetTravel))
		@assert(rNetTravel.isReduced)
		@assert(rNetTravel.modeIndex == i)
	end
	return rNetTravels
end

function readStationsFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["stations"]
	n = size(table.data,1) # number of stations
	@assert(n >= 1)

	# create stations from data in table
	columns = table.columns # shorthand
	attributes = parseAttributesColumn(table)
	stations = Vector{Station}(undef, n)
	for i = 1:n
		stations[i] = Station()
		stations[i].index = columns["index"][i]
		stations[i].location.x = columns["x"][i]
		stations[i].location.y = columns["y"][i]
		stations[i].capacity = columns["capacity"][i]
		stations[i].attributes = attributes[i]

		@assert(stations[i].index == i)
		@assert(stations[i].capacity >= 0)
	end

	return stations
end

function readStatsControlFile(filename::String)
	tables = readTablesFromFile(filename)

	table = tables["params"]
	param(s::String) = table.columns[s][1]
	warmUpDuration = Float(param("warmUpDuration"))
	periodDurationsString = param("periodDurations")
	doCyclePeriodDurations = Bool(param("doCyclePeriodDurations"))

	@assert(warmUpDuration >= 0.0)

	# create period durations iterator
	@assert(isa(periodDurationsString, String) || isa(periodDurationsString, SubString)) # should be string representation of a vector
	periodDurations = periodDurationsString |> Meta.parse |> eval
	@assert(all(x -> x > 0, periodDurations))
	if doCyclePeriodDurations periodDurations = Iterators.cycle(periodDurations) end
	periodDurationsIter = Stateful(periodDurations)

	recordResponseDurationHist = false
	responseDurationHistBinWidth = nullTime
	if haskey(tables, "responseDurationHist_params")
		table = tables["responseDurationHist_params"]
		recordResponseDurationHist = Bool(table.columns["doRecord"][1])
		binWidth = table.columns["binWidth"][1]
		responseDurationHistBinWidth = Float(typeof(binWidth) <: AbstractString ? (binWidth |> Meta.parse |> eval) : binWidth)
	end

	# pack into dict
	d = Dict{String,Any}()
	d["periodDurationsIter"] = periodDurationsIter
	d["warmUpDuration"] = warmUpDuration
	d["recordResponseDurationHist"] = recordResponseDurationHist
	d["responseDurationHistBinWidth"] = responseDurationHistBinWidth

	return d
end

function readStatsDictFile(filename::String)
	table = readTablesFromFile(filename)["statsDict"]
	cols = table.columns # shorthand
	return Dict([string(cols["stat"][i]) => MeanAndHalfWidth(cols["mean"][i], cols["halfWidth"][i]) for i = 1:size(table.data, 1)])
end

function readTravelFile(filename::String)
	tables = readTablesFromFile(filename)

	# create travel object, will populate gradually
	travel = Travel()

	# travel modes
	table = tables["travelModes"]
	n = numModes = travel.numModes = size(table.data,1)
	@assert(n >= 1)
	# create travel modes from data in table (network based data will need to be filled later)
	columns = table.columns # shorthand
	travel.modes = Vector{TravelMode}(undef, n)
	for i = 1:n
		travel.modes[i] = TravelMode()
		travel.modes[i].index = columns["travelModeIndex"][i]
		travel.modes[i].offRoadSpeed = columns["offRoadSpeed"][i]

		@assert(travel.modes[i].index == i)
		@assert(travel.modes[i].offRoadSpeed > 0)
	end

	# travel sets
	table = tables["travelSets"]
	n = size(table.data,1)
	@assert(n >= 1)
	# populate travel object from data in table
	columns = table.columns # shorthand
	numSets = travel.numSets = maximum(columns["travelSetIndex"])
	@assert(sort(unique(columns["travelSetIndex"])) == [1:numSets;]) # all sets from 1:numSets should be used
	@assert(sort(unique(columns["travelModeIndex"])) == [1:numModes;]) # all travel modes should be used
	travel.modeLookup = fill(nullIndex, numSets, numPriorities)
	for i = 1:n
		setIndex = columns["travelSetIndex"][i]
		priority = eval(Meta.parse(columns["priority"][i]))
		modeIndex = columns["travelModeIndex"][i]

		@assert(travel.modeLookup[setIndex, Int(priority)] == nullIndex) # value should not yet be filled
		@assert(1 <= modeIndex <= numModes)
		travel.modeLookup[setIndex, Int(priority)] = modeIndex
	end
	@assert(all(travel.modeLookup .!= nullIndex)) # check that travel.modeLookup is filled correctly

	# travel sets timing
	table = tables["travelSetsTiming"]
	# populate travel object from data in table
	columns = table.columns # shorthand
	travel.setsStartTimes = columns["startTime"]
	@assert(travel.setsStartTimes[1] > nullTime)
	@assert(issorted(travel.setsStartTimes, lt=<=)) # times should be strictly increasing
	travel.setsTimeOrder = columns["travelSetIndex"]
	@assert(sort(unique(travel.setsTimeOrder)) == [1:numSets;]) # each travel set should be used at least once

	# misc
	travel.recentSetsStartTimesIndex = 1

	return travel
end

# read deployments from file
function readDeploymentsFile(filename::String)
	tables = readTablesFromFile(filename)

	# get counts
	table = tables["miscData"]
	numStations = table.columns["numStations"][1]
	numDeployments = table.columns["numDeployments"][1]

	# deployments
	table = tables["deployments"]
	columns = table.columns # shorthand
	# check number of ambulances
	ambIndexCol = columns["ambIndex"]
	@assert(ambIndexCol == collect(1:length(ambIndexCol)))
	# check that table header has all deployments
	for i = 1:numDeployments
		@assert(in("deployment_$i stationIndex", table.header), "Missing deployment $i")
	end
	# create deployments from table
	deployments = Vector{Deployment}(undef, numDeployments)
	for i = 1:numDeployments
		deployments[i] = columns["deployment_$i stationIndex"]

		for value in deployments[i]
			@assert(in(value, 1:numStations))
		end
	end

	return deployments, numStations
end

function readZhangIpParamsFile(filename::String)
	tables = readTablesFromFile(filename)

	zid = ZhangIpData()

	table = tables["miscParams"]
	@assert(size(table.data,1) == 1) # table should have one data row
	for fname in [:travelTimeCost, :onRoadMoveUpDiscountFactor, :regretTravelTimeThreshold]
		# setfield!(zid, fname, table.columns[string(fname)][1])
		setfield!(zid, fname, convert(fieldtype(typeof(zid), fname), table.columns[string(fname)][1]))
	end
	if haskey(table.columns, "expectedHospitalHandoverDuration")
		zid.expectedHospitalHandoverDuration = table.columns["expectedHospitalHandoverDuration"][1]
	else # compat
		@assert(haskey(table.columns, "expectedHospitalTransferDuration"))
		@warn("`expectedHospitalTransferDuration` header in Zhang IP params file is deprecated, use `expectedHospitalHandoverDuration` instead.")
		zid.expectedHospitalHandoverDuration = table.columns["expectedHospitalTransferDuration"][1]
	end

	table = tables["stationCapacities"]
	n = numStations = size(table.data, 1) # number of stations
	for i = 1:n
		@assert(table.columns["stationIndex"][i] == i)
		push!(zid.stationCapacities, table.columns["capacity"][i])
	end

	table = tables["stationMarginalBenefits"]
	n = size(table.data, 1) # number of stations
	@assert(n == numStations)
	for i = 1:n
		@assert(table.columns["stationIndex"][i] == i)
		push!(zid.marginalBenefits, [])
		for j = 1:zid.stationCapacities[i]
			push!(zid.marginalBenefits[i], table.columns["slot_$j"][i])
		end
	end

	return zid
end


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

# UltraNest docstrings are reproduced here under MIT License with the kind
# permission of the original author
# Johannes Buchner <johannes.buchner.acad@gmx.com>.

"""
    struct ReactiveNestedSampling <: AbstractUltraNestAlgorithm

*Experimental feature, not part of stable public API.*

[UltraNest](https://github.com/JohannesBuchner/UltraNest) reactive nested
sampling algorithm with.

Uses the UltraNest Python package, via
[UltraNest.jl](https://github.com/bat/UltraNest.jl) (and PyCall).

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)

!!! note

    This functionality is only available when the
    [UltraNest](https://github.com/bat/UltraNest.jl) package is loaded (e.g. via
    `import UltraNest`).
"""
@with_kw struct ReactiveNestedSampling{TR<:AbstractDensityTransformTarget} <: AbstractSamplingAlgorithm
    trafo::TR = PriorToUniform()

    # "Indicating whether this parameter wraps around (circular parameter)"
    # wrapped_params::Array{Bool}

    "Test transform and likelihood with this number of random points for errors first. Useful to catch bugs."
    num_test_samples::Int = 2

    "If efficiency goes down, dynamically draw more points from the region between ndraw_min and ndraw_max. If set to False, few points are sampled at once."
    draw_multiple::Bool = true

    "Number of logZ estimators and MLFriends region bootstrap rounds."
    num_bootstraps::Int = 30

    "Minimum number of points to simultaneously propose. Increase this if your likelihood makes vectorization very cheap."
    ndraw_min::Int = 128

    "Maximum number of points to simultaneously propose. Increase this if your likelihood makes vectorization very cheap. Memory allocation may be slow for extremely high values."
    ndraw_max::Int = 65536

    "Update region when the volume shrunk by this amount."
    update_interval_volume_fraction::Float64 = 0.8

    "Update stdout status line every log_interval iterations."
    log_interval::Int = -1

    "Show integration progress as a status line."
    show_status::Bool = true

    # "Callback function when region was rebuilt. Allows to show current state of the live points."
    # viz_callback::Function = nop_func

    "Target evidence uncertainty. This is the std between bootstrapped logz integrators."
    dlogz::Float64 = 0.5

    "Target posterior uncertainty. This is the Kullback-Leibler divergence in nat between bootstrapped integrators."
    dKL::Float64 = 0.5

    "Integrate until this fraction of the integral is left in the remainder. Set to a low number (1e-2 … 1e-5) to make sure peaks are discovered. Set to a higher number (0.5) if you know the posterior is simple."
    frac_remain::Float64 = 0.01

    "Terminate when live point likelihoods are all the same, within Lepsilon tolerance. Increase this when your likelihood function is inaccurate, to avoid unnecessary search."
    Lepsilon::Float64 = 0.001

    "Target number of effective posterior samples."
    min_ess::Int = 400

    "maximum number of integration iterations."
    max_iters::Int = -1

    "Stop after this many likelihood evaluations."
    max_ncalls::Int = -1

    "The algorithm tries to assess iteratively where more samples are needed. This number limits the number of improvement loops."
    max_num_improvement_loops::Int = -1

    "Minimum number of live points throughout the run."
    min_num_live_points::Int = 400

    "Require at least this many live points per detected cluster."
    cluster_num_live_points::Int = 40

    "z-score used as a threshold for the insertion order test. Set to infinity to disable."
    insertion_test_window::Float64 = 10.0

    "Number of iterations after which the insertion order test is reset."
    insertion_test_zscore_threshold::Float64 = 2.0
end
export ReactiveNestedSampling

function bat_sample_impl(
    rng::AbstractRNG,
    target::AnyDensityLike,
    algorithm::ReactiveNestedSampling
)
    density_notrafo = convert(AbstractDensity, target)
    shaped_density, trafo = bat_transform(algorithm.trafo, density_notrafo)
    vs = varshape(shaped_density)
    density = unshaped(shaped_density)

    bounds = var_bounds(density)
    if !(all(isequal(0), bounds.vol.lo) && all(isequal(1), bounds.vol.hi))
        throw(ArgumentError("ReactiveNestedSampling only supports (transformed) densities defined on the unit hypercube"))
    end

    function vec_ultranest_logpstr(V_rowwise::AbstractMatrix{<:Real})
        map(logdensityof(density), nestedview(copy(V_rowwise')))
    end

    ndims = totalndof(vs)
    paramnames = all_active_names(varshape(density_notrafo))

    smplr = UltraNest.ultranest.ReactiveNestedSampler(
        paramnames, vec_ultranest_logpstr, vectorized = true,
        num_test_samples = algorithm.num_test_samples,
        draw_multiple = algorithm.draw_multiple,
        num_bootstraps = algorithm.num_bootstraps,
        ndraw_min = algorithm.ndraw_min,
        ndraw_max = algorithm.ndraw_max
    )

    unest_result = smplr.run(
        log_interval = algorithm.log_interval < 0 ? nothing : algorithm.log_interval,
        show_status = algorithm.show_status,
        #viz_callback = algorithm.# viz_callback,
        dlogz = algorithm.dlogz,
        dKL = algorithm.dKL,
        frac_remain = algorithm.frac_remain,
        Lepsilon = algorithm.Lepsilon,
        min_ess = algorithm.min_ess,
        max_iters = algorithm.max_iters < 0 ? nothing : algorithm.max_iters,
        max_ncalls = algorithm.max_ncalls < 0 ? nothing : algorithm.max_ncalls,
        max_num_improvement_loops = algorithm.max_num_improvement_loops,
        min_num_live_points = algorithm.min_num_live_points,
        cluster_num_live_points = algorithm.cluster_num_live_points,
        insertion_test_window = algorithm.insertion_test_window,
        insertion_test_zscore_threshold = algorithm.insertion_test_zscore_threshold
    )

    r = convert(Dict{String, Any}, unest_result)

    unest_wsamples = convert(Dict{String, Any}, r["weighted_samples"])
    v_trafo_us = nestedview(convert(Matrix{Float64}, unest_wsamples["points"]'))
    logvals_trafo = convert(Vector{Float64}, unest_wsamples["logl"])
    weight = convert(Vector{Float64}, unest_wsamples["weights"])
    samples_trafo = DensitySampleVector(vs.(v_trafo_us), logvals_trafo, weight = weight)
    samples_notrafo = inv(trafo).(samples_trafo)

    uwv_trafo_us = nestedview(convert(Matrix{Float64}, r["samples"]'))
    uwlogvals_trafo = map(logdensityof(density), uwv_trafo_us)
    uwsamples_trafo = DensitySampleVector(vs.(uwv_trafo_us), uwlogvals_trafo)
    uwsamples_notrafo = inv(trafo).(uwsamples_trafo)

    logz = convert(BigFloat, r["logz"])::BigFloat
    logzerr = convert(BigFloat, r["logzerr"])::BigFloat
    logintegral = Measurements.measurement(logz, logzerr)

    ess = convert(Float64, r["ess"])

    return (
        result = samples_notrafo, result_trafo = samples_trafo, trafo = trafo,
        uwresult = uwsamples_notrafo, uwresult_trafo = uwsamples_trafo,
        logintegral = logintegral, ess = ess,
        info = r
    )
end

ext: jl


text:
__precompile__()
module TestOriginal
using OpenEphysLoader
using Main.TestUtilities

using Test, Dates

export write_fheader_fun,
    verify_header

### Helper functions ###
function write_bad_fheader(badtype::Symbol, nbytes::Integer = 1024)
    local outfun
    if badtype == :version
        outfun = write_fheader_fun(nbytes, "header_wrongversion.txt")
    elseif badtype == :format
        outfun = write_fheader_fun(nbytes, "header_wrongformat.txt")
    elseif badtype == :noise
        baddata = rand(UInt8, nbytes)
        outfun = io::IOStream -> write(io, baddata)
    else
        error("badtype unrecognized")
    end
    return outfun
end

function write_fheader_fun(
    nbytes::Integer = 1024,
    headerfile::String = "header.txt"
)
    local head
    headerpath = joinpath(dirname(@__FILE__), "data", headerfile)
    @assert isfile(headerpath) "Could not load header file"
    open(headerpath, "r") do readio
        @assert stat(readio).size >= nbytes "Header not long enough"
        head = read(readio, String)
    end
    trunchead = head[1:nbytes]
    return io::IOStream -> write(io, trunchead)
end

function verify_header(header::OriginalHeader)
    @test header.format == "Open Ephys Data Format"
    @test header.version == v"0.4"
    @test header.headerbytes == 1024
    @test header.description == "each record contains one 64-bit timestamp, one 16-bit sample count (N), 1 uint16 recordingNumber, N 16-bit samples, and one 10-byte record marker (0 1 2 3 4 5 6 7 8 255)"
    @test header.created == DateTime("21-Jul-2015 145012", Dates.DateFormat("d-u-y HHMMSS"))
    @test header.channel == "CH30"
    @test header.channeltype == "Continuous"
    @test header.samplerate == 30000
    @test header.blocklength == 1024
    @test header.buffersize == 1024
    @test isapprox(header.bitvolts, 0.195)
end

### Tests ###
@testset "OriginalHeader" begin
    # OriginalHeader constructor
    filecontext(write_fheader_fun()) do io
        header = OriginalHeader(io)
        verify_header(header)
        @test (show(devnull, header); true) # test that it does not error
        @test (OpenEphysLoader.showcompact(devnull, header); true)
    end

    @test (showerror(devnull, CorruptedException("test")); true)

    # truncated header
    filecontext(write_fheader_fun(512)) do io
        @test_throws CorruptedException OriginalHeader(io)
    end

    # Header with bad content
    filecontext(write_bad_fheader(:noise)) do io
        @test_throws CorruptedException OriginalHeader(io)
    end

    filecontext(write_bad_fheader(:version)) do io
        @test_throws CorruptedException OriginalHeader(io)
    end

    filecontext(write_bad_fheader(:format)) do io
        @test_throws CorruptedException OriginalHeader(io)
    end
end

end


text:
attacks(x1, y1, x2, y2) = (x1 == x2 || y1 == y2 || x1 - y1 == x2 - y2 || x1 + y1 == x2 + y2)

function queens(N, xs...)
    lxs = length(xs)
    c = N - lxs
    c == 0 && return reverse(xs)
    for i=1:N
        any(1:lxs) do j
            attacks(xs[j], N-j+1, i, c)
        end && continue
        v = queens(N, xs..., i)
        v !== nothing && return v
    end
    return nothing
end

using WebIO

N, W = 8, 30
allqueens = queens(N)
println(allqueens)

using Blink
w = Window()
body!(w, dom"div"(
    dom"canvas.#ChessBoard[width=$(N*W), height=$(N*W)]"(),
    dom"script"(
        """
            function draw_queen(ctx, sx, sy){
                ctx.fillStyle = "#C0C0C0";
                ctx.fillText("Q", sx+$W/5, sy+4*$W/5);
                ctx.fillStyle = "#000000";
            }

            function draw_board(ctx, queens){
                n = queens.length
                ctx.fillStyle = "#FFFFFF";
                ctx.fillRect(0, 0, n*$W, n*$W);
                ctx.fillStyle = "#000000";
                sx = 0; sy = 0;

                for (i=n-1; i >= 0; i--){
                    for (j=0; j < n; j++){
                        if ((i+j)%2 == 0)
                            ctx.fillRect(sx, sy, $W, $W);
                        if (j == queens[i] - 1)
                            draw_queen(ctx, sx, sy)
                        sx += $W;
                    }
                    sy += $W; sx = 0;
                }
            }

            var c = document.getElementById("ChessBoard");
            var ctx = c.getContext("2d");
            ctx.font = "$(4*div(W,5))px sans-serif";
            draw_board(ctx, $([allqueens...]));
        """
    )
))

max_stars_repo_path: Chapter 13/NQueens/NQueens.jl

max_issues_repo_path: Chapter 13/NQueens/NQueens.jl

max_forks_repo_path: Chapter 13/NQueens/NQueens.jl


text:
function line2segments(points)
    nsegs = 8
    indices = RPR.rpr_int[]
    last_idx = 0
    for seg in 1:nsegs
        append!(indices, last_idx .+ (0:3))
        last_idx += 3
    end
    indices
end

function to_rpr_object(context, matsys, scene, plot::Makie.Lines)
    points = to_value(plot[1])
    if isempty(points)
        return nothing
    end
    segments = TupleView{4, 3}(RPR.rpr_int(0):RPR.rpr_int(length(points)-1))
    indices = reinterpret(RPR.rpr_int, segments)
    radius = [plot.linewidth[]]
    curve = RPR.Curve(context, points, indices, radius, [Vec2f(0.0)], [length(segments)])
    material = RPR.MaterialNode(matsys, RPR.RPR_MATERIAL_NODE_DIFFUSE)
    set!(material, RPR.RPR_MATERIAL_INPUT_COLOR, to_color(plot.color[]))
    set!(curve, material)
    return curve
end

function to_rpr_object(context, matsys, scene, plot::Makie.LineSegments)
    points = to_value(plot[1])
    segments = TupleView{2, 2}(RPR.rpr_int(0):RPR.rpr_int(length(points)-1))
    indices = RPR.rpr_int[]

    for (a, b) in segments
        push!(indices, a, a, b, b)
    end

    radius = if plot.linewidth[] isa AbstractVector
        @show (length(points)÷2) == plot.linewidth[]
        [plot.linewidth[][1]/1000]
    else
        [plot.linewidth[]/1000]
    end
    curve = RPR.Curve(context, points, indices, radius, [Vec2f(0.0)], [length(points)÷2])
    material = RPR.MaterialNode(matsys, RPR.RPR_MATERIAL_NODE_DIFFUSE)
    color = plot.color[] isa AbstractVector ? plot.color[][1] : plot.color[]
    set!(material, RPR.RPR_MATERIAL_INPUT_COLOR, to_color(color))
    set!(curve, material)
    return curve
end


text:
using ..TicToc

"""
    writepvtu(pvtuprefix, vtkprefixes, fieldnames, FT)

Write a pvtu file with the prefix 'pvtuprefix' for the collection of vtk files
given by 'vtkprefixes' using names  of fields 'fieldnames'. The data in the
`vtu` files is of type `FT`.
"""
function writepvtu(pvtuprefix, vtkprefixes, fieldnames, FT)
    open(pvtuprefix * ".pvtu", "w") do pvtufile
        write(
            pvtufile,
            """
            <?xml version="1.0"?>
            <VTKFile type="PUnstructuredGrid" version="0.1" compressor="vtkZLibDataCompressor" byte_order="LittleEndian">
              <PUnstructuredGrid GhostLevel="0">
                <PPoints>
                  <PDataArray type="$FT" Name="Position" NumberOfComponents="3" format="binary"/>
                </PPoints>
                <PPointData>
            """,
        )

        for name in fieldnames
            write(
                pvtufile,
                """
                      <PDataArray type="$FT" Name="$name" format="binary"/>
                """,
            )
        end

        write(
            pvtufile,
            """
                </PPointData>
            """,
        )

        for name in vtkprefixes
            write(
                pvtufile,
                """
                    <Piece Source="$name.vtu"/>
                """,
            )
        end

        write(
            pvtufile,
            """
              </PUnstructuredGrid>
            </VTKFile>
            """,
        )
    end
    return nothing
end


text:
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--part2"
    help = "Do part 2"
    action = :store_true
    "--test"
    help = "Run test data"
    action = :store_true
end
parsed_args = parse_args(ARGS, s)

mutable struct SNumber
    value :: Union{UInt32, Tuple{SNumber, SNumber}}
end

SNumber(n1::SNumber, n2::SNumber) = SNumber((n1,n2))
SNumber(n1::UInt32, n2::UInt32) = SNumber((SNumber(n1), SNumber(n2)))

function iscompound(n::SNumber)
    isa(n.value, Tuple)
end

function Base.show(io::IO, n::SNumber)
    if iscompound(n)
        print(io, "[")
        Base.show(io, n.value[1])
        print(io, ',')
        Base.show(io, n.value[2])
        print(io, ']')
    else
        print(io, n.value::UInt32)
    end
end

function magnitude(n::SNumber)
    if iscompound(n)
        magnitude(n.value[1]) * 3 + magnitude(n.value[2]) * 2
    else
        n.value
    end
end

function parseNumber(s)
    function parseSide(s)
        c = popfirst!(s)
        if c == '['
            v = parseNumber(s)
        else
            v = SNumber(parse(UInt32, c))
        end
        v
    end
    left = parseSide(s)
    @assert popfirst!(s) == ','
    right = parseSide(s)
    @assert popfirst!(s) == ']'
    SNumber((left, right))
end

function parseLine(l)
    s = Iterators.Stateful(l)
    if popfirst!(s) != '['
        throw(ArgumentError("number starts wit '['"))
    end
    value = parseNumber(s)
    @assert isempty(s)
    value
end

test_suffix = parsed_args["test"] ? "_test" : ""
input = [parseLine(l) for l in readlines("input/day_18$test_suffix.txt")]

function prop_right(n::SNumber, val::UInt32)
    if val == 0xffffffff
        return n
    elseif !iscompound(n)
        return SNumber(n.value + val)
    else
        return SNumber(prop_right(n.value[1], val), n.value[2])
    end
end

function prop_left(n::SNumber, val::UInt32)
    if val == 0xffffffff
        return n
    elseif !iscompound(n)
        return SNumber(n.value+val)
    else
        return SNumber(n.value[1],prop_left(n.value[2], val))
    end
end

function explode(n::SNumber) :: Tuple{Bool, SNumber}
    boom = false
    function do_explode(n, dep) :: Tuple{Tuple{UInt32, UInt32}, SNumber}
        if !iscompound(n)
            return ((0xffffffff,0xffffffff), n)
        elseif dep <= 3
            lex, ln = do_explode(n.value[1], dep+1)
            if boom  # something down the left branch exploded
                return ((lex[1], 0xffffffff), SNumber(ln, prop_right(n.value[2], lex[2])))
            else
                rex, rn = do_explode(n.value[2],dep+1)
                return ((0xffffffff, rex[2]), SNumber(prop_left(ln, rex[1]), rn))
           end
        else
            boom = true
            return ((n.value[1].value, n.value[2].value), SNumber(zero(UInt32)))
        end
    end
    ex, res = do_explode(n, 0)
    (boom, res)
end

function split(n::SNumber) :: Tuple{Bool, SNumber}
    splat = false
    function do_split(n::SNumber)
        if iscompound(n)
            ln = do_split(n.value[1])
            if splat
                return SNumber(ln, n.value[2])
            else
                return SNumber(ln, do_split(n.value[2]))
            end
        elseif n.value >= 10
            splat = true
            return SNumber(div(n.value,UInt32(2)), div(n.value+UInt32(1), UInt32(2)))
        else
            return n
        end
    end
    val = do_split(n)
    splat, val
end

function addsnap(n1::SNumber, n2::SNumber) ::SNumber
    result = SNumber((n1, n2))
    while true
        (boom, result) = explode(result)
        if boom
            continue
        end
        (splat, result) = split(result)
        if !splat
            break
        end
    end
    result
end

function part1()
    sum = foldl(addsnap, input)
    magnitude(sum)
end

function part2()
    maxmag = 0
    l = length(input)
    for i in 1:l
        for o in 1:l
            if i == o
                continue
            end
            m = magnitude(addsnap(input[i],input[o]))
            maxmag = max(m, maxmag)
        end
    end
    maxmag
end

p = parsed_args["part2"] ? part2 : part1
answer = p()
println("The answer is: $answer")


text:
include("Tests/Tests.jl")
include("utilities.jl")

using Test, JuMP

# Needs https://github.com/JuliaOpt/MathOptInterface.jl/pull/669
import Pkg
if Pkg.installed()["MathOptInterface"] > v"0.8.2"
    @testset "Term" begin
        config = MOI.Test.TestConfig()
        optimize!(mock) = MOIU.mock_optimize!(mock, [0.0, 0.0],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0]])
        mock = bridged_mock(optimize!)
        Tests.sos_term_test(mock, config)
        Tests.sdsos_term_test(mock, config)
        Tests.dsos_term_test(mock, config)
    end
end
@testset "Bivariate quadratic" begin
    config = MOI.Test.TestConfig()
    optimize!(mock) = MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0, 1.0],
        (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[1.0, 1.0, -√2]],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0, -1.0, 1.0]])
    mock = bridged_mock(optimize!)
    Tests.sos_bivariate_quadratic_test(mock, config)
    Tests.sdsos_bivariate_quadratic_test(mock, config)
    optimize!(mock) = MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0, 1.0, 1.0],
        (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0.0, 2.0, 2.0, 0.0],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[2.0, -1.0, 0.0]])
    Tests.dsos_bivariate_quadratic_test(mock, config)
end
@testset "Concave then convex cubic" begin
    config = MOI.Test.TestConfig()
    optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 6.0; zeros(12); 6.0; zeros(10)])
    mock = bridged_mock(optimize!)
    Tests.dsos_concave_then_convex_cubic_test(mock, config)
    optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 3.0; zeros(3); 3.0; zeros(7); 3.0; zeros(3); 3.0; zeros(5)])
    mock = bridged_mock(optimize!)
    Tests.sdsos_concave_then_convex_cubic_test(mock, config)
    optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 6.0; zeros(8); 6.0; zeros(6)])
    mock = bridged_mock(optimize!)
    Tests.sos_concave_then_convex_cubic_test(mock, config)
end
@testset "Horn" begin
    config = MOI.Test.TestConfig()
    optimize!_inf(mock) = MOIU.mock_optimize!(mock, MOI.INFEASIBLE)
    optimize!(mock) = MOIU.mock_optimize!(mock,
        [1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0,
         -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0,
         1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0,
         1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0,
         1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0,
         1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 0.0, 0.0, 0.0,
         0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 4.0, 0.0, 0.0, 0.0, 0.0,
         0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0,
         0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    mock = bridged_mock(optimize!_inf, optimize!_inf, optimize!)
    Tests.sos_horn_test(mock, config)
    mock = bridged_mock(optimize!_inf)
    Tests.dsos_horn_test(mock, config)
    mock = bridged_mock(optimize!_inf)
    Tests.sdsos_horn_test(mock, config)
end


text:
#=Descrição

Imprimir na saída padrão uma mesa de xadrez de tamanho N,M onde N é o número de colunas e M é o número de linhas.
Você tem que imprimir 0 quando na mesa for branco e 1 quando  não for.
Índice(0,0) é branco e está na céluado do topo no lado esquesdo na saída produzida.

Entrada de Exemplo:
7,4

Saída de Exemplo:
0101010
1010101
0101010
1010101
=#

include("..\\src\\xadrez.jl")

module XadrezTeste
  include("error_handler_pt_br.jl")
end

using Base.Test

import Xadrez

Test.with_handler(custom_handler) do
    @test xadrez(0,0) == ""
    @test xadrez(0,1) == ""
    @test xadrez(1,0) == ""
    @test xadrez(1,1) == "0\n"
    @test xadrez(2,1) == "01\n"
    @test xadrez(1,2) == "0\n1\n"
    @test xadrez(7,4) == "0101010\n1010101\n0101010\n1010101\n"
    @test xadrez(4,7) == "0101\n1010\n0101\n1010\n0101\n1010\n0101\n"
end

max_stars_repo_name: raphaelpanta/julia-lang-exemplos

max_issues_repo_name: raphaelpanta/julia-lang-exemplos

max_forks_repo_name: raphaelpanta/julia-lang-exemplos


text:
include(joinpath(@__DIR__, "intensitymap.jl"))
include(joinpath(@__DIR__, "data.jl"))


text:
module kara_xml

using Test
include("../src/JuliaKara_noGUI.jl"); using .JuliaKara_noGUI

@testset "JuliaKara XML" begin
    path = joinpath(@__DIR__,"..","test","example.world")
    load_world = kara_xml.JuliaKara_noGUI.xml_load_world(path)
    act = load_world.actors

    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(7, 8))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:tree]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(7, 7))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:tree]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(7, 6))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:tree]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(3, 5))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:mushroom]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(3, 4))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:mushroom]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(3, 3))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:mushroom]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(7, 3))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:leaf]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(8, 3))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:leaf]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(9, 3))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:leaf]
    @test JuliaKara_noGUI.get_actors_at_location(load_world,JuliaKara_noGUI.Location(3, 8))[1].actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:kara]

    path_save = joinpath(@__DIR__,"..","test","example_save.world")
    kara_xml.JuliaKara_noGUI.xml_save_world(load_world,path_save)
    @test isfile(path_save) == true
    rm(path_save)
    @test isfile(path_save) == false
end

end


text:
using FileIO

"""
    save(file_name, config)

Save benchmark data
# Examples
```julia
AcuteBenchmark.save("benchmarkdata.jld2", config)
```
"""
function save(file_name::String, config::StructArray{Funb,T1,T2,T3}) where {T1,T2,T3}
    file_name_parts = split(file_name, '.')
    if length(file_name_parts) > 1 && file_name_parts[end] == "jld2"
        # ignore
    else
        # to make sure it is saved with correct format
        file_name = file_name * ".jld2"
    end

    # convert to array before saving
    configArray = Array(config)

    # Always save/load with configArray name
    FileIO.save(file_name, Dict("configArray"=>configArray))
end

"""
    load(file_name)

load benchmark data
# Examples
```julia
configs = AcuteBenchmark.load("benchmarkdata.jld2")
```
"""
function load(file_name::String)
    file_name_parts = split(file_name, '.')
    if length(file_name_parts) > 1 && file_name_parts[end] == "jld2"
        # ignore
    else
        # to make sure it is saved with correct format
        file_name = file_name * ".jld2"
    end
    # Always save/load with configArray name
    return FunbArray(FileIO.load(file_name)["configArray"])
end

################################################################


text:
module GraphPCA

# package code goes here
export GLPCA, reconstruct, reconstruction_error

"""
Retrieves the matrix for finding Y in graph PCA
"""
function getG(X::AbstractMatrix, β::Number, L::AbstractMatrix)
  if β < 0 || β > 1
    throw(ArgumentError("β must be between 0 and 1"))
  end
  n = size(L,2)
  XtX = X'X
  L = L
  #eet = ones(n,n)
  @assert n == size(XtX,2)
  #Largest eigenvalue of XtX to normalize
  λ = eigmax(Symmetric(XtX))
  η = eigmax(Symmetric(L))
  G = (1 - β)*(I - XtX/λ) + β*(L/η) #+ eet/n)
  Symmetric(G)
end

type GLPCA
  U::Matrix{Float64}
  Q::Matrix{Float64}
  X::Matrix{Float64}
  L::Matrix{Float64}
  β::Number
  k::Int
end

"""
Fits a gLPCA model
"""
function GLPCA(X::AbstractMatrix, β::Number, L::AbstractMatrix, k::Int)
  g = getG(X, β, L)
  F = eigfact(g, 1:k)
  Q = F[:vectors]
  U = X*Q
  GLPCA(U, Q, X, L, β, k)
end

"""
Reconstructs the matrix from low-rank components U and Q
"""
function reconstruct(g::GLPCA)
  g.U*g.Q'
end

function reconstruction_error(g::GLPCA)
  g.X - reconstruct(g)
end

end # module

max_stars_repo_name: mihirparadkar/GraphPCA

max_issues_repo_name: mihirparadkar/GraphPCA

max_forks_repo_name: mihirparadkar/GraphPCA


text:
module BioServices

include("eutils.jl")
include("umls.jl")
include("gggenome.jl")

end


text:
module ShallowWaters

    export run_model, Parameter, ∂x, ∂y, Ix, Iy, ∇²

    using NetCDF, Parameters, Printf, Dates, Interpolations

    include("default_parameters.jl")
    include("grid.jl")
    include("constants.jl")
    include("forcing.jl")
    include("model_setup.jl")
    include("initial_conditions.jl")
    include("preallocate.jl")

    include("time_integration.jl")
    include("ghost_points.jl")
    include("rhs.jl")
    include("gradients.jl")
    include("interpolations.jl")
    include("PV_advection.jl")
    include("continuity.jl")
    include("bottom_drag.jl")
    include("diffusion.jl")
    include("tracer_advection.jl")

    include("feedback.jl")
    include("output.jl")
    include("run_model.jl")

end


text:
#include("gerg2008.jl")
using NLopt,ForwardDiff,DiffResults, Plots, StatsBase, Roots
#held implementation of thermo solver, after that, a interfase to build custom solvers

abstract type AbstractStabilityTest end
abstract type AbstractThermoSolver end
abstract type AbstractHelmholtzSolver <: AbstractThermoSolver end

#xw = zeros(2)
#xw[1] = 1.0
#@btime pressure(m,(322.0)u"kg/m^3",647.0u"K",[1.0])

struct HELDSolver
eps_lambda::Float64
eps_composition::Float64
eps_b::Float64
eps_x::Float64
eps_v::Float64
eps_mu::Float64
eps_g::Float64
min_volume::Float64
max_volume::Float64
min_x::Array{Float64,1}
max_x::Array{Float64,1}
end

#volume from packing fraction
_vol(b,eta)=b/eta
_randvol(b) = _vol(b,rand())

#random molar fraction,
function randfrac(N::Integer)
    x = rand(N)
    return x ./sum(x)
end

#random molar fraction, but considering the presence or absence of elements in a vector
function randfrac(x0::Array)
    length(x0)>0 && begin
    x = similar(x0)
    zerox = 0.0
    for i = 1:length(x0)
     x[i]=ifelse(x0[i]==zerox,zerox,rand())
    end
    x /= sum(x)
    return x
end
end

function randpurefrac(x0::Array)
    x = zeros(length(x0))
    i = rand(LinearIndices(findall(x->x!=0.0,x0)))
    x[i] = 1.0
    return x
end

function randshiftfrac(x0::Array)
    x = x0 + 0.1*rand(length(x0))
    x /=sum(x)
    return x
end

function _stability_function(model::AbstractHelmholtzModel,P0,T0,x0,v,x)
    res = 0.0
    (A,dAdx)= _gradientx2(model,v,T0,x);
        res =A+P0*v
        for i = 1:length(dAdx)-1
            x0[i] != 0 && (res+=dAdx[i]*(x0[i]-x[i]))
        end
    return res
end

###here starts the optimization routine

_minvx(model::AbstractHelmholtzModel,P0,T0,x0) = vcat(min_volume(model,P0,T0,x0),(zeros(length(x0))))
_maxvx(model::AbstractHelmholtzModel,P0,T0,x0) =vcat(max_volume(model,P0,T0,x0),map(z->ifelse(z==0,0.0,1.0),x0))
_randvx(model::AbstractHelmholtzModel,P0,T0,x0) = vcat(random_volume(model,P0,T0,x0),randfrac(x0))

minigerg = GERG2008([19,1]) #H2S + Ch4
P0 = 4.53e6
T0 = 190
x0 = [0.05,0.95]
v0 = 5.690353439153631e-5
x01 = [1.731e-2,0.98269]
x02 = [0.06618,0.93382]

#gives a function with their grad, automatically, a la
function nlopt_form(f,xx,gg)
        if length(gg) > 0
            df = DiffResults.GradientResult(xx)
            df = ForwardDiff.gradient!(df,f,xx)
            gg .= DiffResults.gradient(df)
            return DiffResults.value(df)
        else
            return f(xx)
        end
end

function sumx(x::Vector, grad::Vector)
    if length(grad)>0
    grad[:] .= ones(size(grad))
    grad[1] = 0.0
    end
    return (sum(x[2:end])-1.0)
end

function st(model::AbstractHelmholtzModel,P0,T0,x0)

    fmin0 = vec -> _stability_function(model,P0,T0,x0,vec[1],vec[2:end])
    fmin1 = (x,g) -> nlopt_form(fmin0,x,g)

    stability_N = 50*(length(x0)+1)

    vx = _randvx(model,P0,T0,x0)
    lower = _minvx(model,P0,T0,x0)
    upper = _maxvx(model,P0,T0,x0)

    opt = NLopt.Opt(:LD_SLSQP, length(x0)+1)
    opt.lower_bounds = lower
    opt.upper_bounds = upper
    opt.xtol_rel = 1e-6
    opt.stopval = -1e-8
    opt.min_objective = fmin1
    NLopt.equality_constraint!(opt,sumx)
    best_f = Inf

    for i = 1:stability_N
        (min_f,min_vx,status) = NLopt.optimize(opt, vx)
        min_f < 0.0 && (return (min_f,min_vx,i))
        vx[:] .= _randvx(model,P0,T0,x0)
        #min_f<best_f && begin vx[1]= min_vx[1];best_f=min_f end
        #println((min_f,min_vx,status),i)
    end

    return (best_f,vcat(0.0,x0),-1)
end

function st_test(model::AbstractHelmholtzModel,P0,T0,x0,n::Int=100)
    x = 0
    for i = 1:n
    a =  st(model,P0,T0,x0)[3]
    if a > 0
    x+=1
    end
end
    return x/n
    end

####################
##PHASE 2::searching candidate phases,
#variant with support for true zero values in the composition
#######################

function x_set(x0)
    nc = length(x0)
    nc2 = count(z->z>0,x0)
    Mmin = zeros(nc,nc2-1)
    Mmax = zeros(nc,nc2-1)
    set = findall(z->z>0,x0)
    for i = set[1:end-1]
        mmin_i = x0[i]/2.0
        mmax_i = mmin_i+0.5
        Mmin[set,i].=(1.0-mmin_i)/(nc2-1)
        Mmax[set,i].=(1.0-mmax_i)/(nc2-1)
        Mmin[i,i]=mmin_i
        Mmax[i,i]=mmax_i
        end
    return (hcat(Mmin,Mmax))
end

function _gibbs1v(model,P0,T0,x0)
    G = v -> core_helmholtz(model,v[1],T0,x0)+P0*v[1]
    dG = (grad,v) -> begin
        ForwardDiff.gradient!(grad,G,v)
    end
    lower = [min_volume(model,P0,T0,x0)]
    upper = [max_volume(model,P0,T0,x0)]

    opt = NLopt.Opt(:GN_CRS2_LM, 1)
    opt.xtol_rel = 5e-3
    opt.lower_bounds = lower
    opt.upper_bounds = upper
    fmin1 = (x,g) -> nlopt_form(v -> core_helmholtz(model,v[1],T0,x0)+P0*v[1], x,g)
    opt.min_objective = fmin1

    opt2 = NLopt.Opt(:LD_SLSQP, 1)
    opt2.xtol_abs = 1e-8
    opt2.xtol_rel =  1e-8
    opt2.lower_bounds = lower
    opt2.upper_bounds = upper
    opt2.min_objective = fmin1

    initial_x = lower + 0.7*(upper-lower)
    (min_f,min_x) = NLopt.optimize(opt, initial_x)
    (min_f,min_x) = NLopt.optimize(opt2, min_x)

    return min_x[1]
end

###THIS NEEDS TO BE OPTIMIZED, change to NLopt soon
function m_set(model::AbstractHelmholtzModel,P0,T0,x0,vol=random_volume(model,P0,T0,x0))

    x_set0 = x_set(x0)
    v0 = zeros(size(x_set0,2))
    vmin0 = min_volume(model,P0,T0,x0)
    vmax0 = max_volume(model,P0,T0,x0)
    #v00 =_gibbs1v(model,P0,T0,x0)
    for i = 1:size(x_set0,2)
        #vmin[i] = Roots.find_zero(z->_dgibbs_held1(model,z,P0,T0,Mmin[:,i]), (min_volume(model,P0,T0,x0), max_volume(model,P0,T0,x0)), Roots.A42()) #minimizing gibbs energy
        #vmax[i] = Roots.find_zero(z->_dgibbs_held1(model,z,P0,T0,Mmax[:,i]), (min_volume(model,P0,T0,x0), max_volume(model,P0,T0,x0)), Roots.A42())
        #v0[i] = Roots.find_zero(z->_pressure(model,z,T0,x_set0[:,i])-P0, (vmin0,vmax0), Roots.A42())
        v0[i] = vol
        #vmin[i] = Roots.find_zero(z->_pressure(model,z,T0,Mmin[:,i])-P0, random_volume(model,P0,T0,x0))      #dont, work
        #vmax[i] = Roots.find_zero(z->_pressure(model,z,T0,Mmax[:,i])-P0, random_volume(model,P0,T0,x0))       #dont, work

    end
    M = vcat(transpose(v0),x_set0)

    return M
end

##here using NLopt: is faster, but doest work, fix later

#base problem

function _gibbs1(model,P0,T0,x0)
    G = v -> core_helmholtz(model,v[1],T0,x0)+P0*v[1]
    dG = (grad,v) -> begin
        ForwardDiff.gradient!(grad,G,v)
    end
    lower = [min_volume(model,P0,T0,x0)]
    upper = [max_volume(model,P0,T0,x0)]

    opt = NLopt.Opt(:GN_CRS2_LM, 1)
    opt.xtol_rel = 1e-2
    opt.lower_bounds = lower
    opt.upper_bounds = upper
    fmin1 = (x,g) -> nlopt_form(v -> core_helmholtz(model,v[1],T0,x0)+P0*v[1], x,g)
    opt.min_objective = fmin1

    opt2 = NLopt.Opt(:LD_SLSQP, 1)
    opt2.xtol_abs = 1e-8
    opt2.xtol_rel =  1e-8
    opt2.lower_bounds = lower
    opt2.upper_bounds = upper
    opt2.min_objective = fmin1

    initial_x = lower + 0.7*(upper-lower)
    (min_f,min_x) = NLopt.optimize(opt, initial_x)
    #(min_f,min_x) = NLopt.optimize(opt2, min_x)

    return (min_f,min_x[1])
end

## volume using gibbs minimization

#base problem
function L_v(model::AbstractHelmholtzModel,T0,P0,x0,v,x,lambda)
    res = core_helmholtz(model,v,T0,x) + P0*v
    for i = 1:length(x0)
    res+= lambda[i]*(x0[i]-x[i])
    end
    return res
end

function outer_problem(model::AbstractHelmholtzModel,P0,T0,x0,m_set,UBD,GP)

    opt = NLopt.Opt(:LD_SLSQP, length(x0)+1)
    fmin0 = vec -> vec[1]
    fmin1 = (x,g) -> nlopt_form(fmin0,x,g)
    opt.max_objective = fmin1
    bounds = ones(length(x0)+1)*1e16
    opt.lower_bounds = -bounds
    opt.upper_bounds = bounds

    constraint_gp = vec-> (vec[1]-GP)
    NLopt.inequality_constraint!(opt,(x,g) -> nlopt_form(constraint_gp,x,g),1e-8)

    for i in 1:size(m_set,2)
        constraint_i0 = vec-> (vec[1] - L_v(minigerg,T0,P0,x0,m_set[1,i],m_set[2:end,i],vec[2:end]))
        NLopt.inequality_constraint!(opt,(x,g) ->nlopt_form(constraint_i0,x,g),1e-8)
    end

    (min_f,min_vx,status) = NLopt.optimize(opt,vcat(UBD,-100*ones(length(x0))))
    return (min_f,min_vx,opt)
end

function inner_problem(model::AbstractHelmholtzModel,P0,T0,x0,UBD,lambda)

    opt = NLopt.Opt(:LD_SLSQP, length(x0)+1)
    fmin0 = vec -> L_v(model,P0,T0,x0,vec[1],vec[2:end],lambda)
    fmin1 = (x,g) -> nlopt_form(fmin0,x,g)
    opt.min_objective = fmin1

    opt.lower_bounds = _minvx(model,P0,T0,x0)
    opt.upper_bounds = _maxvx(model,P0,T0,x0)

    constraint_gp = vec -> L_v(model,P0,T0,x0,vec[1],vec[2:end],lambda) - UBD
    NLopt.inequality_constraint!(opt,(x,g) -> nlopt_form(constraint_gp,x,g),1e-8)
    NLopt.equality_constraint!(opt,sumx)
    #NLopt.equality_constraint!(opt,fmin1)

    (min_f,min_vx,status) = NLopt.optimize(opt,_randvx(model,P0,T0,x0))
    return (min_f,min_vx,status)
end

function held(model::AbstractHelmholtzModel,P0,T0,x0)

#step 1:: stability test
step1 = st(minigerg,P0,T0,x0)
step1[3] == -1 && begin
(_,v)=_gibbs1(model,P0,T0,x0)
return (v,x0,1.0)
end
#step 2, inicialization
mayor_k = 0
M = m_set(model,P0,T0,x0)
(GP,_)=_gibbs1(model,P0,T0,x0)
#inicialize UBD, lambda
UBD = GP
LBD = 0.0
(vv,lambda,_)=outer_problem(model,P0,T0,x0,M,UBD,GP)
UBD = vv
(vv2,vx,_) = inner_problem(model,P0,T0,x0,UBD,lambda)
return(vx,lambda)
end


text:
function plot(tb::TBModel; plothopping::Bool=false)
	figure()
	if plothopping
		x_cncts = Vector{Float64}[]
		y_cncts = similar(x_cncts)
		for site in eachindex(tb.neighbors)
			for neighbor in tb.neighbors[site]
				push!(x_cncts, tb.pos[[site,neighbor],1])
				push!(y_cncts, tb.pos[[site,neighbor],2])
			end
		end
		# faster to collect all the connections in one array than to issue multiple calls to
		# PyPlot for each connection
		plot(hcat(x_cncts...), hcat(y_cncts...), "-", color="gray")
		# TODO: Avoid plotting same connection twice...
	end

	plot(tb.pos[:,1], tb.pos[:,2], ".k")
	axis("equal")
end

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
module Variables

include("./constants.jl")
using JuMP
using .Constants

export init_variables

function init_variables(m)

    #Dimensions of box
    @variable(m, lob_length_of_box[i=1:length(boxes)] == boxes[i][1])
    @variable(m, wob_width_of_box[i=1:length(boxes)] == boxes[i][2])
    @variable(m, hob_height_of_box[i=1:length(boxes)] == boxes[i][3])

    #Front-left coordinate of each box
    @variable(m, flbx_front_left_bot_x[i=1:length(boxes)] >= 0)
    @variable(m, flby_front_left_bot_y[i=1:length(boxes)] >= 0)
    @variable(m, flbz_front_left_bot_z[i=1:length(boxes)] >= 0)

    #Check if length (pi) is parallel to X, Y or Z-axis
    @variable(m, lpx_length_of_box_parallel_to_x[i=1:length(boxes)], Bin)
    @variable(m, lpy_length_of_box_parallel_to_y[i=1:length(boxes)], Bin)
    @variable(m, lpz_length_of_box_parallel_to_z[i=1:length(boxes)], Bin)

    #Check if width (qi) is parallel to X, Y or Z-axis
    @variable(m, wpx_width_of_box_parallel_to_x[i=1:length(boxes)], Bin)
    @variable(m, wpy_width_of_box_parallel_to_y[i=1:length(boxes)], Bin)
    @variable(m, wpz_width_of_box_parallel_to_z[i=1:length(boxes)], Bin)

    #Check if height (ri) is parallel to X, Y or Z-axis
    @variable(m, hpx_height_of_box_parallel_to_x[i=1:length(boxes)], Bin)
    @variable(m, hpy_height_of_box_parallel_to_y[i=1:length(boxes)], Bin)
    @variable(m, hpz_height_of_box_parallel_to_z[i=1:length(boxes)], Bin)

    # box_j_x_axis_cover = lpx_length_of_box_parallel_to_x[j]*lob_length_of_box[j] + wpx_width_of_box_parallel_to_x[j]*wob_width_of_box[j]
    # box_j_y_axis_cover = lpy_length_of_box_parallel_to_y[j]*lob_length_of_box[j] + wpy_width_of_box_parallel_to_y[j]*wob_width_of_box[j]
    @expressions(m, begin
        box_x_axis_cover[i=1:length(boxes)],
        (
        m[:lpx_length_of_box_parallel_to_x][i]*fix_value(m[:lob_length_of_box][i]) +
        m[:wpx_width_of_box_parallel_to_x][i]*fix_value(m[:wob_width_of_box][i])
        )

        box_y_axis_cover[i=1:length(boxes)],
        (
        m[:lpy_length_of_box_parallel_to_y][i]*fix_value(m[:lob_length_of_box][i]) +
        m[:wpy_width_of_box_parallel_to_y][i]*fix_value(m[:wob_width_of_box][i])
        )
    end)

    #Make objective available for the whole program. Later used in objective function.
    @expression(
        m, heights[i=1:length(boxes)], m[:flbz_front_left_bot_z][i] + fix_value(m[:hob_height_of_box][i])
    )
    @variable(m, H >= maxHeight >= 0)
    @constraint(m, m[:maxHeight] .>= m[:heights])

    return m
end

end


text:
# episodes.jl
## https://developer.spotify.com/documentation/web-api/reference/episodes/

## https://developer.spotify.com/documentation/web-api/reference/episodes/get-an-episode/
@doc """
# Get an Episode
**Summary**: Get Spotify catalog information for a single peisode identified by it's unique Spotify ID.\n

`episode_id` _Required_: The Spotify ID for the episode_id. Up to 50 episodes can be passed, comma delimited.\n
`market` _Optional_: An ISO 3166-1 alpha 2 country code.

[Reference](https://developer.spotify.com/documentation/web-api/reference/episodes/get-an-episode/)
""" ->
function episodes_get(episode_id, market="US")
    return spotify_request("episodes/$episode_id?market=$market")
end


text:
"""
    StrideTransform(x)

Subtype of `IndexTransform` that transforms and index using stride layout information
derived from `x`.
"""
struct StrideTransform{N,S,O} <: IndexTransform{N,1}
    strides::S
    offsets::O

    StrideTransform{N}(s::S, o::O) where {N,S<:Tuple,O<:Tuple} = new{N,S,O}(s, o)
    StrideTransform(s::CanonicalInt, o::CanonicalInt) = StrideTransform{1}((s,), (o,))
    function StrideTransform{N}(a::A) where {N,A}
        StrideTransform{N}(ArrayInterface.strides(a), offsets(a))
    end
    StrideTransform(a::A) where {A} = StrideTransform{ndims(A)}(a)
end

#const OffsetTransform{O} = StrideTransform{1,Tuple{Static{1}},Tuple{O}}
#OffsetTransform(o::CanonicalInt) = StrideTransform(static(1), o)

@generated function _strides2int(o::O, s::S, i::I) where {O,S,I}
    N = known_length(S)
    out = :()
    for i in 1:N
        tmp = :(((getfield(i, $i) - getfield(o, $i)) * getfield(s, $i)))
        out = ifelse(i === 1, tmp, :($out + $tmp))
    end
    return Expr(:block, Expr(:meta, :inline), out)
end
function (t::StrideTransform{N})(i::AbstractCartesianIndex{N}) where {N}
    _strides2int(getfield(t, :offsets), getfield(t, :strides), Tuple(i))
end

@inline function (t::StrideTransform{1})(i::CanonicalInt)
    (i - getfield(getfield(t, :offsets), 1)) * getfield(getfield(t, :strides), 1)
end
@inline function (t::StrideTransform{1})(i::AbstractRange{Int})
    t(static_first(i)):t(static_step(i)):t(static_last(i))
end

(t::StrideTransform{1})(i::AbstractVector{Int}) = TransformedIndex(i, t)

function (t::StrideTransform{N})(i::AbstractArray{T}) where {N,T<:AbstractCartesianIndex{N}}
    TransformedIndex(i, t)
end

struct ReduceTransform{I,F} <: CoordinateTransform{I,1}
    f::F
end

(t::ReduceTransform{I,F})(x::AbstractCartesianIndex) where {I,F} = t(Tuple(x))
@inline (t::ReduceTransform{I,F})(x::Tuple) where {I,F} = Static.reduce_tup(getfield(t, :f), x)

function ReduceTransform(f::F, t::CoordinateTransform{<:Any,O}) where {F,O}
    ComposedTransform(ReduceTransform{O,F}(f), t)
end
function Base.show(io::IO, ::MIME"text/plain", @nospecialize(x::StrideTransform))
    print(io, "StrideTransform{$(ndims(x))}($(strides(x)), $(offsets(x)))")
end


text:
module TestTypeDef
using Compat
using Compat.Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType, level,  reftype, leveltype, catvalue, iscatvalue

@testset "CategoricalPool, a b c order" begin
    pool = CategoricalPool(
        [
            "a",
            "b",
            "c"
        ],
        Dict(
            "a" => DefaultRefType(1),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(3),
        )
    )

    @test iscatvalue(Int) == false
    @test iscatvalue(Any) == false
    @test iscatvalue(Missing) == false

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector)
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict)
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)

    # leveltype() only accepts categorical value type
    @test_throws ArgumentError leveltype("abc")
    @test_throws ArgumentError leveltype(String)
    @test_throws ArgumentError leveltype(1.0)
    @test_throws ArgumentError leveltype(Int)

    for i in 1:3
        x = catvalue(i, pool)

        @test iscatvalue(x)
        @test iscatvalue(typeof(x))
        @test eltype(x) === Char
        @test eltype(typeof(x)) === Char
        @test leveltype(x) === String
        @test leveltype(typeof(x)) === String
        @test reftype(x) === DefaultRefType
        @test reftype(typeof(x)) === DefaultRefType
        @test x isa CategoricalArrays.CategoricalString{DefaultRefType}

        @test isa(level(x), DefaultRefType)
        @test level(x) === DefaultRefType(i)

        @test isa(CategoricalArrays.pool(x), CategoricalPool)
        @test CategoricalArrays.pool(x) === pool

        @test typeof(x)(x) === x
    end
end

@testset "CategoricalPool, c b a order" begin
    pool = CategoricalPool(
        [
            "a",
            "b",
            "c"
        ],
        Dict(
            "a" => DefaultRefType(1),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(3),
        ),
        [
            DefaultRefType(3),
            DefaultRefType(2),
            DefaultRefType(1),
        ]
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector)
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict)
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(3)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(1)

    for i in 1:3
        y = catvalue(i, pool)

        @test iscatvalue(y)

        @test isa(level(y), DefaultRefType)
        @test level(y) === DefaultRefType(i)

        @test isa(CategoricalArrays.pool(y), CategoricalPool)
        @test CategoricalArrays.pool(y) === pool

        @test typeof(y)(y) === y
    end
end

end


text:
s = open("input.txt") do file
    readlines(file);
end

startTime = parse(Int, s[1]);

nums = split(s[2], ',');
offsets = findall(x->x!="x", nums) .- 1;
busNums = map(x->parse(Int,x), nums[offsets.+1]);

function FindNextMatchup(startTime, iter, offset, num)
    curTime = startTime;
    while (curTime+offset)%num != 0
        curTime += iter;
    end
    println((curTime, iter*num));
    return (curTime, iter*num);
end

function GetPatternTime(start, startIter, indOrder, offsets, nums)
    iter = startIter;
    time = start;
    for i = 1:length(indOrder)
        curInd = indOrder[i];
        time, iter = FindNextMatchup(time, iter, offsets[curInd], nums[curInd]);
    end

    return (time, iter);
end

function CheckTime(time, offsets, nums)
    println(time);
    times = offsets .+ time;
    mods = times.%nums;
    println(mods);

    return all(mods.==0);
end

startTime = 100000000000000;
while startTime % busNums[1] != 0
    startTime += 1
end

# The start time came from past brute force efforts,
# But you could just use the above
sortedIndices = sortperm(busNums[2:end]);
theTime, bigIter = GetPatternTime(141883788274894, busNums[1], sortedIndices, offsets[2:end], busNums[2:end]);
CheckTime(theTime, offsets, busNums)


text:
@testset "Lists" begin

@test filter(iseven, List(1, 2, 3)) == List(2)
@test filter(!iseven, List(1, 2, 3)) == List(1, 3)
@test List(1, 2, 3) ++ List(4, 5, 6) == List(1, 2, 3, 4, 5, 6)
@test append(List(1), List(2), List(3)) == List(1, 2, 3)

@testset "repr" begin
    @test repr(list(1, 2, 3)) == "(1 2 3)"
    @test repr(list(1, true, nothing)) == "(1 #t #<void>)"
end

@testset "map" begin
    @test map(-, List(1, 2, 3)) == List(-1, -2, -3)
    @test map(+, List(1, 2), List(3, 4)) == List(4, 6)
end

@testset "broadcast" begin
    α = List(1, 2, 3)
    β = List(3, 2, 1)
    @test (x -> x^2).(α) == List(1, 4, 9)
    @test α .+ β == List(4, 4, 4)
    @test 3 .* α == List(3, 6, 9)
    @test α .- (α .- 3 .* β) == List(9, 6, 3)

    # This is somewhat hard to get right, because Lists can’t efficiently act like
    # one-dimensional arrays.
    @test_broken α .+ [2, 2, 2] == [3, 4, 5]
end

end