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

# initialise travel data for fGraph
# arcTravelTimes[i,j] gives travel time along arc j (in fGraph) for travel mode i
function initFNetTravels!(net::Network, arcTravelTimes::Array{Float,2})
	numModes = size(arcTravelTimes, 1) # shorthand
	@assert(length(net.fGraph.arcs) == size(arcTravelTimes, 2)) # check number of arcs

	# create fNetTravels
	net.fNetTravels = [NetTravel(false) for i = 1:numModes]
	for i = 1:numModes
		fNetTravel = net.fNetTravels[i]
		fNetTravel.modeIndex = i
		fNetTravel.arcTimes = arcTravelTimes[i,:]
	end
end

# return true if fNode i is in rGraph, false otherwise
function isFNodeInRGraph(net::Network, i::Int)
	return (net.fNodeRNode[i] != nullIndex)
end

# create the reduced graph from the full graph, for the given network
function createRGraphFromFGraph!(net::Network)
	# assumptions:
	# - no nodes are connected to themselves
	# - fGraph has no duplicate arcs

	net.rGraph = Graph(true)

	# shorthand names:
	fGraph = net.fGraph
	fNodes = fGraph.nodes
	fArcs = fGraph.arcs
	numFNodes = length(fNodes)
	rGraph = net.rGraph

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

	# some useful fGraph connectivity data

	# shorthand:
	fNodesOut = fGraph.light.fadjlist # fNodesOut[i] gives indices of fNodes incident to arcs outgoing from fNode i
	fNodesIn = fGraph.light.badjlist # fNodesIn[i] gives indices of fNodes incident to arcs incoming to fNode i

	fNodeNumAdjacent = [length(union(fNodesIn[i], fNodesOut[i])) for i = 1:numFNodes] # number of fNodes adjacent to each fNode

	# check that no node is connected to itself
	for i = 1:numFNodes
		@assert(!in(i, fNodesIn[i]))
	end

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

	# find intersection/decision nodes in fGraph, add to rGraph

	# only keep nodes with >= 3 adjacent nodes, or 2 adjacent nodes if the incoming and outgoing arcs cannot be replaced
	# Also, for easier programming, put leaf nodes in fGraph into rGraph, hopefully there are not too many...
	keepNode = Vector{Bool}(undef, numFNodes) # find which fNodes belong in rGraph
	keepNode[:] .= false
	for i = 1:numFNodes
		if fNodeNumAdjacent[i] >= 3
			keepNode[i] = true
		elseif fNodeNumAdjacent[i] == 2 && (length(fNodesIn[i]) + length(fNodesOut[i]) == 3)
			keepNode[i] = true
		elseif fNodeNumAdjacent[i] == 1
			keepNode[i] = true
		end
	end
	if !any(keepNode) keepNode[1] = true end # if there are still no rNodes, make one (can happen if network is just a loop)
	rGraph.nodes = deepcopy(fNodes[keepNode]) # need to renumber rNode indices later
	net.rNodeFNode = findall(keepNode)
	net.fNodeRNode = [nullIndex for i = 1:numFNodes]
	net.fNodeRNode[keepNode] = [1:sum(keepNode);]

	# special consideration: need to add one decision node to any loop with currently only one rNode
	# Otherwise, fNodes in the loop still have a decision as to how to reach the rNode
	for ri = 1:length(rGraph.nodes)
		fi = net.rNodeFNode[ri]
		for fj = fNodesOut[fi] # fNodes outgoing from fNode[fi]
			firstFNode = fj
			prevFNode = fi
			fNode = fj
			# traverse fGraph from rNode[ri], until another rNode or dead end (leaf node?) reached
			while !isFNodeInRGraph(net, fNode) # !(isFNodeInRGraph(net, fNode) || fNodeNumAdjacent[fNode] == 1)
				# traverse along an arc to next node
				nextFNode = first(setdiff(fNodesOut[fNode], prevFNode)) # fNodesOut[fNode] = indices of fNodes outgoing from fNode[fNode]
				prevFNode = fNode
				fNode = nextFNode
			end
			# if we have looped back to same rNode, add firstFNode as decision nodes
			if ri == net.fNodeRNode[fNode]
				keepNode[firstFNode] = true
				net.fNodeRNode[firstFNode] = numFNodes + 1 # a temporary value that will make isFNodeInRGraph(net, firstFNode) return true
				@assert(isFNodeInRGraph(net, firstFNode))
			end
		end
	end

	rGraph.nodes = deepcopy(fNodes[keepNode])
	# renumber rNode indices
	for i = 1:length(rGraph.nodes)
		rGraph.nodes[i].index = i
	end
	net.rNodeFNode = findall(keepNode)
	net.fNodeRNode = [nullIndex for i = 1:numFNodes]
	net.fNodeRNode[keepNode] = [1:sum(keepNode);]

	# shorthand names:
	rNodes = rGraph.nodes
	fNodeRNode = net.fNodeRNode
	rNodeFNode = net.rNodeFNode
	numRNodes = length(rNodes)

	@assert(numRNodes > 0) # should not be possible at this point, and would cause problems below if false

	# check index conversion between fNodes and rNodes
	@assert(all([i == fNodeRNode[rNodeFNode[rNodes[i].index]] for i = 1:numRNodes]))
	# @assert(all([i == fNodeRNode[rNodeFNode[i]] for i = 1:numRNodes]))

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

	# calculate rArcs, rArcFNodes, fNodeRArcs
	# find outgoing arcs from each rNode to other rNodes
	net.fNodeFromRNodes = [Vector{Int}() for i = 1:numFNodes]
	net.rArcFNodes = []
	net.fNodeRArcs = [Vector{Int}() for i = 1:numFNodes]
	rArcFNodes = net.rArcFNodes # shorthand
	fNodeRArcs = net.fNodeRArcs # shorthand
	rArcs = rGraph.arcs # shorthand
	numRArcs = 0
	for ri = 1:numRNodes
		fi = rNodeFNode[ri]
		for fj = fNodesOut[fi] # fNodes outgoing from fNode[fi]
			numRArcs += 1 # will create another rArc
			# traverse fGraph from rNode[ri], until another rNode or dead end (leaf node?) reached
			fNodesOnRArc = [fi] # indices of fNodes on rArc
			prevFNode = fi
			fNode = fj
			while !isFNodeInRGraph(net, fNode)
				push!(fNodeRArcs[fNode], numRArcs)
				push!(fNodesOnRArc, fNode)
				# traverse along an arc to next node
				nextFNode = first(setdiff(fNodesOut[fNode], prevFNode)) # fNodesOut[fNode] = indices of fNodes outgoing from fNode
				prevFNode = fNode
				fNode = nextFNode
			end
			push!(fNodesOnRArc, fNode)
			rj = fNodeRNode[fNode] # rNode at end of rArc, fNodeRNode[fNode] = nullIndex if fNode[fNode] is not in rGraph
			@assert(rj != nullIndex) # should not happen; earlier we put all leaf nodes in fGraph in rGraph
			@assert(ri != rj) # ri == rj would mean that arc is a loop, should have added enough rNodes to remove this problem earlier

			# add rArc
			arc = Arc()
			arc.index = numRArcs
			arc.fromNodeIndex = ri
			arc.toNodeIndex = rj
			push!(rArcs, arc)
			push!(rArcFNodes, fNodesOnRArc)
		end
	end

	@assert(numRArcs > 0)

	# calculate rArcFNodeIndex from rArcFNodes
	net.rArcFNodeIndex = [Dict{Int,Int}() for i = 1:numRArcs]
	for ri = 1:numRArcs
		for (i, fi) in enumerate(rArcFNodes[ri])
			net.rArcFNodeIndex[ri][fi] = i
		end
	end

	# set fNodeRArcs for each rNode (for corresponding fNode)
	for i = 1:numRArcs
		fromRNode = rArcs[i].fromNodeIndex
		toRNode = rArcs[i].toNodeIndex
		push!(fNodeRArcs[rNodeFNode[fromRNode]], i)
		push!(fNodeRArcs[rNodeFNode[toRNode]], i)
	end

	# check that fNodeRArcs is not empty for rNodes, and other nodes have 1-2 entries
	for fi = 1:numFNodes
		if isFNodeInRGraph(net, fi)
			@assert(!isempty(fNodeRArcs[fi]))
		else
			@assert(length(fNodeRArcs[fi]) >= 1 && length(fNodeRArcs[fi]) <= 2)
		end
	end

	# check that start and end fNodes on each rArc are also in rGraph, and arcs are numbered 1 to n
	for i = 1:numRArcs
		@assert(isFNodeInRGraph(net, rArcFNodes[i][1]))
		@assert(isFNodeInRGraph(net, rArcFNodes[i][end]))
		@assert(rArcs[i].index == i)
	end

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

	# calculate fNodeToRNodes and fNodeFromRNodes

	net.fNodeToRNodes = [Vector{Int}() for i = 1:numFNodes]
	net.fNodeFromRNodes = [Vector{Int}() for i = 1:numFNodes]
	net.fNodeToRNodeNextFNode = [Dict{Int,Int}() for i = 1:numFNodes]
	@unpack fNodeToRNodes, fNodeFromRNodes, fNodeToRNodeNextFNode = net

	# traverse each rArc
	for rArc = rArcs
		fNodesOnRArc = rArcFNodes[rArc.index]
		@assert(isFNodeInRGraph(net, fNodesOnRArc[1]))
		@assert(isFNodeInRGraph(net, fNodesOnRArc[end]))
		startRNode = fNodeRNode[fNodesOnRArc[1]]
		endRNode = fNodeRNode[fNodesOnRArc[end]]
		# for fNodes on arc (excluding rNodes):
		for i = 2:length(fNodesOnRArc)-1
			fNode = fNodesOnRArc[i]
			push!(fNodeToRNodes[fNode], endRNode)
			push!(fNodeFromRNodes[fNode], startRNode)
			fNodeToRNodeNextFNode[fNode][endRNode] = fNodesOnRArc[i+1]
		end
	end
	# for rNodes, set corresponding fNode to have fNodeToRNodes and fNodeFromRNodes return the rNode
	for ri = 1:numRNodes
		fi = rNodeFNode[ri]
		fNodeFromRNodes[fi] = [ri]
		fNodeToRNodes[fi] = [ri]
	end

	# check that fNodeFromRNodes and fNodeToRNodes have been filled
	for i = 1:numFNodes
		@assert(!isempty(fNodeFromRNodes[i]))
		@assert(!isempty(fNodeToRNodes[i]))
	end

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

	# arc distances

	# set rArcFArcs
	rArcFArcs = net.rArcFArcs = [Vector{Int}() for i = 1:numRArcs]
	for rArc in rArcs
		fNodesOnRArc = rArcFNodes[rArc.index]
		for i = 1:length(fNodesOnRArc)-1
			fNode = fNodesOnRArc[i]
			nextFNode = fNodesOnRArc[i+1]
			push!(rArcFArcs[rArc.index], fGraph.nodePairArcIndex[fNode, nextFNode])
		end
		@assert(length(rArcFArcs[rArc.index]) >= 1) # at least one fArc per rArc
	end
	# @assert(sort(vcat(rArcFArcs...)) == 1:length(fArcs)) # rArcFArcs contains each fArc exactly once

	# set arc distances
	@assert(all(arc -> !isnan(arc.distance) && arc.distance >= 0, fArcs))
	for rArc in rArcs
		rArc.distance = 0.0
		for fArcIndex in rArcFArcs[rArc.index]
			rArc.distance += fArcs[fArcIndex].distance
		end
	end

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

	# set fNodeToRNodeDist, fNodeFromRNodeDist
	fNodeToRNodeDist = net.fNodeToRNodeDist = [Dict{Int, Float}() for i = 1:numFNodes]
	fNodeFromRNodeDist = net.fNodeFromRNodeDist = [Dict{Int, Float}() for i = 1:numFNodes]

	for rArc in rArcs
		# traverse arc, calculate fNodeFromRNodeDist and fNodeToRNodeDist
		startRNode = rArc.fromNodeIndex
		endRNode = rArc.toNodeIndex
		dist = 0.0
		for fArcIndex in rArcFArcs[rArc.index]
			fArc = fArcs[fArcIndex]
			fNodeFromRNodeDist[fArc.fromNodeIndex][startRNode] = dist
			dist += fArc.distance
			fNodeToRNodeDist[fArc.toNodeIndex][endRNode] = rArc.distance - dist
		end
	end

	# for rNodes, check that corresponding fNode has fNodeFromRNodeDist and fNodeToRNodeDist return 0
	for ri = 1:numRNodes
		fi = rNodeFNode[ri]
		@assert(fNodeFromRNodeDist[fi][ri] == 0.0)
		@assert(fNodeToRNodeDist[fi][ri] == 0.0)
	end

	# for fNodes along each rArc, fNodeFromRNodeDist and fNodeToRNodeDist should add to rArc distance
	for rArc in rArcs
		fNodesOnRArc = rArcFNodes[rArc.index]
		startRNode = fNodeRNode[fNodesOnRArc[1]]
		endRNode = fNodeRNode[fNodesOnRArc[end]]
		for i = 2:length(fNodesOnRArc)-1
			fNode = fNodesOnRArc[i]
			@assert(isapprox(fNodeFromRNodeDist[fNode][startRNode] + fNodeToRNodeDist[fNode][endRNode], rArc.distance; rtol = eps(Float)))
		end
	end

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

	# calculate fields for rGraph: arcDists, light
	initGraph!(rGraph)
end

# calculate arcTimes, spTimes, and spSuccs, for rNetTravels
# also, calculate fNodeToRNodeTime, and fNodeFromRNodeTime, for fNetTravels...
# can only call after rGraph is created
function createRNetTravelsFromFNetTravels!(net::Network;
	rNetTravelsLoaded::Vector{NetTravel} = NetTravel[])

	# shorthand names:
	@unpack fGraph, rGraph, fNetTravels, rArcFNodes, fNodeRNode, rNodeFNode, fNodeToRNodes = net
	numTravelModes = length(fNetTravels)
	numFNodes = length(fGraph.nodes)
	numRNodes = length(rGraph.nodes)
	numRArcs = length(rGraph.arcs)

	rNetTravels = net.rNetTravels = [NetTravel(true) for i = 1:numTravelModes]

	@assert(length(fGraph.arcs) == length(fNetTravels[1].arcTimes))

	for travelModeIndex = 1:numTravelModes
		# shorthand:
		fNetTravel = fNetTravels[travelModeIndex]
		rNetTravel = rNetTravels[travelModeIndex]

		rNetTravel.arcTimes = Vector{Float}(undef, numRArcs)
		fNetTravel.fNodeToRNodeTime = [Dict{Int, Float}() for i = 1:numFNodes]
		fNetTravel.fNodeFromRNodeTime = [Dict{Int, Float}() for i = 1:numFNodes]
		fNetTravel.rArcFNodesTimes = [Vector{Float}() for i = 1:numRArcs]
		fNodeToRNodeTime = fNetTravel.fNodeToRNodeTime # shorthand
		fNodeFromRNodeTime = fNetTravel.fNodeFromRNodeTime # shorthand
		rArcFNodesTimes = fNetTravel.rArcFNodesTimes

		rNetTravel.modeIndex = fNetTravel.modeIndex
		rNetTravel.arcDists = rGraph.arcDists
		rNetTravel.fadjList = rGraph.fadjList

		for rArcIndex = 1:numRArcs
			fNodesOnRArc = rArcFNodes[rArcIndex]
			startRNode = fNodeRNode[fNodesOnRArc[1]]
			endRNode = fNodeRNode[fNodesOnRArc[end]]

			# traverse arc, calculate rArcTime and fNodeFromRNodeTime
			rArcTime = 0.0
			for i = 1:length(fNodesOnRArc)-1
				fNode = fNodesOnRArc[i]
				nextFNode = fNodesOnRArc[i+1]
				fNodeFromRNodeTime[fNode][startRNode] = rArcTime
				push!(rArcFNodesTimes[rArcIndex], rArcTime)
				fArcIndex = fGraph.nodePairArcIndex[fNode, nextFNode]
				rArcTime += fNetTravel.arcTimes[fArcIndex]
			end
			push!(rArcFNodesTimes[rArcIndex], rArcTime)
			rNetTravel.arcTimes[rArcIndex] = rArcTime

			# calculate fNodeToRNodeTime
			for i = 2:length(fNodesOnRArc)-1
				fNode = fNodesOnRArc[i]
				fNodeToRNodeTime[fNode][endRNode] = rArcTime - fNodeFromRNodeTime[fNode][startRNode]
			end
			fNodeToRNodeTime[fNodesOnRArc[end]][endRNode] = 0.0 # time to self is 0

			# check length of rArcFNodesTimes is same as length of rArcFNodes for current arc
			@assert(length(fNodesOnRArc) == length(rArcFNodesTimes[rArcIndex]))
		end

		# # for rNodes, set corresponding fNode to have fNodeFromRNodeTime and fNodeToRNodeTime return 0
		# for ri = 1:numRNodes
			# fi = rNodeFNode[ri]
			# # fNodeToRNodeTime[fi][ri] = 0.0 # should already be true?
			# # fNodeFromRNodeTime[fi][ri] = 0.0 # should already be true?
		# end

		# for rNodes, check that corresponding fNode has fNodeFromRNodeTime and fNodeToRNodeTime return 0
		for ri = 1:numRNodes
			fi = rNodeFNode[ri]
			@assert(fNodeToRNodeTime[fi][ri] == 0.0)
			@assert(fNodeFromRNodeTime[fi][ri] == 0.0)
		end

		# for fNodes along each rArc, fNodeFromRNodeTime and fNodeToRNodeTime should add to rArc time
		for rArcIndex = 1:numRArcs
			rArcTime = rNetTravel.arcTimes[rArcIndex]
			fNodesOnRArc = rArcFNodes[rArcIndex]
			startRNode = fNodeRNode[fNodesOnRArc[1]]
			endRNode = fNodeRNode[fNodesOnRArc[end]]
			for i = 2:length(fNodesOnRArc)-1
				fNode = fNodesOnRArc[i]
				@assert(isapprox(fNodeFromRNodeTime[fNode][startRNode] + fNodeToRNodeTime[fNode][endRNode], rArcTime; rtol = eps(Float)))
			end
		end

		if rNetTravelsLoaded == []
			# calculate shortest path data for rNetTravel
			calcRNetTravelShortestPaths!(net, rNetTravel)
		else
			rNetTravelLoaded = rNetTravelsLoaded[travelModeIndex]
			@assert(rNetTravel.arcTimes == rNetTravelLoaded.arcTimes)
			@assert(rNetTravel.arcDists == rNetTravelLoaded.arcDists)
			@assert(rNetTravel.fadjList == rNetTravelLoaded.fadjList)

			# set rNetTravel values by rNetTravelLoaded
			for fname in [:spFadjIndex, :spNodePairArcIndex, :spFadjArcList]
				setfield!(rNetTravel, fname, getfield(rNetTravelLoaded, fname))
			end

			rNetTravel.spTimes = calcRNetTravelShortestPathTimes(net, rNetTravel)
			rNetTravel.spDists = calcRNetTravelShortestPathDists(net, rNetTravel)
		end

	end
end

# For a given rNetTravel, calculate and store shortest path data
# Uses repeated Dijkstra's shortest path algorithm,
# for sparse graphs, this is faster than Floyd-Warshall algorithm
function calcRNetTravelShortestPaths!(net::Network, rNetTravel::NetTravel)
	@assert(rNetTravel.isReduced)

	# shorthand
	rGraph = net.rGraph
	rArcs = rGraph.arcs
	rArcTimes = rNetTravel.arcTimes
	fadjList = rGraph.fadjList
	numRNodes = length(rGraph.nodes)

	# create reverse graph
	revRGraph = LightGraphs.reverse(rGraph.light)

	# create arc time matrix for reverse graph, for use by LightGraphs Dijkstra SP function
	revRArcsTimes = spzeros(FloatSpTime, numRNodes, numRNodes)
	for rArc in rArcs
		# add time of this arc to revRArcsTimes if first/shortest arc connecting the two nodes
		# (there can be multiple arcs connecting a pair of nodes)
		t1 = revRArcsTimes[rArc.toNodeIndex, rArc.fromNodeIndex] # 0 if not yet filled
		t2 = rArcTimes[rArc.index]
		revRArcsTimes[rArc.toNodeIndex, rArc.fromNodeIndex] = (t1 == 0 ? t2 : min(t1, t2))
	end

	# calculate shortest paths for each possible destination node
	spData = LightGraphs.DijkstraState{Float} # init
	spTimes = rNetTravel.spTimes = Array{FloatSpTime,2}(undef, numRNodes, numRNodes)
	spSuccs = Array{IntRNode,2}(undef, numRNodes, numRNodes)
	for i = 1:numRNodes
		spData = LightGraphs.dijkstra_shortest_paths(revRGraph, i, revRArcsTimes)
		spTimes[:,i] = spData.dists
		spSuccs[:,i] = spData.parents
	end

	# change spSuccs where needed
	# where a path i -> l has spSuccs[i,l] = k, but spSuccs[i,k] = j,
	# need to change all cases where spSuccs[i,-] = k to instead have spSuccs[i,-] = j
	# this case can happen when an arc (i,k) has the same travel time as a path (i,j,k)
	for rArc in rArcs
		i = rArc.fromNodeIndex
		j = k = rArc.toNodeIndex
		while j != spSuccs[i,j]
			j = spSuccs[i,j]
		end
		if j != k
			# for any paths from i that have successor node k, change successor to be j
			spSuccs[i, spSuccs[i,:] .== k] .= j
		end
	end

	# set spNodePairArcIndex and spFadjArcList
	spNodePairArcIndex = rNetTravel.spNodePairArcIndex = spzeros(Int, numRNodes, numRNodes)
	spFadjArcList = rNetTravel.spFadjArcList = deepcopy(fadjList)
	for rArc in rArcs
		i = rArc.fromNodeIndex
		j = rArc.toNodeIndex
		t = FloatSpTime(rArcTimes[rArc.index])
		@assert(spTimes[i,j] <= t)
		if spSuccs[i,j] == j && spTimes[i,j] == t
			spNodePairArcIndex[i,j] = rArc.index
			spFadjArcList[i][findfirst(isequal(j), fadjList[i])] = rArc.index
		end
	end

	# set spFadjIndex
	spFadjIndex = rNetTravel.spFadjIndex = Array{IntFadj,2}(undef, numRNodes, numRNodes)
	for i = 1:numRNodes, j = 1:numRNodes
		spFadjIndex[i,j] = (i == j ? nullIndex : findfirst(isequal(spSuccs[i,j]), fadjList[i]))
	end

	# check stored shortest path times are same as those from traversing stored shortest paths
	# this can be slow, and should not be necessary unless shortest path code above has changed
	if checkMode && false
		checkRNetTravelShortestPathTimes(net, rNetTravel)
	end

	# check assumptions:
	for i = 1:numRNodes, j = [1:i-1;i+1:numRNodes]
		@assert(0 < spTimes[i,j] < Inf)
	end

	rNetTravel.spDists = calcRNetTravelShortestPathDists(net, rNetTravel)
end

# For a given travel mode, check that the stored shortest paths data
# gives paths with the same travel time as stored in rNetTravel.spTimes
function checkRNetTravelShortestPathTimes(net::Network, rNetTravel::NetTravel)
	# shorthand:
	numRNodes = length(net.rGraph.nodes)
	travelModeIndex = rNetTravel.modeIndex
	@unpack spTimes, spNodePairArcIndex = rNetTravel

	@assert(size(spTimes) == (numRNodes, numRNodes))

	# check each pair of start and end nodes
	for startRNode = 1:numRNodes, endRNode = 1:numRNodes
		if startRNode == endRNode
			@assert(spTimes[startRNode, endRNode] == 0.0)
			@assert(spNodePairArcIndex[startRNode, endRNode] == 0)
		else
			# follow path from start to end node, calculate and compare travel time
			i = startRNode
			t = 0.0
			while i != endRNode
				j = shortestPathNextRNode(net, travelModeIndex, i, endRNode)
				rArcIndex = shortestPathNextRArc(net, travelModeIndex, i, endRNode)
				@assert(rArcIndex == rNetTravel.spNodePairArcIndex[i,j])
				# @assert(j == net.rGraph.arcs[rArcIndex].toNodeIndex)
				t += rNetTravel.arcTimes[rArcIndex]
				i = j # go to next node
			end
			@assert(isapprox(spTimes[startRNode, endRNode], t; rtol = eps(FloatSpTime)))
		end
	end
end

# Calculate time of shortest path between each pair of nodes for rNetTravel.
# Requires some shortest path data (spFadjIndex, spNodePairArcIndex) to be calculated already.
# See also: calcRNetTravelShortestPathDists
function calcRNetTravelShortestPathTimes(net::Network, rNetTravel::NetTravel)
	@assert(rNetTravel.isReduced)

	# shorthand
	rGraph = net.rGraph
	fadjList = rGraph.fadjList
	n = length(rGraph.nodes)
	@unpack spFadjIndex, spFadjArcList, spNodePairArcIndex, arcTimes = rNetTravel

	spTimes = fill(FloatSpTime(Inf), n, n) # spTimes[i,j] = shortest path time from node i to j
	for i = 1:n spTimes[i,i] = 0 end
	for arc in rGraph.arcs
		i = arc.fromNodeIndex
		j = arc.toNodeIndex
		if spNodePairArcIndex[i,j] == arc.index
			spTimes[i,j] = arcTimes[arc.index]
		end
	end

	spNextNode = zeros(Int, n)
	spNextArcTime = zeros(Float, n)
	visited = fill(false, n) # visited[i] = true if node i has been visited and so had time populated
	path = zeros(Int, n) # path node indices
	for k = 1:n # root node of shortest path tree, find time from each node to this node
		# get data for paths towards node k
		spNextNode[k] = 0
		for j = 1:n
			if j == k continue end
			fadjIndex = spFadjIndex[j,k]
			spNextNode[j] = fadjList[j][fadjIndex]
			spNextArcTime[j] = arcTimes[spFadjArcList[j][fadjIndex]]
		end

		fill!(visited, false)
		visited[k] = true
		for i = 1:n
			if i == k || visited[i] continue end

			# trace path from i to k, or to node that already has data for path to k
			j = i
			p = 0
			while j != k && !visited[j]
				p += 1
				path[p] = j
				j = spNextNode[j]
			end

			# backtrack path, setting time from node j to k
			t = spTimes[j,k]
			while p > 0
				j = path[p]
				t += spNextArcTime[j]
				spTimes[j,k] = t
				visited[j] = true
				p -= 1
			end
		end
	end

	@assert(all(t -> t >= 0 && t != FloatSpTime(Inf), spTimes))

	return spTimes
end

# Calculate distance of shortest path between each pair of nodes for rNetTravel.
# Requires some shortest path data (spFadjIndex, spNodePairArcIndex) to be calculated already.
# See also: calcRNetTravelShortestPathTimes
function calcRNetTravelShortestPathDists(net::Network, rNetTravel::NetTravel)
	@assert(rNetTravel.isReduced)

	# shorthand
	rGraph = net.rGraph
	fadjList = rGraph.fadjList
	n = length(rGraph.nodes)
	@unpack spFadjIndex, spFadjArcList, spNodePairArcIndex, arcDists = rNetTravel

	spDists = fill(FloatSpDist(Inf), n, n) # spDists[i,j] = shortest path distance from node i to j
	for i = 1:n spDists[i,i] = 0 end
	for arc in rGraph.arcs
		i = arc.fromNodeIndex
		j = arc.toNodeIndex
		if spNodePairArcIndex[i,j] == arc.index
			spDists[i,j] = arc.distance
		end
	end

	spNextNode = zeros(Int, n)
	spNextArcDist = zeros(Float, n)
	visited = fill(false, n) # visited[i] = true if node i has been visited and so had distance populated
	path = zeros(Int, n) # path node indices
	for k = 1:n # root node of shortest path tree, find distance from each node to this node
		# get data for paths towards node k
		spNextNode[k] = 0
		for j = 1:n
			if j == k continue end
			fadjIndex = spFadjIndex[j,k]
			spNextNode[j] = fadjList[j][fadjIndex]
			spNextArcDist[j] = arcDists[spFadjArcList[j][fadjIndex]]
		end

		fill!(visited, false)
		visited[k] = true
		for i = 1:n
			if i == k || visited[i] continue end

			# trace path from i to k, or to node that already has data for path to k
			j = i
			p = 0
			while j != k && !visited[j]
				p += 1
				path[p] = j
				j = spNextNode[j]
			end

			# backtrack path, setting distance from node j to k
			d = spDists[j,k]
			while p > 0
				j = path[p]
				d += spNextArcDist[j]
				spDists[j,k] = d
				visited[j] = true
				p -= 1
			end
		end
	end

	@assert(all(d -> d >= 0 && d != FloatSpDist(Inf), spDists))

	return spDists
end

# for the shortest path from startRNode to endRNode,
# travelling with the given travel mode,
# return the index of the successor rNode of startRNode
function shortestPathNextRNode(net::Network, travelModeIndex::Int, startRNode::Int, endRNode::Int)
	@assert(startRNode != endRNode)

	rNetTravel = net.rNetTravels[travelModeIndex] # shorthand
	fadjIndex = rNetTravel.spFadjIndex[startRNode, endRNode]
	return net.rGraph.fadjList[startRNode][fadjIndex]
end

# for the shortest path from startRNode to endRNode,
# travelling with the given travel mode,
# return the index of the rArc to travel along from startRNode
function shortestPathNextRArc(net::Network, travelModeIndex::Int, startRNode::Int, endRNode::Int)
	@assert(startRNode != endRNode)

	rNetTravel = net.rNetTravels[travelModeIndex] # shorthand
	fadjIndex = rNetTravel.spFadjIndex[startRNode, endRNode]
	return rNetTravel.spFadjArcList[startRNode][fadjIndex]
end

# Given the travel mode, and indices of start and end node (in full graph),
# returns travel time for shortest path, and the first and last rNodes in the path (if any)
function shortestPathData(net::Network, travelModeIndex::Int, startFNode::Int, endFNode::Int)

	# shorthand:
	fNetTravel = net.fNetTravels[travelModeIndex]
	rNetTravel = net.rNetTravels[travelModeIndex]
	@unpack fNodeFromRNodeTime, fNodeToRNodeTime = fNetTravel
	@unpack fNodeRArcs, fNodeToRNodes, fNodeFromRNodes, isFNodeCommon, fNodeCommonFNodeIndex = net

	# return values:
	shortestTravelTime = Inf
	rNodes = (nullIndex, nullIndex) # first and last rNodes in path

	if startFNode == endFNode
		shortestTravelTime = 0.0
		rNodes = (nullIndex, nullIndex)
		return shortestTravelTime, rNodes
	end

	# if startFNode or endFNode are common, can look up stored shortest path data
	if isFNodeCommon[startFNode]
		i = fNodeCommonFNodeIndex[startFNode]
		shortestTravelTime = fNetTravel.commonFNodeToFNodeTime[i, endFNode]
		rNodes = fNetTravel.commonFNodeToFNodeRNodes[i, endFNode]
		return shortestTravelTime, rNodes
	elseif isFNodeCommon[endFNode]
		i = fNodeCommonFNodeIndex[endFNode]
		shortestTravelTime = fNetTravel.fNodeToCommonFNodeTime[startFNode, i]
		rNodes = fNetTravel.fNodeToCommonFNodeRNodes[startFNode, i]
		return shortestTravelTime, rNodes
	end

	# see if travel is possible without using any nodes in rGraph
	if !isFNodeInRGraph(net, startFNode) && !isFNodeInRGraph(net, endFNode) && fNodeRArcs[startFNode] == fNodeRArcs[endFNode]
		rArc = findRArcFromFNodeToFNode(net, startFNode, endFNode)
		if rArc != nullIndex
			# calculate travel time on rArc
			fromRNode = net.rGraph.arcs[rArc].fromNodeIndex
			toRNode = net.rGraph.arcs[rArc].toNodeIndex
			travelTime = fNodeFromRNodeTime[endFNode][fromRNode] - fNodeFromRNodeTime[startFNode][fromRNode]
			@assert(travelTime > 0)
			@assert(isapprox(travelTime + fNodeToRNodeTime[endFNode][toRNode], fNodeToRNodeTime[startFNode][toRNode]; rtol = eps(Float)))

			shortestTravelTime = travelTime # not necessarily shortest travel time, still need to check paths that use rNodes
			rNodes = (nullIndex, nullIndex)
		end
	end

	# for each pair of starting and ending rNodes (near given fNodes), calculate travel time
	for startRNode = fNodeToRNodes[startFNode], endRNode = fNodeFromRNodes[endFNode]
		travelTime = rNetTravel.spTimes[startRNode, endRNode]
		travelTime += fNodeToRNodeTime[startFNode][startRNode]
		travelTime += fNodeFromRNodeTime[endFNode][endRNode]
		if travelTime < shortestTravelTime
			shortestTravelTime = travelTime
			rNodes = (startRNode, endRNode)
		end
	end

	return shortestTravelTime, rNodes
end

# Returns the travel time of the shortest path between two nodes; see shortestPathData function.
function shortestPathTravelTime(net::Network, travelModeIndex::Int, startFNode::Int, endFNode::Int)
	return shortestPathData(net, travelModeIndex, startFNode, endFNode)[1]
end

# Returns distance of shortest path (shortest by time) between two nodes.
# The startRNode and endRNode of shortest path can be given as kwargs if known.
function shortestPathDistance(net::Network, travelModeIndex::Int, startFNode::Int, endFNode::Int;
	startRNode::Int = nullIndex - 1, endRNode::Int = nullIndex - 1)
	# note that the default value of startRNode and endRNode are nullIndex - 1 and not just nullIndex;
	# this is to allow them to be set to nullIndex which indicates that the shortest path has no rNodes.

	# shorthand:
	rArcs = net.rGraph.arcs
	@unpack fNodeToRNodeDist, fNodeFromRNodeDist, isFNodeCommon, fNodeCommonFNodeIndex = net
	fNetTravel = net.fNetTravels[travelModeIndex]
	rNetTravel = net.rNetTravels[travelModeIndex]

	dist = 0.0
	if startFNode == endFNode return dist end

	# if startFNode or endFNode are common, can look up stored shortest path distance
	if isFNodeCommon[startFNode]
		i = fNodeCommonFNodeIndex[startFNode]
		return fNetTravel.commonFNodeToFNodeDist[i, endFNode]
	elseif isFNodeCommon[endFNode]
		i = fNodeCommonFNodeIndex[endFNode]
		return fNetTravel.fNodeToCommonFNodeDist[startFNode, i]
	end

	if startRNode == nullIndex - 1 || endRNode == nullIndex - 1
		@assert(startRNode == nullIndex - 1 && endRNode == nullIndex - 1)
		# get start and end rNodes
		(shortestTravelTime, rNodes) = shortestPathData(net, travelModeIndex, startFNode, endFNode)
		(startRNode, endRNode) = rNodes
	elseif startRNode == nullIndex || endRNode == nullIndex
		@assert(startRNode == nullIndex && endRNode == nullIndex)
	else
		@assert(in(startRNode, net.fNodeToRNodes[startFNode]))
		@assert(in(endRNode, net.fNodeFromRNodes[endFNode]))
	end

	if startRNode == nullIndex
		# shortest path does not use any rNodes
		@assert(endRNode == nullIndex)
		rArcIndex = findRArcFromFNodeToFNode(net, startFNode, endFNode)
		@assert(rArcIndex != nullIndex) # otherwise would need to use an rNode
		rArc = rArcs[rArcIndex]
		fromRNode = rArc.fromNodeIndex
		toRNode = rArc.toNodeIndex
		dist = fNodeFromRNodeDist[endFNode][fromRNode] - fNodeFromRNodeDist[startFNode][fromRNode]
		@assert(dist >= 0)
		@assert(isapprox(dist + fNodeToRNodeDist[endFNode][toRNode], fNodeToRNodeDist[startFNode][toRNode]; rtol = eps(Float)))
	else
		# shortest path uses at least one rNode
		@assert(startRNode != nullIndex && endRNode != nullIndex)
		dist += net.fNodeToRNodeDist[startFNode][startRNode]
		dist += rNetTravel.spDists[startRNode, endRNode]
		dist += net.fNodeFromRNodeDist[endFNode][endRNode]
	end

	return dist
end

# Find and return shortest path (as represented by list of nodes) from startFNode to endFNode,
# starting at startFNode at startTime, with a given travel mode
# Also return time at which each node in path will be reached
function shortestPath(net::Network, travelModeIndex::Int, startFNode::Int, endFNode::Int, startTime::Float)

	# shorthand:
	fNetTravel = net.fNetTravels[travelModeIndex]
	rNetTravel = net.rNetTravels[travelModeIndex]
	@unpack fNodeToRNodeTime, fNodeFromRNodeTime, rArcFNodesTimes = fNetTravel
	@unpack fGraph, rGraph, rArcFNodes, rNodeFNode, fNodeRArcs, fNodeToRNodeNextFNode = net
	rArcs = rGraph.arcs

	# return values:
	fNodeList = Vector{Int}() # list of indices of nodes to travel through, in traversal order
	fNodeTimeList = Vector{Float}() # time at which each node in fNodeList will be reached

	if startFNode == endFNode
		fNodeList = [startFNode]
		fNodeTimeList = [startTime]
		return fNodeList, fNodeTimeList
	end

	# find which rNodes are used in the shortest path
	(travelTime, rNodes) = shortestPathData(net, travelModeIndex, startFNode, endFNode)
	startRNode = rNodes[1]
	endRNode = rNodes[2]

	if startRNode == nullIndex
		@assert(endRNode == nullIndex)
		# no nodes from rGraph are used, travel along a single rArc from startFNode to endFNode
		rArcIndices = fNodeRArcs[startFNode]
		@assert(length(rArcIndices) <= 2)
		toRNode = rArcs[rArcIndices[1]].toNodeIndex
		if fNodeToRNodeTime[startFNode][toRNode] - fNodeToRNodeTime[endFNode][toRNode] < 0 # negative travel time
			toRNode = rArcs[rArcIndices[2]].toNodeIndex
		end
		# create fNodeList, fNodeTimeList
		fNodeList = [startFNode]
		fNodeTimeList = [startTime]
		fNode = startFNode
		currentTime = startTime
		while fNode != endFNode
			# find next arc, add next node and arc time
			nextFNode = fNodeToRNodeNextFNode[fNode][toRNode]
			arcTime = fNetTravel.arcTimes[fGraph.nodePairArcIndex[fNode, nextFNode]]
			currentTime += arcTime
			push!(fNodeList, nextFNode)
			push!(fNodeTimeList, currentTime)
			fNode = nextFNode # go to next node
		end
		return fNodeList, fNodeTimeList
	end

	# travel using rGraph
	if startRNode != nullIndex
		@assert(endRNode != nullIndex)
		currentTime = startTime

		# find path from startFNode to startRNode
		fNodeList = [startFNode]
		fNodeTimeList = [startTime]
		fNode = startFNode
		startRNodeTime = startTime + fNodeToRNodeTime[startFNode][startRNode] # time to reach first rNode
		while fNode != rNodeFNode[startRNode]
			# find next arc, add next node and arc time
			nextFNode = fNodeToRNodeNextFNode[fNode][startRNode]
			# arcTime = fNodeToRNodeTime[fNode][startRNode] - fNodeToRNodeTime[nextFNode][startRNode]
			# arcTime = fNetTravel.arcTimes[fGraph.nodePairArcIndex[fNode, nextFNode]]
			# currentTime += arcTime
			currentTime = startRNodeTime - fNodeToRNodeTime[nextFNode][startRNode]
			push!(fNodeList, nextFNode)
			push!(fNodeTimeList, currentTime)
			fNode = nextFNode # go to next node
		end

		# find path from after startRNode to endRNode (rNodes[1] to rNodes[2])
		rNode = startRNode
		while rNode != endRNode
			# find next rArc, add fNodes
			rArc = shortestPathNextRArc(net, travelModeIndex, rNode, endRNode)
			nextRNode = rArcs[rArc].toNodeIndex
			fNodesTimesForRArc = rArcFNodesTimes[rArc]
			rNodeTime = currentTime
			fNodesOnRArc = rArcFNodes[rArc]
			i = 2
			while i <= length(fNodesOnRArc)
				fNode = fNodesOnRArc[i]
				currentTime = rNodeTime + fNodesTimesForRArc[i]
				push!(fNodeList, fNode)
				push!(fNodeTimeList, currentTime)
				i += 1 # go to next node
			end
			rNode = nextRNode # go to next node
		end

		# find path from after endRNode to endFNode
		fNode = rNodeFNode[endRNode]
		if fNode != endFNode
			rArcIndices = fNodeRArcs[endFNode]
			@assert(length(rArcIndices) <= 2) # endFNode should only be on at most two rArcs, as it is not in rGraph
			rArc = nullIndex
			for rArcIndex in rArcIndices
				if rArcFNodes[rArcIndex][1] == fNode
					rArc = rArcIndex
				end
			end
			fNodesOnRArc = rArcFNodes[rArc]
			fNodeTime = currentTime
			i = 1
			while fNodesOnRArc[i] != endFNode
				i += 1 # go to next node on rArc
				fNode = fNodesOnRArc[i]
				currentTime = fNodeTime + fNodeFromRNodeTime[fNode][endRNode]
				push!(fNodeList, fNode)
				push!(fNodeTimeList, currentTime)
			end
		end
	end

	return fNodeList, fNodeTimeList
end

# set the common fNodes of the network,
# and the shortest path travel data between each fNode and commonFNodes
function setCommonFNodes!(net::Network, commonFNodes::Vector{Int})
	numFNodes = length(net.fGraph.nodes) # shorthand

	# add common fNodes to net
	commonFNodes = sort(unique(commonFNodes))
	@assert(all(commonFNode -> 1 <= commonFNode <= numFNodes, commonFNodes))
	numCommonFNodes = length(commonFNodes) # shorthand
	net.isFNodeCommon = [false for i = 1:numFNodes]
	fNodeCommonFNodeIndex = [nullIndex for i = 1:numFNodes]
	for (i, commonFNode) in enumerate(commonFNodes)
		fNodeCommonFNodeIndex[commonFNode] = i
	end

	# calculate and store shortest path travel data between all fNodes and commonFNodes
	for fNetTravel in net.fNetTravels
		fNetTravel.commonFNodeToFNodeTime = Array{Float,2}(undef, numCommonFNodes, numFNodes)
		fNetTravel.commonFNodeToFNodeDist = Array{Float,2}(undef, numCommonFNodes, numFNodes)
		fNetTravel.commonFNodeToFNodeRNodes = Array{Tuple{Int,Int},2}(undef, numCommonFNodes, numFNodes)

		fNetTravel.fNodeToCommonFNodeTime = Array{Float,2}(undef, numFNodes, numCommonFNodes)
		fNetTravel.fNodeToCommonFNodeDist = Array{Float,2}(undef, numFNodes, numCommonFNodes)
		fNetTravel.fNodeToCommonFNodeRNodes = Array{Tuple{Int,Int},2}(undef, numFNodes, numCommonFNodes)

		for commonFNode in commonFNodes, fNode = 1:numFNodes
			i = fNodeCommonFNodeIndex[commonFNode]

			(travelTime, rNodes) = shortestPathData(net, fNetTravel.modeIndex, commonFNode, fNode)
			fNetTravel.commonFNodeToFNodeTime[i,fNode] = travelTime
			fNetTravel.commonFNodeToFNodeRNodes[i,fNode] = rNodes
			fNetTravel.commonFNodeToFNodeDist[i,fNode] = shortestPathDistance(net, fNetTravel.modeIndex, commonFNode, fNode; startRNode = rNodes[1], endRNode = rNodes[2])

			(travelTime, rNodes) = shortestPathData(net, fNetTravel.modeIndex, fNode, commonFNode)
			fNetTravel.fNodeToCommonFNodeTime[fNode,i] = travelTime
			fNetTravel.fNodeToCommonFNodeRNodes[fNode,i] = rNodes
			fNetTravel.fNodeToCommonFNodeDist[fNode,i] = shortestPathDistance(net, fNetTravel.modeIndex, fNode, commonFNode; startRNode = rNodes[1], endRNode = rNodes[2])
		end
	end

	# add common fNodes to net
	# need to do this after storing results from shortest path calculations above,
	# otherwise the shortestPathData function may try to use values before they are initialised.
	net.commonFNodes = commonFNodes
	net.fNodeCommonFNodeIndex = fNodeCommonFNodeIndex
	for commonFNode in commonFNodes
		net.isFNodeCommon[commonFNode] = true
	end
end

# Return index of the unique rArc that two fNodes have in common,
# where the rArc contains fromFNode before toFNode
# return nullIndex if no such arc found
function findRArcFromFNodeToFNode(net::Network, fromFNode::Int, toFNode::Int)
	@assert(fromFNode != toFNode)
	@assert(isFNodeInRGraph(net, fromFNode) + isFNodeInRGraph(net, toFNode) < 2) # otherwise result may not be unique, as multiple rArcs can connect two rNodes

	rArc = nullIndex
	numRArcsFound = 0
	for rArcIndex in net.fNodeRArcs[fromFNode], rArcIndex2 in net.fNodeRArcs[toFNode]
		if rArcIndex == rArcIndex2 && net.rArcFNodeIndex[rArcIndex][fromFNode] < net.rArcFNodeIndex[rArcIndex][toFNode]
			rArc = rArcIndex
			numRArcsFound += 1
		end
	end
	@assert(numRArcsFound <= 1) # result should be unique

	return rArc
end


text:
function pixelshufflehead(k_in, k_out; n_upscale = 3, k_mid = 64, σ = relu)
    layers = []

    push!(layers, Chain(
        PixelShuffle(2, k_in, k_mid),
        BatchNorm(k_mid, σ),
    ))
    for i = 2:n_upscale
        push!(layers, Chain(
            PixelShuffle(2, k_mid),
            BatchNorm(k_mid, σ),

        ))
    end

    push!(layers, Conv((1, 1), k_mid => k_out))

    return Chain(layers...)
end

ext: jl


text:
using LULESH
using MPI
# using Enzyme

# Enzyme.API.printperf!(true)
# Enzyme.API.printall!(true)
# Enzyme.API.instname!(true)

# Enzyme.API.inlineall!(true)
# Enzyme.API.maxtypeoffset!(32)
# ccall((:EnzymeSetCLInteger, Enzyme.API.libEnzyme), Cvoid, (Ptr{Cvoid}, Int64), cglobal((:MaxTypeOffset, Enzyme.API.libEnzyme)), 32)

function main(nx, structured, num_iters, mpi)
    # TODO: change default nr to 11
    nr = 1
    balance = 1
    cost = 1
    floattype = Float64
    devicetype = Vector

    if mpi
        !MPI.Initialized() && MPI.Init()
        comm = MPI.COMM_WORLD
    else
        comm = nothing
    end

    prob = LuleshProblem(num_iters, structured, nx, nr, balance, cost, devicetype, floattype, comm)

    # Set up the mesh and decompose. Assumes regular cubes for now
    # TODO: modify this constructor to account for new fields
    # TODO: setup communication buffers

    domain = Domain(prob)
    shadowDomain = Domain(prob)

    # getnodalMass = nodalMass(domain)

    # Initial domain boundary communication
    # commRecv(domain, MSG_COMM_SBN, 1,
    #             domain.sizeX + 1, domain.sizeY + 1, domain.sizeZ + 1,
    #             true, false, prob.comm)
    #    CommSend<&Domain::nodalMass>(*domain, MSG_COMM_SBN,
    #             domain->sizeX() + 1, domain->sizeY() + 1, domain->sizeZ() +  1,
    #             true, false)
    #    CommSBN<&Domain::nodalMass>(*domain)

    # End initialization
    if mpi
        MPI.Barrier(prob.comm::MPI.Comm)
    end

    if getMyRank(prob.comm) == 0
        if (structured)
            @info "Running" until=domain.stoptime domain=(nx,nx,nx)
        else
            @info "Running" until=domain.stoptime domain=domain.numElem
            @warn "Unstructured setup not supported"
        end
    end

    # timestep to solution
    start = getWtime(prob.comm)
    while domain.time < domain.stoptime
        # this has been moved after computation of volume forces to hide launch latencies
        timeIncrement!(domain)
        lagrangeLeapFrog(domain)

        # checkErrors(domain, its, myRank)
        if getMyRank(prob.comm) == 0
            @info "Completed" cycle=domain.cycle time=domain.time dt=domain.deltatime
        end
        if domain.cycle >= num_iters
            break
        end
    end

    # Use reduced max elapsed time
    elapsed_time = getWtime(prob.comm) - start
    elapsed_timeG = comm_max(elapsed_time, prob.comm)

    if getMyRank(prob.comm) == 0
        verifyAndWriteFinalOutput(elapsed_timeG, domain, nx, getNumRanks(prob.comm))
    end

    if mpi
        MPI.Finalize()
    end
end

if !isinteractive()
    args = LULESH.parse_cmd()

    # assume cube subdomain geometry for now (nx)
    nx = args["N"]
    structured = args["s"]
    num_iters = args["num_iters"]
    mpi = args["mpi"]
    main(nx, structured, num_iters, mpi)
end


text:
using ProcessMining
using Test

# @time ProcessMining.read_xes(joinpath("data", "Performance.xes"))


text:
using Test, DataLoaders

include("dataset.jl")

@testset "RandomSampler" begin
  dataset = getdataset((10,10), 100)
  sampler = RandomSampler(dataset)

  @test sampler.dataset     == dataset
  @test length(sampler)     == length(dataset)
  @test IndexStyle(sampler) == IndexLinear()
end


text:
function  householder{T<:FloatingPoint}(a::Array{T,1})

  n = length(a)
  m = 1

  q = eye(n)
  r = copy(a)

  for i = 1:m

    x = r[i:n,i]
    e = zeros(n-i+1)
    e[1] = norm(x)
    u = x-e
    u = u/norm(u)
    qq = eye(n-i+1)-2*u*u'
    r[i:n,:] = qq*r[i:n,:]
    q[i:n,:] = qq*q[i:n,:]

  end

  return q',r

end

function  householder{T<:FloatingPoint}(a::Array{Complex{T},1})

  n = length(a)
  m = 1

  q = complex(eye(n))
  r = copy(a)

  for i = 1:m

    x = r[i:n,i]
    e = complex(zeros(n-i+1))
    e[1] = norm(x)
    u = x-e
    u = u/norm(u)
    qq = complex(eye(n-i+1))-(complex(1.0)+(x'u)/(u'x)).*u*u'
    r[i:n,:] = qq*r[i:n,:]
    q[i:n,:] = qq*q[i:n,:]

  end

  return q',r

end

function  householder{T<:FloatingPoint}(a::Array{T,2})

  (n,m) = size(a)

  q = eye(n)
  r = copy(a)

  for i = 1:m

    x = r[i:n,i]
    e = zeros(n-i+1)
    e[1] = norm(x)
    u = x-e
    u = u/norm(u)
    qq = eye(n-i+1)-2*u*u'
    r[i:n,:] = qq*r[i:n,:]
    q[i:n,:] = qq*q[i:n,:]

  end

  return q',r

end

function  householder{T<:FloatingPoint}(a::Array{Complex{T},2})

  (n,m) = size(a)

  q = complex(eye(n))
  r = copy(a)

  for i = 1:m

    x = r[i:n,i]
    e = complex(zeros(n-i+1))
    e[1] = norm(x)
    u = x-e
    u = u/norm(u)
    qq = complex(eye(n-i+1))-(complex(1.0)+(x'u)/(u'x)).*u*u'
    r[i:n,:] = qq*r[i:n,:]
    q[i:n,:] = qq*q[i:n,:]

  end

  return q',r

end

max_stars_repo_name: dingliumath/SolveDSGE

max_issues_repo_name: dingliumath/SolveDSGE

max_forks_repo_name: dingliumath/SolveDSGE


text:
module CasesLH

using DocStringExtensions

# include("types.jl");
include("modifier.jl");
include("base_case.jl"); # wrong name +++

end


text:
using Abaco
using DelimitedFiles
using Test

datadir=joinpath(@__DIR__,"data")

expected_triggers = 6
actual_triggers = 0

ts = 1635601415
sn1 = "mi01"
sn2 = "mi02"

values = Dict(
    "x" => 2.0,
    "y" => 4.0,
    "z" => 8.0
    )

df = readdlm(joinpath(datadir,"formulas.csv"), ';', header=true)

function onresult(ts, ne, name, value, inputs)
    global actual_triggers
    @debug "age [$ts]: ne: [$ne] $name = $value"
    actual_triggers += 1
end

abaco = abaco_init(onresult, interval=-1)

setup_settings(abaco, "sensor", oncomplete=onresult)
setup_settings(abaco, "hub", oncomplete=onresult)

add_formulas(abaco, df[1])

city = node(abaco, "milano", "hub")

#sensor = node(abaco, ne, "sensor")
node(abaco, city, sn1, "sensor")
node(abaco, city, sn2, "sensor")

ingest(abaco, ts, sn1, values)
ingest(abaco, ts, sn2, values)

#ingest(abaco, ts, sn1, Dict("x"=>4))

#@info "elements: $(abaco.node)"
#@info "hub: $(abaco.node["milano"])"
#@info "hub settings: $(abaco.cfg["hub"])"
#@info "mi01: $(abaco.node["mi01"])"
#@info "origins: $(abaco.origins)"
#@info "target: $(abaco.target)"

vals = Abaco.origins_vals(abaco, "milano", ts)
@debug "origin_vals: $vals"

@test actual_triggers === expected_triggers


text:
# ---
# title: 54. Spiral Matrix
# id: problem54
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/spiral-matrix/description/>
# hidden: true
# ---
#
# Given an `m x n` `matrix`, return _all elements of the_ `matrix` _in spiral
# order_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/13/spiral1.jpg)
#
#
#
#     Input: matrix = [[1,2,3],[4,5,6],[7,8,9]]
#     Output: [1,2,3,6,9,8,7,4,5]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/13/spiral.jpg)
#
#
#
#     Input: matrix = [[1,2,3,4],[5,6,7,8],[9,10,11,12]]
#     Output: [1,2,3,4,8,12,11,10,9,5,6,7]
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[i].length`
#   * `1 <= m, n <= 10`
#   * `-100 <= matrix[i][j] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
struct 𝕊{n} <: Manifold end

dim(::Type{𝕊{n}}) where n = n
embdim(::Type{𝕊{n}}) where n = n+1

# Creates aliases S1, TS1, etc.
for n = 1:3
    @eval const $(Symbol("S$n")) = 𝕊{$n}
    @eval const $(Symbol("TS$n")) = T{𝕊{$n}}
    @eval export $(Symbol("S$n"))
    @eval export $(Symbol("TS$n"))
end
TS{n} = T{𝕊{n}}

# Sterographic Projection
abstract type SterProj <: ChartID end
isglobal(::Chart{<:SterProj,𝕊{n}}) where n = false
struct SterProjSouth <: SterProj end  # South pole
struct SterProjNorth <: SterProj end  # North pole
emb_to_chart(pe, ::Chart{SterProjSouth,𝕊{n}}) where n = SVector{n}(pe[1:n])./(1+pe[n+1])
emb_to_chart(pe, ::Chart{SterProjNorth,𝕊{n}}) where n = SVector{n}(pe[1:n])./(1-pe[n+1])
function chart_to_emb(pn, ::Chart{SterProjSouth,S2})
    den = 1 + pn[1]^2 + pn[2]^2
    SA[2*pn[1]/den, 2*pn[2]/den, (2-den)/den]
end
function chart_to_emb(ps, ::Chart{SterProjNorth,S2})
    den = 1 + ps[1]^2 + ps[2]^2
    SA[2*ps[1]/den, 2*ps[2]/den, (den-2)/den]
end
chart_choice_coord_rep(::Chart{<:SterProj,S2}) = EmbRep()
choose_chart_emb(::EmbRep, pe, ::Chart{<:SterProj,S2}) =
    (pe[3] < 0)[1] ? Chart{SterProjSouth,S2}() : Chart{SterProjNorth,S2}()

# Same both ways
chart_transition(ps, ::Chart{SterProjSouth,𝕊{n}}, ::Chart{SterProjNorth,𝕊{n}}) where n =
    ps/norm(ps)^2
chart_transition(pn, ::Chart{SterProjNorth,𝕊{n}}, ::Chart{SterProjSouth,𝕊{n}}) where n =
    pn/norm(pn)^2

# Angle representation for S1
abstract type AngleChart <: ChartID end
isglobal(::Chart{<:AngleChart,S1}) = false
struct Angleπ <: AngleChart end    # Angles in (-π, π)
struct Angle2π <: AngleChart end   # Angles in (0, 2π)
default_chart(::Type{S1}) = Chart{Angleπ,S1}()
emb_to_chart(pe, ::Chart{Angleπ,S1}) = SA[atan(pe[2],pe[1])]
function emb_to_chart(pe, C2π::Chart{Angle2π,S1})
    Cπ = Chart{Angleπ,S1}()
    pπ = emb_to_chart(pe, Cπ) # Between -π and π
    chart_transition(pπ, Cπ, C2π)
end

chart_to_emb(p, ::Chart{<:AngleChart,S1}) = [cos.(p); sin.(p)]

chart_transition(p2π, ::Chart{Angle2π,S1}, ::Chart{Angleπ,S1}) = (p2π .< π)[1] ? p2π : p2π .- 2π
chart_transition(pπ, ::Chart{Angleπ,S1}, ::Chart{Angle2π,S1}) = (pπ .> 0.)[1] ? pπ : pπ .+ 2π

chart_choice_coord_rep(::Chart{<:AngleChart,S1}) = EmbRep()
choose_chart_emb(::EmbRep, pe, ::Chart{<:AngleChart,S1}) =
    (pe[1] .< 0.)[1] ? Chart{Angle2π,S1}() : Chart{Angleπ,S1}()


text:
#####################################################################
#                                                                   #
# This file provides functions for interfacing with pti .raw files  #
#                                                                   #
#####################################################################

"""
_get_pti_sections()

Internal function. Returns `Array` of the names of the sections, in the order
that they appear in a PTI file, v33
"""
function _get_pti_sections()::Array
    return ["CASE IDENTIFICATION", "BUS", "LOAD", "FIXED SHUNT", "GENERATOR", "BRANCH", "TRANSFORMER",
            "AREA INTERCHANGE", "TWO-TERMINAL DC", "VOLTAGE SOURCE CONVERTER", "IMPEDANCE CORRECTION",
            "MULTI-TERMINAL DC", "MULTI-SECTION LINE", "ZONE", "INTER-AREA TRANSFER", "OWNER",
            "FACTS CONTROL DEVICE", "SWITCHED SHUNT", "GNE DEVICE", "INDUCTION MACHINE"]
end

"""
    _get_pti_dtypes(field_name)

Internal function. Returns array of data types for PTI file section given by
`field_name`, as enumerated by PSS/E Program Operation Manual.
"""
function _get_pti_dtypes(field_name::AbstractString)::Array
    transaction_dtypes = [("IC", Int64), ("SBASE", Float64), ("REV", Int64),
                          ("XFRRAT", Float64), ("NXFRAT", Float64), ("BASFRQ", Float64)]

    bus_dtypes = [("I", Int64), ("NAME", String), ("BASKV", Float64),
                  ("IDE", Int64), ("AREA", Int64), ("ZONE", Int64), ("OWNER", Int64),
                  ("VM", Float64), ("VA", Float64), ("NVHI", Float64), ("NVLO", Float64),
                  ("EVHI", Float64), ("EVLO", Float64)]

    load_dtypes = [("I", Int64), ("ID", String), ("STATUS", Int64), ("AREA", Int64),
                   ("ZONE", Int64), ("PL", Float64), ("QL", Float64), ("IP", Float64),
                   ("IQ", Float64), ("YP", Float64), ("YQ", Float64), ("OWNER", Int64),
                   ("SCALE", Int64), ("INTRPT", Int64)]

    fixed_shunt_dtypes = [("I", Int64), ("ID", String), ("STATUS", Int64), ("GL", Float64),
                           ("BL", Float64)]

    generator_dtypes = [("I", Int64), ("ID", String), ("PG", Float64), ("QG", Float64),
                        ("QT", Float64), ("QB", Float64), ("VS", Float64), ("IREG", Int64),
                        ("MBASE", Float64), ("ZR", Float64), ("ZX", Float64), ("RT", Float64),
                        ("XT", Float64), ("GTAP", Float64), ("STAT", Int64), ("RMPCT", Float64),
                        ("PT", Float64), ("PB", Float64), ("O1", Int64), ("F1", Float64),
                        ("O2", Int64), ("F2", Float64), ("O3", Int64), ("F3", Float64),
                        ("O4", Int64), ("F4", Float64), ("WMOD", Int64), ("WPF", Float64)]

    branch_dtypes = [("I", Int64), ("J", Int64), ("CKT", String), ("R", Float64), ("X", Float64),
                     ("B", Float64), ("RATEA", Float64), ("RATEB", Float64), ("RATEC", Float64),
                     ("GI", Float64), ("BI", Float64), ("GJ", Float64), ("BJ", Float64),
                     ("ST", Int64), ("MET", Int64), ("LEN", Float64), ("O1", Int64), ("F1", Float64),
                     ("O2", Int64), ("F2", Float64), ("O3", Int64), ("F3", Float64),
                     ("O4", Int64), ("F4", Float64)]

                     transformer_dtypes = [("I", Int64), ("J", Int64), ("K", Int64), ("CKT", String),
                     ("CW", Int64), ("CZ", Int64), ("CM", Int64), ("MAG1", Float64),
                     ("MAG2", Float64), ("NMETR", Int64), ("NAME", String),
                     ("STAT", Int64),
                     ("O1", Int64), ("F1", Float64),
                     ("O2", Int64), ("F2", Float64),
                     ("O3", Int64), ("F3", Float64),
                     ("O4", Int64), ("F4", Float64),
                     ("VECGRP", String)]

transformer_3_1_dtypes = [("R1-2", Float64), ("X1-2", Float64), ("SBASE1-2", Float64),
                         ("R2-3", Float64), ("X2-3", Float64), ("SBASE2-3", Float64),
                         ("R3-1", Float64), ("X3-1", Float64), ("SBASE3-1", Float64),
                         ("VMSTAR", Float64), ("ANSTAR", Float64)]

transformer_3_2_dtypes = [("WINDV1", Float64), ("NOMV1", Float64), ("ANG1", Float64),
                         ("RATA1", Float64), ("RATB1", Float64), ("RATC1", Float64),
                         ("COD1", Int64), ("CONT1", Int64), ("RMA1", Float64), ("RMI1", Float64),
                         ("VMA1", Float64), ("VMI1", Float64), ("NTP1", Float64), ("TAB1", Int64),
                         ("CR1", Float64), ("CX1", Float64), ("CNXA1", Float64)]

transformer_3_3_dtypes = [("WINDV2", Float64), ("NOMV2", Float64), ("ANG2", Float64),
                         ("RATA2", Float64), ("RATB2", Float64), ("RATC2", Float64),
                         ("COD2", Int64), ("CONT2", Int64), ("RMA2", Float64), ("RMI2", Float64),
                         ("VMA2", Float64), ("VMI2", Float64), ("NTP2", Float64), ("TAB2", Int64),
                         ("CR2", Float64), ("CX2", Float64), ("CNXA2", Float64)]

transformer_3_4_dtypes = [("WINDV3", Float64), ("NOMV3", Float64), ("ANG3", Float64),
                         ("RATA3", Float64), ("RATB3", Float64), ("RATC3", Float64),
                         ("COD3", Int64), ("CONT3", Int64), ("RMA3", Float64), ("RMI3", Float64),
                         ("VMA3", Float64), ("VMI3", Float64), ("NTP3", Float64), ("TAB3", Int64),
                         ("CR3", Float64), ("CX3", Float64), ("CNXA3", Float64)]

transformer_2_1_dtypes = [("R1-2", Float64), ("X1-2", Float64), ("SBASE1-2", Float64)]

transformer_2_2_dtypes = [("WINDV1", Float64), ("NOMV1", Float64), ("ANG1", Float64),
                         ("RATA1", Float64), ("RATB1", Float64), ("RATC1", Float64),
                         ("COD1", Int64), ("CONT1", Int64), ("RMA1", Float64), ("RMI1", Float64),
                         ("VMA1", Float64), ("VMI1", Float64), ("NTP1", Float64), ("TAB1", Int64),
                         ("CR1", Float64), ("CX1", Float64), ("CNXA1", Float64)]

transformer_2_3_dtypes = [("WINDV2", Float64), ("NOMV2", Float64)]

    area_interchange_dtypes = [("I", Int64), ("ISW", Int64),
                               ("PDES", Float64), ("PTOL", Float64),
                               ("ARNAME", String)]

    two_terminal_line_dtypes = [("NAME", String), ("MDC", Int64), ("RDC", Float64),
                                ("SETVL", Float64), ("VSCHD", Float64), ("VCMOD", Float64),
                                ("RCOMP", Float64), ("DELTI", Float64), ("METER", String),
                                ("DCVMIN", Float64), ("CCCITMX", Int64), ("CCCACC", Float64),
                                ("IPR", Int64), ("NBR", Int64), ("ANMXR", Float64),
                                ("ANMNR", Float64), ("RCR", Float64), ("XCR", Float64),
                                ("EBASR", Float64), ("TRR", Float64), ("TAPR", Float64),
                                ("TMXR", Float64), ("TMNR", Float64), ("STPR", Float64),
                                ("ICR", Int64), ("IFR", Int64), ("ITR", Int64), ("IDR", String),
                                ("XCAPR", Float64), ("IPI", Int64), ("NBI", Int64),
                                ("ANMXI", Float64), ("ANMNI", Float64), ("RCI", Float64),
                                ("XCI", Float64), ("EBASI", Float64), ("TRI", Float64),
                                ("TAPI", Float64), ("TMXI", Float64), ("TMNI", Float64),
                                ("STPI", Float64), ("ICI", Int64), ("IFI", Int64),
                                ("ITI", Int64), ("IDI", String), ("XCAPI", Float64)]

    vsc_line_dtypes = [("NAME", String), ("MDC", Int64), ("RDC", Float64),
                       ("O1", Int64), ("F1", Float64),
                       ("O2", Int64), ("F2", Float64),
                       ("O3", Int64), ("F3", Float64),
                       ("O4", Int64), ("F4", Float64)]

    vsc_subline_dtypes = [("IBUS", Int64), ("TYPE", Int64), ("MODE", Int64),
                          ("DCSET", Float64), ("ACSET", Float64),
                          ("ALOSS", Float64), ("BLOSS", Float64), ("MINLOSS", Float64),
                          ("SMAX", Float64), ("IMAX", Float64), ("PWF", Float64),
                          ("MAXQ", Float64), ("MINQ", Float64),
                          ("REMOT", Int64), ("RMPCT", Float64)]

    impedance_correction_dtypes = [("I", Int64), ("T1", Float64), ("F1", Float64),
                                   ("T2", Float64), ("F2", Float64), ("T3", Float64), ("F3", Float64),
                                   ("T4", Float64), ("F4", Float64), ("T5", Float64), ("F5", Float64),
                                   ("T6", Float64), ("F6", Float64), ("T7", Float64), ("F7", Float64),
                                   ("T8", Float64), ("F8", Float64), ("T9", Float64), ("F9", Float64),
                                   ("T10", Float64), ("F10", Float64), ("T11", Float64), ("F11", Float64)]

    multi_term_main_dtypes = [("NAME", String), ("NCONV", Int64), ("NDCBS", Int64), ("NDCLN", Int64),
                              ("MDC", Int64), ("VCONV", Int64), ("VCMOD", Float64), ("VCONVN", Float64)]

    multi_term_nconv_dtypes = [("IB", Int64), ("N", Int64), ("ANGMX", Float64), ("ANGMN", Float64),
                               ("RC", Float64), ("XC", Float64), ("EBAS", Float64), ("TR", Float64),
                               ("TAP", Float64), ("TPMX", Float64), ("TPMN", Float64), ("TSTP", Float64),
                               ("SETVL", Float64), ("DCPF", Float64), ("MARG", Float64), ("CNVCOD", Int64)]

    multi_term_ndcbs_dtypes = [("IDC", Int64), ("IB", Int64), ("AREA", Int64), ("ZONE", Int64),
                               ("DCNAME", String), ("IDC2", Int64), ("RGRND", Float64), ("OWNER", Int64)]

    multi_term_ndcln_dtypes = [("IDC", Int64), ("JDC", Int64), ("DCCKT", String), ("MET", Int64),
                               ("RDC", Float64), ("LDC", Float64)]

    multi_section_dtypes = [("I", Int64), ("J", Int64), ("ID", String), ("MET", Int64), ("DUM1", Int64),
                            ("DUM2", Int64), ("DUM3", Int64), ("DUM4", Int64), ("DUM5", Int64),
                            ("DUM6", Int64), ("DUM7", Int64), ("DUM8", Int64), ("DUM9", Int64)]

    zone_dtypes = [("I", Int64), ("ZONAME", String)]

    interarea_dtypes = [("ARFROM", Int64), ("ARTO", Int64), ("TRID", String), ("PTRAN", Float64)]

    owner_dtypes = [("I", Int64), ("OWNAME", String)]

    FACTS_dtypes = [("NAME", String), ("I", Int64), ("J", Int64), ("MODE", Int64), ("PDES", Float64),
                    ("QDES", Float64), ("VSET", Float64), ("SHMX", Float64), ("TRMX", Float64),
                    ("VTMN", Float64), ("VTMX", Float64), ("VSMX", Float64), ("IMX", Float64),
                    ("LINX", Float64), ("RMPCT", Float64), ("OWNER", Int64), ("SET1", Float64),
                    ("SET2", Float64), ("VSREF", Int64), ("REMOT", Int64), ("MNAME", String)]

    switched_shunt_dtypes = [("I", Int64), ("MODSW", Int64), ("ADJM", Int64), ("STAT", Int64),
                             ("VSWHI", Float64), ("VSWLO", Float64), ("SWREM", Int64),
                             ("RMPCT", Float64), ("RMIDNT", String), ("BINIT", Float64),
                             ("N1", Int64), ("B1", Float64), ("N2", Int64), ("B2", Float64),
                             ("N3", Int64), ("B3", Float64), ("N4", Int64), ("B4", Float64),
                             ("N5", Int64), ("B5", Float64), ("N6", Int64), ("B6", Float64),
                             ("N7", Int64), ("B7", Float64), ("N8", Int64), ("B8", Float64)]

    # TODO: Account for multiple lines in GNE Device entries
    gne_device_dtypes = [("NAME", String), ("MODEL", String), ("NTERM", Int64), ("BUSi", Int64),
                         ("NREAL", Int64), ("NINTG", Int64), ("NCHAR", Int64),

                         ("STATUS", Int64), ("OWNER", Int64), ("NMETR", Int64),

                         ("REALi", Float64),
                         ("INTGi", Int64),
                         ("CHARi", String)]

    induction_machine_dtypes = [("I", Int64), ("ID", String), ("STAT", Int64), ("SCODE", Int64),
                                ("DCODE", Int64), ("AREA", Int64), ("ZONE", Int64), ("OWNER", Int64),
                                ("TCODE", Int64), ("BCODE", Int64), ("MBASE", Float64),
                                ("RATEKV", Float64), ("PCODE", Int64), ("PSET", Float64),
                                ("H", Float64), ("A", Float64), ("B", Float64), ("D", Float64),
                                ("E", Float64), ("RA", Float64), ("XA", Float64), ("XM", Float64),
                                ("R1", Float64), ("X1", Float64), ("R2", Float64), ("X2", Float64),
                                ("X3", Float64), ("E1", Float64), ("SE1", Float64), ("E2", Float64),
                                ("SE2", Float64), ("IA1", Float64), ("IA2", Float64), ("XAMULT", Float64)]

                                dtypes = Dict{String,Array}("BUS" => bus_dtypes,
                                "LOAD" => load_dtypes,
                                "FIXED SHUNT" => fixed_shunt_dtypes,
                                "GENERATOR" => generator_dtypes,
                                "BRANCH" => branch_dtypes,
                                "TRANSFORMER" => transformer_dtypes,
                                "TRANSFORMER TWO-WINDING LINE 1" => transformer_2_1_dtypes,
                                "TRANSFORMER TWO-WINDING LINE 2" => transformer_2_2_dtypes,
                                "TRANSFORMER TWO-WINDING LINE 3" => transformer_2_3_dtypes,
                                "TRANSFORMER THREE-WINDING LINE 1" => transformer_3_1_dtypes,
                                "TRANSFORMER THREE-WINDING LINE 2" => transformer_3_2_dtypes,
                                "TRANSFORMER THREE-WINDING LINE 3" => transformer_3_3_dtypes,
                                "TRANSFORMER THREE-WINDING LINE 4" => transformer_3_4_dtypes,
                                "AREA INTERCHANGE" => area_interchange_dtypes,
                                "TWO-TERMINAL DC" => two_terminal_line_dtypes,
                                "VOLTAGE SOURCE CONVERTER" => vsc_line_dtypes,
                                "VOLTAGE SOURCE CONVERTER SUBLINES" => vsc_subline_dtypes,
                                "IMPEDANCE CORRECTION" => impedance_correction_dtypes,
                                "MULTI-TERMINAL DC" => multi_term_main_dtypes,
                                "MULTI-TERMINAL DC NCONV" => multi_term_nconv_dtypes,
                                "MULTI-TERMINAL DC NDCBS" => multi_term_ndcbs_dtypes,
                                "MULTI-TERMINAL DC NDCLN" => multi_term_ndcln_dtypes,
                                "MULTI-SECTION LINE" => multi_section_dtypes,
                                "ZONE" => zone_dtypes,
                                "INTER-AREA TRANSFER" => interarea_dtypes,
                                "OWNER" => owner_dtypes,
                                "FACTS CONTROL DEVICE" => FACTS_dtypes,
                                "SWITCHED SHUNT" => switched_shunt_dtypes,
                                "CASE IDENTIFICATION" => transaction_dtypes,
                                "GNE DEVICE" => gne_device_dtypes,
                                "INDUCTION MACHINE" => induction_machine_dtypes)

    return dtypes[field_name]
end

"""
    _get_component_property(section, ret, search_field, search_value)

Internal function. Finds a component in `section` where `search_field` ==
`search_value` and returns `ret` from that component.
"""
function _get_component_property(section, ret, search_field, search_value)
    for component in section
        if component[search_field] == search_value
            return component[ret]
        end
    end
    Memento.warn(_LOGGER, "cannot find $search_field = $search_value, no return for $ret")
end

"""
    _get_pti_default(section, field, data, component)

Internal function. Returns a default value in `section` for `field` in
`component` from `data`.
"""
function _get_pti_default(section::AbstractString, field::AbstractString, data::Dict, component::Dict; sub_field=nothing)
    case_identification = Dict{String,Any}("IC" => 0, "SBASE" => 100.0, "REV" => 33,
                                           "XFRRAT" => 0, "NXFRAT" => 0, "BASFRQ" => 60)

    bus = Dict{String,Any}("BASKV" => 0.0, "IDE" => 1, "AREA" => 1, "ZONE" => 1,
                           "OWNER" => 1, "VM" => 1.0, "VA" => 0.0, "NVHI" => 1.1,
                           "NVLO" => 0.9, "EVHI" => 1.1, "EVLO" => 0.9,
                           "NAME" => "            ")

    load = Dict{String,Any}("ID" => 1, "STATUS" => 1, "PL" => 0.0, "QL" => 0.0,
                            "IP" => 0.0, "IQ" => 0.0, "YP" => 0.0, "YQ" => 0.0,
                            "SCALE" => 1, "INTRPT" => 0,
                            "AREA" => Expr(:call, :_get_component_property, data["BUS"], "AREA", "I", get(component, "I", 0)),
                            "ZONE" => Expr(:call, :_get_component_property, data["BUS"], "ZONE", "I", get(component, "I", 0)),
                            "OWNER" => Expr(:call, :_get_component_property, data["BUS"], "OWNER", "I", get(component, "I", 0)))

    fixed_shunt = Dict{String,Any}("ID" => 1, "STATUS" => 1, "GL" => 0.0, "BL" => 0.0)

    generator = Dict{String,Any}("ID" => 1, "PG" => 0.0, "QG" => 0.0, "QT" => 9999.0,
                                 "QB" => -9999.0, "VS" => 1.0, "IREG" => 0,
                                 "MBASE" => data["CASE IDENTIFICATION"][1]["SBASE"],
                                 "ZR" => 0.0, "ZX" => 1.0, "RT" => 0.0, "XT" => 0.0,
                                 "GTAP" => 1.0, "STAT" => 1, "RMPCT" => 100.0,
                                 "PT" => 9999.0, "PB" => -9999.0,
                                 "O1" => Expr(:call, :_get_component_property, data["BUS"], "OWNER", "I", get(component, "I", 0)),
                                 "O2" => 0, "O3" => 0, "O4" => 0, "F1" => 1.0,"F2" => 1.0, "F3" => 1.0,
                                 "F4" => 1.0, "WMOD" => 0, "WPF" => 1.0)

    branch = Dict{String,Any}("CKT" => 1, "B" => 0.0, "RATEA" => 0.0, "RATEB" => 0.0,
                              "RATEC" => 0.0, "GI" => 0.0, "BI" => 0.0, "GJ" => 0.0,
                              "BJ" => 0.0, "ST" => 1, "MET" => 1, "LEN" => 0.0,
                              "O1" => Expr(:call, :_get_component_property, data["BUS"], "OWNER", "I", get(component, "I", 0)),
                              "O2" => 0, "O3" => 0, "O4" => 0, "F1" => 1.0,
                              "F2" => 1.0, "F3" => 1.0, "F4" => 1.0)

    transformer = Dict{String,Any}("K" => 0, "CKT" => 1, "CW" => 1, "CZ" => 1, "CM" => 1,
                                   "MAG1" => 0.0, "MAG2" => 0.0, "NMETR" => 2,
                                   "NAME" => "            ", "STAT" => 1,
                                   "O1" => Expr(:call, :_get_component_property, data["BUS"], "OWNER", "I", get(component, "I", 0)),
                                   "O2" => 0, "O3" => 0, "O4" => 0,
                                   "F1" => 1.0, "F2" => 1.0, "F3" => 1.0, "F4" => 1.0,
                                   "VECGRP" => "            ", "R1-2" => 0.0,
                                   "SBASE1-2" => data["CASE IDENTIFICATION"][1]["SBASE"], "R2-3" => 0.0,
                                   "SBASE2-3" => data["CASE IDENTIFICATION"][1]["SBASE"], "R3-1" => 0.0,
                                   "SBASE3-1" => data["CASE IDENTIFICATION"][1]["SBASE"], "VMSTAR" => 1.0,
                                   "ANSTAR" => 0.0,
                                   "WINDV1" => get(component, "CW", 1) == 2 ? Expr(:call, :_get_component_property, data["BUS"], "BASKV", "I", get(component, "I", 0)) : 1.0,
                                   "NOMV1" => 0.0, "ANG1" => 0.0, "RATA1" => 0.0,
                                   "RATB1" => 0.0, "RATC1" => 0.0, "COD1" => 0,
                                   "CONT1" => 0, "RMA1" => 1.1, "RMI1" => 0.9,
                                   "VMA1" => 1.1, "VMI1" => 0.9, "NTP1" => 33,
                                   "TAB1" => 0, "CR1" => 0.0, "CX1" => 0.0, "CNXA1" => 0.0,
                                   "WINDV2" => get(component, "CW", 1) == 2 ? Expr(:call, :_get_component_property, data["BUS"], "BASKV", "I", get(component, "I", 0)) : 1.0,
                                   "NOMV2" => 0.0, "ANG2" => 0.0, "RATA2" => 0.0,
                                   "RATB2" => 0.0, "RATC2" => 0.0, "COD2" => 0,
                                   "CONT2" => 0, "RMA2" => 1.1, "RMI2" => 0.9,
                                   "VMA2" => 1.1, "VMI2" => 0.9, "NTP2" => 33,
                                   "TAB2" => 0, "CR2" => 0.0, "CX2" => 0.0,
                                   "CNXA2" => 0.0,
                                   "WINDV3" => get(component, "CW", 1) == 2 ? Expr(:call, :_get_component_property, data["BUS"], "BASKV", "I", get(component, "I", 0)) : 1.0,
                                   "NOMV3" => 0.0, "ANG3" => 0.0, "RATA3" => 0.0,
                                   "RATB3" => 0.0, "RATC3" => 0.0, "COD3" => 0,
                                   "CONT3" => 0, "RMA3" => 1.1, "RMI3" => 0.9,
                                   "VMA3" => 1.1, "VMI3" => 0.9, "NTP3" => 33,
                                   "TAB3" => 0, "CR3" => 0.0, "CX3" => 0.0,
                                   "CNXA3" => 0.0)

    area_interchange = Dict{String,Any}("ISW" => 0, "PDES" => 0.0, "PTOL" => 10.0,
                                        "ARNAME" => "            ")

    two_terminal_dc = Dict{String,Any}("MDC" => 0, "VCMOD" => 0.0, "RCOMP" => 0.0,
                                       "DELTI" => 0.0, "METER" => "I", "DCVMIN" => 0.0,
                                       "CCCITMX" => 20, "CCCACC" => 1.0, "TRR" => 1.0,
                                       "TAPR" => 1.0, "TMXR" => 1.5, "TMNR" => 0.51,
                                       "STPR" => 0.00625, "ICR" => 0, "IFR" => 0,
                                       "ITR" => 0, "IDR" => "1", "XCAPR" => 0.0,
                                       "TRI" => 1.0,
                                       "TAPI" => 1.0, "TMXI" => 1.5, "TMNI" => 0.51,
                                       "STPI" => 0.00625, "ICI" => 0, "IFI" => 0,
                                       "ITI" => 0, "IDI" => "1", "XCAPI" => 0.0)

    vsc_dc = Dict{String,Any}("MDC" => 1,
                              "O1" => Expr(:call, :_get_component_property, data["BUS"], "OWNER", "I", get(get(component, "CONVERTER BUSES", [Dict()])[1], "IBUS", 0)),
                              "O2" => 0, "O3" => 0, "O4" => 0,
                              "F1" => 1.0, "F2" => 1.0, "F3" => 1.0, "F4" => 1.0,
                              "CONVERTER BUSES" => Dict{String,Any}("MODE" => 1, "ACSET" => 1.0, "ALOSS" => 1.0, "BLOSS" => 0.0, "MINLOSS" => 0.0,
                                                                    "SMAX" => 0.0, "IMAX" => 0.0, "PWF" => 1.0, "MAXQ" => 9999.0, "MINQ" => -9999.0,
                                                                    "REMOT" => 0, "RMPCT" => 100.0))

    impedance_correction = Dict{String,Any}("T1" => 0.0, "T2" => 0.0, "T3" => 0.0,
                                            "T4" => 0.0, "T5" => 0.0, "T6" => 0.0,
                                            "T7" => 0.0, "T8" => 0.0, "T9" => 0.0,
                                            "T10" => 0.0, "T11" => 0.0, "F1" => 0.0,
                                            "F2" => 0.0, "F3" => 0.0, "F4" => 0.0,
                                            "F5" => 0.0, "F6" => 0.0, "F7" => 0.0,
                                            "F8" => 0.0, "F9" => 0.0, "F10" => 0.0,
                                            "F11" => 0.0)

    multi_term_dc = Dict{String,Any}("MDC" => 0, "VCMOD" => 0.0, "VCONVN" => 0,
                                     "CONV" => Dict{String,Any}("TR" => 1.0, "TAP" => 1.0, "TPMX" => 1.5, "TPMN" => 0.51, "TSTP" => 0.00625,
                                                                "DCPF" => 1, "MARG" => 0.0, "CNVCOD" => 1),
                                     "DCBS" => Dict{String,Any}("IB" => 0.0, "AREA" => 1, "ZONE" => 1, "DCNAME" => "            ", "IDC2" => 0,
                                                                "RGRND" => 0.0, "OWNER" => 1),
                                     "DCLN" => Dict{String,Any}("DCCKT" => 1, "MET" => 1, "LDC" => 0.0))

    multi_section = Dict{String,Any}("ID" => "&1", "MET" => 1)

    zone = Dict{String,Any}("ZONAME" => "            ")

    interarea = Dict{String,Any}("TRID" => 1, "PTRAN" => 0.0)

    owner = Dict{String,Any}("OWNAME" => "            ")

    facts = Dict{String,Any}("J" => 0, "MODE" => 1, "PDES" => 0.0, "QDES" => 0.0,
                             "VSET" => 1.0, "SHMX" => 9999.0, "TRMX" => 9999.0,
                             "VTMN" => 0.9, "VTMX" => 1.1, "VSMX" => 1.0,
                             "IMX" => 0.0, "LINX" => 0.05, "RMPCT" => 100.0,
                             "OWNER" => 1, "SET1" => 0.0, "SET2" => 0.0,
                             "VSREF" => 0, "REMOT" => 0, "MNAME" => "")

    switched_shunt = Dict{String,Any}("MODSW" => 1, "ADJM" => 0, "STAT" => 1,
                                      "VSWHI" => 1.0, "VSWLO" => 1.0, "SWREM" => 0,
                                      "RMPCT" => 100.0, "RMIDNT" => "", "BINIT" => 0.0,
                                      "N1" => 0.0, "N2" => 0.0, "N3" => 0.0,
                                      "N4" => 0.0, "N5" => 0.0, "N6" => 0.0,
                                      "N7" => 0.0, "N8" => 0.0, "B1" => 0.0,
                                      "B2" => 0.0, "B3" => 0.0, "B4" => 0.0,
                                      "B5" => 0.0, "B6" => 0.0, "B7" => 0.0,
                                      "B8" => 0.0)

    gne_device = Dict{String,Any}("NTERM" => 1, "NREAL" => 0, "NINTG" => 0,
                                  "NCHAR" => 0, "STATUS" => 1,
                                  "OWNER" => Expr(:call, :_get_component_property, data["BUS"], "OWNER", "I", get(component, "BUS1", 0)),
                                  "NMETR" => get(component, "NTERM", 1), "REAL" => 0,
                                  "INTG" => nothing,
                                  "CHAR" => "1")

    induction_machine = Dict{String,Any}("ID" => 1, "STAT" => 1, "SCODE" => 1, "DCODE" => 2,
                                         "AREA" => Expr(:call, :_get_component_property, data["BUS"], "AREA", "I", get(component, "I", 0)),
                                         "ZONE" => Expr(:call, :_get_component_property, data["BUS"], "ZONE", "I", get(component, "I", 0)),
                                         "OWNER" => Expr(:call, :_get_component_property, data["BUS"], "OWNER", "I", get(component, "I", 0)),
                                         "TCODE" => 1, "BCODE" => 1,
                                         "MBASE" => data["CASE IDENTIFICATION"][1]["SBASE"], "RATEKV" => 0.0,
                                         "PCODE" => 1, "H" => 1.0, "A" => 1.0,
                                         "B" => 1.0, "D" => 1.0, "E" => 1.0,
                                         "RA" => 0.0, "XA" => 0.0, "XM" => 2.5,
                                         "R1" => 999.0, "X1" => 999.0, "R2" => 999.0,
                                         "X2" => 999.0, "X3" => 0.0, "E1" => 1.0,
                                         "SE1" => 0.0, "E2" => 1.2, "SE2" => 0.0,
                                         "IA1" => 0.0, "IA2" => 0.0, "XAMULT" => 1)

    defaults = Dict{String,Dict}("BUS" => bus,
                                  "LOAD" => load,
                                  "FIXED SHUNT" => fixed_shunt,
                                  "GENERATOR" => generator,
                                  "BRANCH" => branch,
                                  "TRANSFORMER" => transformer,
                                  "AREA INTERCHANGE" => area_interchange,
                                  "TWO-TERMINAL DC" => two_terminal_dc,
                                  "VOLTAGE SOURCE CONVERTER" => vsc_dc,
                                  "IMPEDANCE CORRECTION" => impedance_correction,
                                  "MULTI-TERMINAL DC" => multi_term_dc,
                                  "MULTI-SECTION LINE" => multi_section,
                                  "ZONE" => zone,
                                  "INTER-AREA TRANSFER" => interarea,
                                  "OWNER" => owner,
                                  "FACTS CONTROL DEVICE" => facts,
                                  "SWITCHED SHUNT" => switched_shunt,
                                  "CASE IDENTIFICATION" => case_identification,
                                  "GNE DEVICE" => gne_device,
                                  "INDUCTION MACHINE" => induction_machine)

    if sub_field != nothing
        return eval(defaults[section][field][sub_field])
    else
        return eval(defaults[section][field])
    end
end

"""
_parse_line_element!(data, elements, section)

Internal function. Parses a single "line" of data elements from a PTI file, as
given by `elements` which is an array of the line, typically split at `,`.
Elements are parsed into data types given by `section` and saved into `data::Dict`.
"""
function _parse_line_element!(data::Dict, elements::Array, section::AbstractString)
    missing_fields = []
    for (field, dtype) in _get_pti_dtypes(section)
        try
            element = popfirst!(elements)
        catch message
            if isa(message, ArgumentError)
                @debug "Have run out of elements in $section at $field"
                push!(missing_fields, field)
                continue
            end
        end

        if startswith(strip(element), "'") && endswith(strip(element), "'")
            dtype = String
            element = chop(reverse(chop(reverse(strip(element)))))
        end

        try
            if dtype != String && element != ""
                data[field] = parse(dtype, element)
            else
                data[field] = element
            end
        catch message
            if isa(message, Meta.ParseError)
                data[field] = element
            else
                error("value '$element' for $field in section $section is not of type $dtype.")
            end
        end
    end

    if length(missing_fields) > 0
        for field in missing_fields
            data[field] = ""
        end
        missing_str = join(missing_fields, ", ")
        if !(section == "SWITCHED SHUNT" && startswith(missing_str, "N")) &&
            !(section == "MULTI-SECTION LINE" && startswith(missing_str, "DUM")) &&
            !(section == "IMPEDANCE CORRECTION" && startswith(missing_str, "T"))
            @info("The following fields in $section are missing: $missing_str")
        end
    end
end

"""
    add_section_data!(pti_data, section_data, section)

Internal function. Adds `section_data::Dict`, which contains all parsed
elements of a PTI file section given by `section`, into the parent
`pti_data::Dict`
"""
function _add_section_data!(pti_data::Dict, section_data::Dict, section::AbstractString)
    try
        pti_data[section] = append!(pti_data[section], [deepcopy(section_data)])
    catch message
        if isa(message, KeyError)
            pti_data[section] = [deepcopy(section_data)]
        else
            @error(sprint(showerror, message))
        end
    end
end

"""
    _get_line_elements(line)

Internal function. Uses regular expressions to extract all separate data
elements from a line of a PTI file and populate them into an `Array{String}`.
Comments, typically indicated at the end of a line with a `'/'` character,
are also extracted separately, and `Array{Array{String}, String}` is returned.
"""
function _get_line_elements(line::AbstractString)::Array
    if length(collect(eachmatch(r"'", line))) % 2 == 1
        throw(DataFormatError("There are an uneven number of single-quotes in \"{line}\", the line cannot be parsed."))
    end

    comment_split = r"(?!\B[\'][^\']*)[\/](?![^\']*[\']\B)"
    line_comment = split(line, comment_split, limit=2)
    line = strip(line_comment[1])
    comment = length(line_comment) > 1 ? strip(line_comment[2]) : ""

    split_string = r",(?=(?:[^']*'[^']*')*[^']*$)"
    elements = [strip(element) for element in split(line, split_string)]

    @debug("$line")
    @debug("$comment")
    @debug("$elements")

    return [elements, comment]
end

"""
    _parse_pti_data(data_string, sections)

Internal function. Parse a PTI raw file into a `Dict`, given the
`data_string` of the file and a list of the `sections` in the PTI
file (typically given by default by `get_pti_sections()`.
"""
function _parse_pti_data(data_io::IO, sections::Array)
    data_lines = readlines(data_io)
    skip_lines = 0
    skip_sublines = 0
    subsection = ""

    pti_data = Dict{String,Array{Dict}}()

    section = popfirst!(sections)
    section_data = Dict{String,Any}()

    for (line_number, line) in enumerate(data_lines)
        @debug "$line_number: $line"

        (elements, comment) = _get_line_elements(line)

        if length(elements) != 0 && elements[1] == "Q" && line_number > 3
            break

        elseif length(elements) != 0 && elements[1] == "0" && line_number > 3
            if line_number == 4
                section = popfirst!(sections)
            end

            if length(elements) > 1
                @info("At line $line_number, new section started with '0', but additional non-comment data is present. Pattern '^\\s*0\\s*[/]*.*' is reserved for section start/end.")
            elseif length(comment) > 0
                @info("At line $line_number, unexpected section: expected: $section, comment specified: $(section)")
            end

            if !isempty(sections)
                section = popfirst!(sections)
            end

            continue
        else
            if line_number == 4
                section = popfirst!(sections)
                section_data = Dict{String,Any}()
            end

            if skip_lines > 0
                skip_lines -= 1
                continue
            end

            @debug join(["Section:", section], " ")
            if !(section in ["CASE IDENTIFICATION","TRANSFORMER","VOLTAGE SOURCE CONVERTER","MULTI-TERMINAL DC","TWO-TERMINAL DC","GNE DEVICE"])
                section_data = Dict{String,Any}()
                try
                    _parse_line_element!(section_data, elements, section)
                catch message
                    throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                end

            elseif section == "CASE IDENTIFICATION"
                if line_number == 1
                    try
                        _parse_line_element!(section_data, elements, section)
                    catch message
                        throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                    end

                    if section_data["REV"] != "" && section_data["REV"] < 33
                        @info("Version $(section_data["REV"]) of PTI format is unsupported, parser may not function correctly.")
                    end
                else
                    section_data["Comment_Line_$(line_number - 1)"] = line
                end

                if line_number < 3
                    continue
                end

            elseif section == "TRANSFORMER"
                section_data = Dict{String,Any}()
                if parse(Int64, _get_line_elements(line)[1][3]) == 0 # two winding transformer
                    winding = "TWO-WINDING"
                    skip_lines = 3
                elseif parse(Int64, _get_line_elements(line)[1][3]) != 0 # three winding transformer
                    winding = "THREE-WINDING"
                    skip_lines = 4
                else
                    @error("Cannot detect type of Transformer")
                end

                try
                    for transformer_line in 0:4
                        if transformer_line == 0
                            temp_section = section
                        else
                            temp_section = join([section, winding, "LINE", transformer_line], " ")
                        end

                        if winding == "TWO-WINDING" && transformer_line == 4
                            break
                        else
                            elements = _get_line_elements(data_lines[line_number + transformer_line])[1]
                            _parse_line_element!(section_data, elements, temp_section)
                        end
                    end
                catch message
                    throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                end

            elseif section == "VOLTAGE SOURCE CONVERTER"
                if length(_get_line_elements(line)[1]) == 11
                    section_data = Dict{String,Any}()
                    try
                        _parse_line_element!(section_data, elements, section)
                    catch message
                        throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                    end
                    skip_sublines = 2
                    continue

                elseif skip_sublines > 0
                    skip_sublines -= 1

                    (elements, comment) = _get_line_elements(line)
                    subsection_data = Dict{String,Any}()

                    for (field, dtype) in _get_pti_dtypes("$section SUBLINES")
                        element = popfirst!(elements)
                        if element != ""
                            subsection_data[field] = parse(dtype, element)
                        else
                            subsection_data[field] = ""
                        end
                    end

                    try
                        section_data["CONVERTER BUSES"] = append!(section_data["CONVERTER BUSES"], [deepcopy(subsection_data)])
                    catch message
                        if isa(message, KeyError)
                            section_data["CONVERTER BUSES"] = [deepcopy(subsection_data)]
                            continue
                        else
                            @error(message)
                        end
                    end
                end

            elseif section == "TWO-TERMINAL DC"
                section_data = Dict{String,Any}()
                if length(_get_line_elements(line)[1]) == 12
                    (elements, comment) = _get_line_elements(join(data_lines[line_number:line_number + 2], ','))
                    skip_lines = 2
                end

                try
                    _parse_line_element!(section_data, elements, section)
                catch message
                    throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                end

            elseif section == "MULTI-TERMINAL DC"
                if skip_sublines == 0
                    section_data = Dict{String,Any}()
                    try
                        _parse_line_element!(section_data, elements, section)
                    catch message
                        throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                    end

                    if section_data["NCONV"] > 0
                        skip_sublines = section_data["NCONV"]
                        subsection = "NCONV"
                        continue
                    elseif section_data["NDCBS"] > 0
                        skip_sublines = section_data["NDCBS"]
                        subsection = "NDCBS"
                        continue
                    elseif section_data["NDCLN"] > 0
                        skip_sublines = section_data["NDCLN"]
                        subsection = "NDCLN"
                        continue
                    end
                end

                if skip_sublines > 0
                    skip_sublines -= 1

                    subsection_data = Dict{String,Any}()

                    try
                        _parse_line_element!(subsection_data, elements, "$section $subsection")
                    catch message
                        throw(error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                    end

                    try
                        section_data["$(subsection[2:end])"] = push!(section_data["$(subsection[2:end])"], deepcopy(subsection_data))
                        if skip_sublines > 0 && subsection != "NDCLN"
                            continue
                        end
                    catch message
                        if isa(message, KeyError)
                            section_data["$(subsection[2:end])"] = [deepcopy(subsection_data)]
                            if skip_sublines > 0 && subsection != "NDCLN"
                                continue
                            end
                        else
                            @error(sprint(showerror, message))
                        end
                    end

                    if skip_sublines == 0 && subsection != "NDCLN"
                        if subsection == "NDCBS"
                            skip_sublines = section_data["NDCLN"]
                            subsection = "NDCLN"
                            continue
                        elseif subsection == "NCONV"
                            skip_sublines = section_data["NDCBS"]
                            subsection = "NDCBS"
                            continue
                        end
                    elseif skip_sublines == 0 && subsection == "NDCLN"
                        subsection = ""
                    else
                        continue
                    end
                end

            elseif section == "GNE DEVICE"
                # TODO: handle multiple lines of GNE Device
                @info("GNE DEVICE parsing is not supported.")
            end
        end
        if subsection != ""
            @debug "appending data"
        end
        _add_section_data!(pti_data, section_data, section)
    end

    _populate_defaults!(pti_data)

    return pti_data
end

"""
    parse_pti(filename::String)

Open PTI raw file given by `filename`, returning a `Dict` of the data parsed
into the proper types.
"""
function parse_pti(filename::String)::Dict
    pti_data = open(filename) do f
        parse_pti(f)
    end

    return pti_data
end

"""
    parse_pti(io::IO)

Reads PTI data in `io::IO`, returning a `Dict` of the data parsed into the
proper types.
"""
function parse_pti(io::IO)::Dict
    pti_data = _parse_pti_data(io, _get_pti_sections())
    try
        pti_data["CASE IDENTIFICATION"][1]["NAME"] = match(r"^\<file\s[\/\\]*(?:.*[\/\\])*(.*)\.raw\>$", lowercase(io.name)).captures[1]
    catch
        throw(error( "This file is unrecognized and cannot be parsed"))
    end

    return pti_data
end

"""
    _populate_defaults!(pti_data)

Internal function. Populates empty fields with PSS(R)E PTI v33 default values
"""
function _populate_defaults!(data::Dict)
    for section in _get_pti_sections()
        if haskey(data, section)
            section_components = []
            for component in data[section]
                new_component = deepcopy(component)
                for (field, field_value) in component
                    if isa(field_value, Array)
                        new_array = []
                        for sub_component in field_value
                            new_sub_component = deepcopy(sub_component)
                            for (sub_field, sub_field_value) in sub_component
                                if sub_field_value == ""
                                    try
                                        new_sub_component[sub_field] = _get_pti_default(section, field, data, sub_component; sub_field=sub_field)
                                    catch msg
                                        if isa(msg, KeyError)
                                            @warn( "'$sub_field' in '$field' in '$section' has no default value")
                                        else
                                            rethrow(msg)
                                        end
                                    end
                                end
                            end
                            push!(new_array, new_sub_component)
                        end
                        new_component[field] = new_array
                    elseif field_value == "" && !(field in ["Comment_Line_1", "Comment_Line_2"]) && !startswith(field, "DUM")
                        try
                            new_component[field] = _get_pti_default(section, field, data, component)
                        catch msg
                            if isa(msg, KeyError)
                                @warn("'$field' in '$section' has no default value")
                            else
                                rethrow(msg)
                            end
                        end
                    end
                end
                push!(section_components, new_component)
            end
        data[section] = section_components
        end
    end
end


text:
#=
    AbstractGUI
      AbsractContext, AbstractContent
         AbstractConstrain, AbstractConstruct, AbstractContemper

     AbstractGUI
       AbstractContent
          AbstractConstrainContent
          AbstractConstructContent
          AbstractContemperContent
       AbstractContext
          AbstractConstrainContext
          AbstractConstructContext
          AbstractContemperContext
=#

"""
    AbstractGUI

The root for type abstractions involved in realizing a GUI.
"""
abstract type AbstractGUI end

"""
    AbstractContext

The abstraction for types that provide, characterize,
organize, or reflect context involved in realizing a GUI.
"""
abstract type AbstractContext  <: AbstractGUI end

"""
    AbstractContent

The abstraction for types that develop, bring, carry,
or hold content involved in realizing a GUI.
"""
abstract type AbstractContent <: AbstractGUI end

"""
    AbstractConstrain

The abstraction for types constraining aspects of elements
involved in realizing a GUI.
"""
abstract type AbstractConstrain <: AbstractGUI end

"""
    AbstractConstrainContext

The abstraction for types constraining aspects of elements
involved in realizing the Context of a GUI.
"""
abstract type AbstractConstrainContext <: AbstractConstrain end

"""
    AbstractConstrainContent

The abstraction for types constraining aspects of elements
involved in realizing the Content of a GUI.
"""
abstract type AbstractConstrainContent <: AbstractConstrain end

"""
    AbstractConstruct

The abstraction for types constructing aspects of elements
involved in realizing a GUI.
"""
abstract type AbstractConstruct <: AbstractGUI end

"""
    AbstractConstructContext

The abstraction for types constructing aspects of elements
involved in realizing the Context of a GUI.
"""
abstract type AbstractConstructContext <: AbstractConstruct end

"""
    AbstractConstructContent

The abstraction for types constructing aspects of elements
involved in realizing the Content of a GUI.
"""
abstract type AbstractConstructContent <: AbstractConstruct end

"""
    AbstractContemper

The abstraction for types contempering aspects of elements
involved in realizing a GUI.

_contemper_ "to gain pliable strength through tempering"
- [_with_][_tempering_]
- tempering through measured admixture of beneficials
"""
abstract type AbstractContemper <: AbstractGUI end

"""
    AbstractContemperContext

The abstraction for types constructing aspects of elements
involved in realizing the Context of a GUI.
"""
abstract type AbstractContemperContext <: AbstractContemper end

"""
    AbstractContemperContent

The abstraction for types constructing aspects of elements
involved in realizing the Content of a GUI.
"""
abstract type AbstractContemperContent <: AbstractContemper end


text:
# utils

macro embed(typename)
    t = getfield(@__MODULE__, typename)
    e = Expr(:block, [])
    for (x, y) in zip(fieldnames(t), fieldtypes(t))
        push!(e.args, :($x::$y))
    end
    return esc(e)
end

module CcallMacro

function parsecall(expr::Expr)
    # setup and check for errors
    if !Meta.isexpr(expr, :(::))
        throw(ArgumentError("@ccall needs a function signature with a return type"))
    end
    rettype = expr.args[2]

    call = expr.args[1]
    if !Meta.isexpr(call, :call)
        throw(ArgumentError("@ccall has to take a function call"))
    end

    # get the function symbols
    func = let f = call.args[1]
        f isa Expr ? :(($(f.args[2]), $(f.args[1]))) : QuoteNode(f)
    end

    # detect varargs
    varargs = nothing
    argstart = 2
    callargs = call.args
    if length(callargs) >= 2 && Meta.isexpr(callargs[2], :parameters)
        argstart = 3
        varargs = callargs[2].args
    end

    # collect args and types
    args = []
    types = []

    function pusharg!(arg)
        if !Meta.isexpr(arg, :(::))
            throw(ArgumentError("args in @ccall need type annotations. '$(repr(arg))' doesn't have one."))
        end
        push!(args, arg.args[1])
        push!(types, arg.args[2])
    end

    for i in argstart:length(callargs)
        pusharg!(callargs[i])
    end
    # add any varargs if necessary
    nreq = 0
    if !isnothing(varargs)
        nreq = length(args)
        for a in varargs
            pusharg!(a)
        end
    end

    return func, rettype, types, args, nreq
end

function lower(convention, func, rettype, types, args, nreq)
    lowering = []
    realargs = []
    gcroots = []
    for (i, (arg, type)) in enumerate(zip(args, types))
        sym = Symbol(string("arg", i, "root"))
        sym2 = Symbol(string("arg", i, ))
        earg, etype = esc(arg), esc(type)
        push!(lowering, :($sym = Base.cconvert($etype, $earg)))
        push!(lowering, :($sym2 = Base.unsafe_convert($etype, $sym)))
        push!(realargs, sym2)
        push!(gcroots, sym)
    end
    etypes = Expr(:call, Expr(:core, :svec), types...)
    exp = Expr(:foreigncall,
               esc(func),
               esc(rettype),
               esc(etypes),
               nreq,
               QuoteNode(convention),
               realargs..., gcroots...)
    push!(lowering, exp)

    return Expr(:block, lowering...)
end

"""
    @ccall some_c_function(arg::Type [...])::ReturnType
    @ccall calling_convetion some_c_function(arg::Type)::ReturnType
convert a julia-style function definition to a ccall:
    @ccall link(source::Cstring, dest::Cstring)::Cint
same as:
    ccall(:link, Cint, (Cstring, Cstring), source, dest)
All arguments must have type annotations and the return type must also
be annotated.
varargs are supported with the following convention:
    @ccall printf("%s = %d"::Cstring ; "foo"::Cstring, foo::Cint)::Cint
Mind the semicolon.
Using functions from other libraries is supported by prefixing
the function name with the name of the C library, like this:
    const glib = "libglib-2.0"
    @ccall glib.g_uri_escape_string(uri::Cstring, ":/"::Cstring, true::Cint)::Cstring
The string literal could also be used directly before the symbol of
the function name, if desired `"libglib-2.0".g_uri_escape_string(...`
"""
macro ccall(convention, expr)
    return lower(convention, parsecall(expr)...)
end

macro ccall(expr)
    return lower(:ccall, parsecall(expr)...)
end

end # CcallMacro

using .CcallMacro: @ccall