text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Nemo

import Base: Array, abs, asin, asinh, atan, atanh, base, bin, checkbounds,
             conj, convert, cmp, contains, cos, cosh, dec, deepcopy,
             deepcopy_internal, deserialize, det, div, divrem, eye,
             gcd, gcdx, getindex, hash, hcat, hex, intersect, inv, invmod,
             isequal, isfinite, isless, isqrt, isreal, iszero, lcm,
             ldexp, length, log, lufact, lufact!, mod, ndigits, nextpow2, norm,
             nullspace, numerator, oct, one, parent, parse, precision, prevpow2,
             rand, rank, Rational, rem, reverse, serialize, setindex!, show,
             similar, sign, sin, sinh, size, sqrt, string, tan, tanh, trace,
             trailing_zeros, transpose, transpose!, truncate, typed_hvcat,
             typed_hcat, var, vcat, zero, zeros, +, -, *, ==, ^, &, |, $, <<,
             >>, ~, <=, >=, <, >, //, /, !=

if VERSION >= v"0.7.0-DEV.1144"
import Base: isone
end

import Base: floor, ceil, hypot, sqrt, log, log1p, expm1, sin, cos, sinpi,
             cospi, tan, cot, sinh, cosh, tanh, coth, atan, asin, acos, atanh,
             asinh, acosh, gamma, lgamma, sinpi, cospi, atan2

export elem_type, parent_type

export SetElem, GroupElem, RingElem, FieldElem, AccessorNotSetError

export PolyElem, SeriesElem, AbsSeriesElem, RelSeriesElem, ResElem, FracElem,
       MatElem, FinFieldElem, MPolyElem

export PolyRing, SeriesRing, AbsSeriesRing, ResRing, FracField, MatSpace,
       FinField, MPolyRing

export JuliaZZ, JuliaQQ, zz, qq, JuliaRealField, RDF

export PermutationGroup, ZZ, QQ, PadicField, FiniteField, RealField, ComplexField,
       CyclotomicField, MaximalRealSubfield, NumberField

export create_accessors, get_handle, package_handle, zeros,
       Array, sig_exists

export flint_cleanup, flint_set_num_threads

export error_dim_negative, ErrorConstrDimMismatch

export is_windows64

if VERSION >= v"0.6.0-dev.2024" # julia started exporting iszero (again?)
   import Base: iszero
end

if VERSION < v"0.6-"
   import Base: isprime, factor, parity, sub, call
end

if VERSION >= v"0.7.0-DEV.264" # julia started exporting sincos
   import Base: sincos
end

if VERSION >= v"0.7.0-DEV.1144"
    import Base: isone
end

include("AbstractTypes.jl")

###############################################################################
#
#   Set up environment / load libraries
#
###############################################################################

is_windows64() = (Sys.iswindows() ? true : false) && (Int == Int64)

const pkgdir = realpath(joinpath(dirname(@__FILE__), ".."))
const libdir = joinpath(pkgdir, "local", "lib")
if Sys.iswindows()
   const libgmp = joinpath(pkgdir, "local", "lib", "libgmp-16")
else
   const libgmp = joinpath(pkgdir, "local", "lib", "libgmp")
end
const libmpfr = joinpath(pkgdir, "local", "lib", "libmpfr")
const libflint = joinpath(pkgdir, "local", "lib", "libflint")
const libarb = joinpath(pkgdir, "local", "lib", "libarb")

function flint_abort()
  error("Problem in the Flint-Subsystem")
end

function __init__()

   if "HOSTNAME" in keys(ENV) && ENV["HOSTNAME"] == "juliabox"
       push!(Libdl.DL_LOAD_PATH, "/usr/local/lib")
   elseif Sys.islinux()
       push!(Libdl.DL_LOAD_PATH, libdir)
       Libdl.dlopen(libgmp)
       Libdl.dlopen(libmpfr)
       Libdl.dlopen(libflint)
       Libdl.dlopen(libarb)
   else
      push!(Libdl.DL_LOAD_PATH, libdir)
   end

   if !Sys.iswindows()
      ccall((:__gmp_set_memory_functions, libgmp), Void,
         (Ref{Void},Ref{Void},Ref{Void}),
         cglobal(:jl_gc_counted_malloc),
         cglobal(:jl_gc_counted_realloc_with_old_size),
         cglobal(:jl_gc_counted_free))

      ccall((:__flint_set_memory_functions, libflint), Void,
         (Ref{Void},Ref{Void},Ref{Void},Ref{Void}),
         cglobal(:jl_malloc),
         cglobal(:jl_calloc),
         cglobal(:jl_realloc),
         cglobal(:jl_free))
   end

   ccall((:flint_set_abort, libflint), Void,
      (Ref{Void},), cfunction(flint_abort, Void, ()))

   println("")
   println("Welcome to Nemo version 0.7.0")
   println("")
   println("Nemo comes with absolutely no warranty whatsoever")
   println("")
end

function flint_set_num_threads(a::Int)
   ccall((:flint_set_num_threads, libflint), Void, (Int,), a)
end

function flint_cleanup()
   ccall((:flint_cleanup, libflint), Void, ())
end

###############################################################################
#
#  Version information
#
################################################################################

function versioninfo()
  print("Nemo version 0.7.0 \n")
  nemorepo = dirname(dirname(@__FILE__))

  print("Nemo: ")
  prepo = Base.LibGit2.GitRepo(nemorepo)
  Base.LibGit2.with(LibGit2.head(prepo)) do phead
    print("commit: ")
    print(string(LibGit2.Oid(phead))[1:8])
    print(" date: ")
    commit = Base.LibGit2.get(Base.LibGit2.GitCommit, prepo, LibGit2.Oid(phead))
    print(Base.Dates.unix2datetime(Base.LibGit2.author(commit).time))
    print(")\n")
  end

  finalize(prepo)

  for deps in ["flint2", "arb", "antic"]
    if ispath(joinpath(nemorepo, "deps", deps))
      print("$deps: ")
      repo = joinpath(nemorepo, "deps", deps)

      prepo = Base.LibGit2.GitRepo(repo)

      Base.LibGit2.with(LibGit2.head(prepo)) do phead
        print("commit: ")
        print(string(LibGit2.Oid(phead))[1:8])
        print(" date: ")
        commit = Base.LibGit2.get(Base.LibGit2.GitCommit, prepo, LibGit2.Oid(phead))
        print(Base.Dates.unix2datetime(Base.LibGit2.author(commit).time))
        print(")\n")
      end
      finalize(prepo)
    end
  end

  return nothing
end

###############################################################################
#
#   Julia types
#
###############################################################################

include("julia/JuliaTypes.jl")

###############################################################################
#
#   Generic submodule
#
###############################################################################

include("Generic.jl")

import .Generic: add!, addeq!, addmul!, base_ring, canonical_unit, character,
                 characteristic, charpoly, charpoly_danilevsky!,
                 charpoly_danilevsky_ff!, charpoly_hessenberg!, chebyshev_t,
                 chebyshev_u, _check_dim, check_parent, coeff, cols, compose,
                 content, cycles, data, degree, denominator, derivative, det, det_clow,
                 det_df, det_fflu, det_popov, dim, discriminant, divexact,
                 divides, divrem, elem_type, elements, evaluate,
                 extended_weak_popov, extended_weak_popov_with_trafo, fflu!,
                 fflu, find_pivot_popov, fit!, gcd, gen, gens, gcdinv, gcdx,
                 gram, has_left_neighbor, has_bottom_neighbor, hash,
                 hessenberg!, hessenberg, hnf, hnf_cohen, hnf_cohen_with_trafo,
                 hnf_kb, hnf_kb_with_trafo, hnf_minors, hnf_minors_with_trafo,
                 hnf_with_trafo, hnf_via_popov, hnf_via_popov_with_trafo,
                 hooklength, identity_matrix, inskewdiag, integral,
                 interpolate, inv, inv!, invmod, isconstant, isdegree,
                 isdomain_type, isexact_type, isgen, ishessenberg, ismonomial,
                 isnegative, isone, isreverse, isrimhook, isrref, isterm,
                 isunit, lcm, lead, length, leglength, main_variable,
                 main_variable_extract, main_variable_insert, matrix,
                 matrix_repr, max_degrees, max_precision, minpoly, mod,
                 modulus, monomial_iszero, monomial_set!, monomial_to_newton!,
                 mul!, mul_classical, mul_karatsuba, mul_ks, mullow, mulmod,
                 needs_parentheses, newton_to_monomial!, normalise, nvars, numerator,
                 O, one, order, ordering, parent_type, parity, partitionseq,
                 polcoeff, pol_length, powmod, pow_multinomial, popov, powers,
                 precision, primpart, pseudodivrem, pseudorem, rand_ordering,
                 rank_profile_popov, remove, renormalize!, resultant,
                 resultant_ducos, resultant_euclidean, resultant_subresultant,
                 resultant_sylvester, resx, reverse, rows, rref, rref!,
                 setcoeff!, set_length!, setpermstyle, set_prec!, set_val!,
                 shift_left, shift_right, show_minus_one, similarity!, snf,
                 snf_kb, snf_kb_with_trafo, snf_with_trafo, solve,
                 solve_rational, solve_triu, subst, swap_rows, swap_rows!,
                 trail, truncate, typed_hcat, typed_hvcat, valuation, var,
                 vars, weak_popov, weak_popov_with_trafo, zero, zero!,
                 zero_matrix

export add!, addeq!, addmul!, base_ring, canonical_unit, character,
                 characteristic, charpoly, charpoly_danilevsky!,
                 charpoly_danilevsky_ff!, charpoly_hessenberg!, chebyshev_t,
                 chebyshev_u, _check_dim, check_parent, coeff, cols, compose,
                 content, cycles, data, degree, denominator, derivative, det, det_clow,
                 det_df, det_fflu, det_popov, dim, discriminant, divexact,
                 divides, divrem, elem_type, elements, evaluate,
                 extended_weak_popov, extended_weak_popov_with_trafo, fflu!,
                 fflu, find_pivot_popov, fit!, gcd, gen, gens, gcdinv, gcdx,
                 gram, has_left_neighbor, has_bottom_neighbor, hash,
                 hessenberg!, hessenberg, hnf, hnf_cohen, hnf_cohen_with_trafo,
                 hnf_kb, hnf_kb_with_trafo, hnf_minors, hnf_minors_with_trafo,
                 hnf_with_trafo, hnf_via_popov, hnf_via_popov_with_trafo,
                 hooklength, identity_matrix, inskewdiag, integral,
                 interpolate, inv, inv!, invmod, isconstant, isdegree,
                 isdomain_type, isexact, isexact_type, isgen, ishessenberg,
                 ismonomial, isnegative, isone, isreverse, isrimhook, isrref,
                 isterm, isunit, iszero, lcm, lead, leglength, length,
                 main_variable, main_variable_extract, main_variable_insert,
                 matrix, matrix_repr, max_degrees, max_precision, minpoly, mod,
                 modulus, monomial_iszero, monomial_set!, monomial_to_newton!,
                 mul!, mul_classical, mul_karatsuba, mul_ks, mullow, mulmod,
                 needs_parentheses, newton_to_monomial!, normalise, nvars, numerator,
                 O, one, order, ordering, parent_type, parity, partitionseq,
                 polcoeff, pol_length, powmod, pow_multinomial, popov, powers,
                 ppio, precision, primpart, pseudodivrem, pseudorem,
                 rand_ordering, rank_profile_popov, remove, renormalize!,
                 resultant, resultant_ducos, resultant_euclidean,
                 resultant_subresultant, resultant_sylvester, resx, reverse,
                 rows, rref, rref!, setcoeff!, set_length!, setpermstyle,
                 set_prec!, set_val!, shift_left, shift_right, show_minus_one,
                 similarity!, snf, snf_kb, snf_kb_with_trafo, snf_with_trafo,
                 solve, solve_rational, solve_triu, subst, swap_rows,
                 swap_rows!, trail, truncate, typed_hcat, typed_hvcat,
                 valuation, var, vars, weak_popov, weak_popov_with_trafo, zero,
                 zero!, zero_matrix

function exp(a::T) where T
   return Base.exp(a)
end

function PermGroup(n::Int, cached=true)
   Generic.PermGroup(n, cached)
end

function AllPerms(n::Int)
   Generic.AllPerms(n)
end

function perm(n::Int)
   Generic.perm(n)
end

function perm(a::Array{Int, 1})
   Generic.perm(a)
end

function Partition(part::Vector{Int}, check::Bool=true)
   Generic.Partition(part, check)
end

function Partitions(n::Int)
   Generic.Partitions(n)
end

function SkewDiagram(lambda::Generic.Partition, mu::Generic.Partition)
   Generic.SkewDiagram(lambda, mu)
end

function SkewDiagram(lambda::Vector{Int}, mu::Vector{Int})
   Generic.SkewDiagram(lambda, mu)
end

function YoungTableau(part::Generic.Partition, tab::Array{Int, 2})
   Generic.YoungTableau(part, tab)
end

function YoungTableau(part::Generic.Partition, fill::Vector{Int}=collect(1:part.n))
   Generic.YoungTableau(part, fill)
end

function YoungTableau(p::Vector{Int})
   Generic.YoungTableau(p)
end

function PowerSeriesRing(R::Ring, prec::Int, s::AbstractString; cached=true, model=:capped_relative)
   Generic.PowerSeriesRing(R, prec, s; cached=cached, model=model)
end

function LaurentSeriesRing(R::Ring, prec::Int, s::AbstractString; cached=true)
   Generic.LaurentSeriesRing(R, prec, s; cached=cached)
end

function LaurentSeriesRing(R::Field, prec::Int, s::AbstractString; cached=true)
   Generic.LaurentSeriesField(R, prec, s; cached=cached)
end

function LaurentSeriesField(R::Field, prec::Int, s::AbstractString; cached=true)
   Generic.LaurentSeriesField(R, prec, s; cached=cached)
end

function PolynomialRing(R::Ring, s::AbstractString; cached::Bool = true)
   Generic.PolynomialRing(R, s; cached=cached)
end

function PolynomialRing(R::Ring, s::Array{String, 1}; cached::Bool = true, ordering::Symbol = :lex)
   Generic.PolynomialRing(R, s; cached=cached, ordering=ordering)
end

function SparsePolynomialRing(R::Ring, s::String; cached::Bool = true)
   Generic.SparsePolynomialRing(R, s; cached=cached)
end

function MatrixSpace(R::Ring, r::Int, c::Int, cached::Bool = true)
   Generic.MatrixSpace(R, r, c, cached)
end

function FractionField(R::Ring; cached=true)
   Generic.FractionField(R; cached=cached)
end

function ResidueRing(R::Ring, a::Union{RingElement, Integer}; cached::Bool = true)
   Generic.ResidueRing(R, a; cached=cached)
end

export PowerSeriesRing, PolynomialRing, SparsePolynomialRing, MatrixSpace,
       FractionField, ResidueRing, Partition, PermGroup, YoungTableau,
       Partitions, SkewDiagram, AllPerms, perm, LaurentSeriesRing,
       LaurentSeriesField

export Generic

###############################################################################
#
#   Polynomial Ring S, x = R["x"] syntax
#
###############################################################################

getindex(R::Ring, s::String) = PolynomialRing(R, s)

getindex(R::Tuple{Ring, T}, s::String) where {T} = PolynomialRing(R[1], s)

###############################################################################
#
#   Matrix M = R[...] syntax
#
################################################################################

function typed_hvcat(R::Ring, dims, d...)
   T = elem_type(R)
   r = length(dims)
   c = dims[1]
   A = Array{T}(r, c)
   for i = 1:r
      dims[i] != c && throw(ArgumentError("row $i has mismatched number of columns (expected $c, got $(dims[i]))"))
      for j = 1:c
         A[i, j] = R(d[(i - 1)*c + j])
      end
   end
   S = matrix(R, A)
   return S
end

function typed_hcat(R::Ring, d...)
   T = elem_type(R)
   r = length(d)
   A = Array{T}(1, r)
   for i = 1:r
      A[1, i] = R(d[i])
   end
   S = matrix(R, A)
   return S
end

###############################################################################
#
#   Load error objects
#
###############################################################################

include("error.jl")

###############################################################################
#
#   Load Nemo Rings/Fields/etc
#
###############################################################################

include("flint/FlintTypes.jl")

include("antic/AnticTypes.jl")

include("arb/ArbTypes.jl")

#include("ambiguities.jl") # remove ambiguity warnings

include("Groups.jl")

include("flint/adhoc.jl")

###########################################################
#
#   Package handle creation
#
###########################################################

const package_handle = [1]

function get_handle()
   package_handle[1] += 1
   return package_handle[1] - 1
end

###############################################################################
#
#   Auxilliary data accessors
#
###############################################################################

mutable struct AccessorNotSetError <: Exception
end

function create_accessors(T, S, handle)
   get = function(a)
      if handle > length(a.auxilliary_data) ||
         !isassigned(a.auxilliary_data, handle)
        throw(AccessorNotSetError())
      end
      return a.auxilliary_data[handle]
   end
   set = function(a, b)
      if handle > length(a.auxilliary_data)
         resize!(a.auxilliary_data, handle)
      end
      a.auxilliary_data[handle] = b
   end
   return get, set
end

###############################################################################
#
#   Promote rule helpers
#
###############################################################################

if VERSION >= v"0.5.0-dev+3171"

function sig_exists(T::Type{Tuple{U, V, W}}, sig_table::Array{X, 1}) where {U, V, W, X}
   for s in sig_table
      if s === T
         return true
      end
   end
   return false
end

else

function sig_exists(T::Type{Tuple{U, V, W}}, sig_table::Array{X, 1}) where {U, V, W, X}
   return false
end

end # if VERSION

###############################################################################
#
#   Array creation functions
#
###############################################################################

Array(R::Ring, r::Int...) = Array{elem_type(R)}(r)

function zeros(R::Ring, r::Int...)
   T = elem_type(R)
   A = Array{T}(r)
   for i in eachindex(A)
      A[i] = R()
   end
   return A
end

###############################################################################
#
#   Set domain for PermutationGroup to Flint
#
###############################################################################

PermutationGroup = PermGroup

###############################################################################
#
#   Set domain for ZZ, QQ, PadicField, FiniteField to Flint
#
###############################################################################

ZZ = FlintZZ
QQ = FlintQQ
PadicField = FlintPadicField
FiniteField = FlintFiniteField

###############################################################################
#
#   Set domain for RR, CC to Arb
#
###############################################################################

RealField = ArbField
ComplexField = AcbField

###############################################################################
#
#   Set domain for NumberField to Antic
#
###############################################################################

NumberField = AnticNumberField
CyclotomicField = AnticCyclotomicField
MaximalRealSubfield = AnticMaximalRealSubfield

###############################################################################
#
#   Test code
#
###############################################################################

include("../benchmarks/runbenchmarks.jl")

function test_module(x, y)
   julia_exe = Base.julia_cmd()
   test_file = joinpath(pkgdir, "test/$x/")
   test_file = test_file * "$y-test.jl";
   test_function_name = "test_"

   if x in ["flint", "arb", "antic"]
     test_function_name *= y
   else x == "generic"
     if y == "RelSeries"
       test_function_name *= "gen_rel_series"
     elseif y == "AbsSeries"
       test_function_name *= "gen_abs_series"
     elseif y == "Matrix"
       test_function_name *= "gen_mat"
     elseif y == "Fraction"
       test_function_name *= "gen_frac"
     elseif y == "Residue"
       test_function_name *= "gen_res"
     else
       test_function_name *= "gen_$(lowercase(y))"
     end
   end

   cmd = "using Base.Test; using Nemo; include(\"$test_file\"); $test_function_name();"
   info("spawning ", `$julia_exe -e \"$cmd\"`)
   run(`$julia_exe -e $cmd`)
end

end # module


text:
abstract type Quadrature end
abstract type Euler <: Quadrature end
struct euler <: Euler end
abstract type RK4 <: Quadrature end
struct rk4 <: RK4 end

function integrate_dynamics(model, x0, us, ts, quadrature = euler())
    xs = Array{typeof(x0[1]),1}(undef, length(x0) * length(ts))
    integrate_dynamics!(xs, model, x0, us, ts, quadrature)
    return xs
end

"""
    integrate_dyanmics!(xs, model, x0, us, ts, quadrature)

Integrates the trajectory provided control inputs. Modifies `xs` in place.

# Inputs

* `xs::Vector{Float64}` - preallocated memory consisting of concatenated state vectors to be populated with simulation states for each element of `ts`
* `model<:AbstractModel` - model to be simulated
* `x0::Vector{Float64}` - initial state
* `us::Vector{Float64}` - concatenated control inputs prescribed for the simulation for each element of `ts`
* `ts::Vector{Float64}` - vector of times at which dynamics are evaluated
* `quadrature<:Quadrature` - what kind of quadrature is to be used for integration; `euler()` and `rk4()` currently supported

"""
function integrate_dynamics!(xs, model, x0, us, ts, quadrature::Euler)
    nx = model.nx # number of states
    @assert nx == length(x0) "length of x0 and model.nx inconsistent"
    nu = model.nu # number of controls
    @assert nu == Int(length(us) / length(ts)) "length of us and length of ts and model.nu inconsistent"
    @assert length(xs) == length(x0) * length(ts) "length of states x not consistent with timesteps and x0"

    # initial knot is the initial guess
    xs[1:nx] .= x0

    # loop through time
    for i=1:length(ts)-1
        Δt = ts[i+1] - ts[i]
        x = xs[(i-1) * nx + 1 : i * nx]
        u = us[(i-1) * nu + 1 : i * nu]
        xdot = dynamics(model, x, u)
        xs[i * nx + 1 : (i+1) * nx] .= x + xdot * Δt
    end

    return nothing
end

function integrate_dyanmics!(xs, model, x0, us, ts, quadrature::RK4)
    nx = model.nx # number of states
    @assert nx == length(x0) "length of x0 and model.nx inconsistent"
    nu = model.nu # number of controls
    @assert nu == Int(length(us) / length(ts)) "length of us and length of ts and model.nu inconsistent"
    @assert length(xs) == length(x0) * length(ts) "length of states x not consistent with timesteps and x0"

    # initial knot is the initial guess
    xs[1:nx] .= x0

    # loop through time
    for i=1:length(ts)-1
        Δt = ts[i+1] - ts[i]
        x = xs[(i-1) * nx + 1 : i * nx]
        u = us[(i-1) * nu + 1 : i * nu]
        u_next = us[i * nu + 1 : (i+1) * nu]
        F1 = dynamics(model, x, u)
        F2 = dynamics(model, x + .5 * Δt .* F1, u)
        F3 = dynamics(model, x + .5 * Δt .* F2, u_next)
        F4 = dynamics(model, x +  Δt .* F3, u_next)
        xs[i * nx + 1 : (i+1) * nx] .= x + Δt/6 .* (F1 + 2 .* F2 + 2 .* F3 + F4)
    end

    return nothing
end

ext: jl

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
"""
    coefficients(setup; D, T₂, ρ, κ, γ)

Prepare PDE compartments.
"""
function coefficients end

coefficients(::PlateSetup{T}; D, T₂, ρ, κ, γ) where {T} = (;
    D = SMatrix{3,3,T,9}.(D),
    T₂ = T.(T₂),
    κ = [T.(κ.interfaces); T.(κ.boundaries)],
    ρ = Complex{T}.(ρ),
    γ = T(γ),
)

function coefficients(setup::NeuronSetup{T}; D, T₂, ρ, κ, γ) where {T}
    (; ecs_shape, ecs_ratio) = setup

    include_ecs = ecs_shape != :no_ecs

    @assert !include_ecs || 0 < ecs_ratio

    # Determine number of compartments and boundaries
    ncompartment = 1 + include_ecs
    nboundary = 1 + include_ecs
    if include_ecs
        compartments = ["neuron", "ecs"]
        boundaries = ["neuron,ecs", "ecs"]
    else
        compartments = ["neuron"]
        boundaries = ["neuron"]
    end

    # Initialize output arrays
    coeffs = (;
        D = [zeros(SMatrix{3,3,T}) for _ = 1:ncompartment],
        T₂ = zeros(T, ncompartment),
        κ = zeros(T, nboundary),
        ρ = zeros(Complex{T}, ncompartment),
        γ = T(γ),
    )

    # Distribute material properties to compartments and boundaries
    coeffs.ρ[compartments.=="neuron"] .= ρ.neuron
    coeffs.ρ[compartments.=="ecs"] .= ρ.ecs
    coeffs.D[compartments.=="neuron"] .= [D.neuron]
    coeffs.D[compartments.=="ecs"] .= [D.ecs]
    coeffs.T₂[compartments.=="neuron"] .= T₂.neuron
    coeffs.T₂[compartments.=="ecs"] .= T₂.ecs
    coeffs.κ[boundaries.=="neuron,ecs"] .= κ.neuron_ecs
    coeffs.κ[boundaries.=="neuron"] .= κ.neuron
    coeffs.κ[boundaries.=="ecs"] .= κ.ecs

    coeffs
end

function coefficients(
    setup::Union{CylinderSetup{T},SphereSetup{T}};
    D,
    T₂,
    ρ,
    κ,
    γ,
) where {T}
    (; ncell, include_in, in_ratio, ecs_shape, ecs_ratio) = setup

    include_ecs = ecs_shape != :no_ecs

    # Check for correct radius ratios and that neurons do not have in-compartments
    @assert !include_in || 0 < in_ratio && in_ratio < 1
    @assert !include_ecs || 0 < ecs_ratio

    # Determine number of compartments and boundaries
    ncompartment = (1 + include_in) * ncell + include_ecs
    if isa(setup, SphereSetup)
        # For a sphere, there is one interface
        nboundary = (include_in + 1) * ncell + include_ecs
    elseif isa(setup, CylinderSetup)
        # An axon has a side interface, and a top-bottom boundary
        nboundary = (2 * include_in + 1 + include_ecs) * ncell + include_ecs
    end

    boundaries = String[]
    boundary_markers = fill(false, ncompartment, nboundary)
    if include_in
        # Add in-compartments and in-out-interfaces
        compartments = repeat(["in"], ncell)
        boundaries = repeat(["in,out"], ncell)
        boundary_markers[CartesianIndex.(1:2*ncell, repeat(1:ncell, 2))] .= true
        ncompartment_old = ncell
        nboundary_old = ncell
    else
        compartments = []
        boundaries = []
        ncompartment_old = 0
        nboundary_old = 0
    end

    # Add out-compartments
    append!(compartments, repeat(["out"], ncell))

    if include_ecs
        # Add ecs-compartment and out-ecs interfaces
        push!(compartments, "ecs")
        append!(boundaries, repeat(["out,ecs"], ncell))
        boundary_markers[CartesianIndex.(
            [ncompartment_old+1:ncompartment_old+ncell; fill(ncompartment, ncell)],
            repeat(nboundary_old+1:nboundary_old+ncell, 2),
        )] .= true
        nboundary_old = nboundary_old + ncell
    end

    if isa(setup, CylinderSetup)
        if include_in
            # Add in boundary
            append!(boundaries, repeat(["in"], ncell))
            ncompartment_old = ncell
            boundary_markers[CartesianIndex.(
                1:ncell,
                nboundary_old+1:nboundary_old+ncell,
            )] .= true
            nboundary_old = nboundary_old + ncell
        else
            ncompartment_old = 0
        end

        # Add out boundary
        append!(boundaries, repeat(["out"], ncell))
        boundary_markers[CartesianIndex.(
            ncompartment_old+1:ncompartment_old+ncell,
            nboundary_old+1:nboundary_old+ncell,
        )] .= true

        if include_ecs
            # Add ecs boundary
            push!(boundaries, "ecs")
            boundary_markers[end, end] = true
        end
    elseif include_ecs
        # Add ecs boundary
        push!(boundaries, "ecs")
        boundary_markers[end, end] = true
    else
        # Add out boundary
        push!(boundaries, "out")
        boundary_markers[end, end] = true
    end

    # Initialize output arrays
    coeffs = (;
        D = [zeros(SMatrix{3,3,T}) for _ = 1:ncompartment],
        T₂ = zeros(T, ncompartment),
        κ = zeros(T, nboundary),
        ρ = zeros(Complex{T}, ncompartment),
        γ = T(γ),
    )

    # Distribute material properties to compartments and boundaries
    coeffs.ρ[compartments.=="in"] .= ρ.in
    coeffs.ρ[compartments.=="out"] .= ρ.out
    coeffs.ρ[compartments.=="ecs"] .= ρ.ecs
    coeffs.D[compartments.=="in"] .= [D.in]
    coeffs.D[compartments.=="out"] .= [D.out]
    coeffs.D[compartments.=="ecs"] .= [D.ecs]
    coeffs.T₂[compartments.=="in"] .= T₂.in
    coeffs.T₂[compartments.=="out"] .= T₂.out
    coeffs.T₂[compartments.=="ecs"] .= T₂.ecs
    coeffs.κ[boundaries.=="in,out"] .= κ.in_out
    coeffs.κ[boundaries.=="out,ecs"] .= κ.out_ecs
    coeffs.κ[boundaries.=="in"] .= κ.in
    coeffs.κ[boundaries.=="out"] .= κ.out
    coeffs.κ[boundaries.=="ecs"] .= κ.ecs

    coeffs
end


text:
function getfield_elim_pass!(ir::IRCode)
    compact = IncrementalCompact(ir)
    insertions = Vector{Any}()
    for (idx, stmt) in compact
        is_call(stmt, :getfield) || continue
        isa(stmt.args[2], SSAValue) || continue
        field = stmt.args[3]
        isa(field, QuoteNode) && (field = field.value)
        isa(field, Union{Int, Symbol}) || continue
        orig_defidx = defidx = stmt.args[2].id
        def = compact[defidx]
        typeconstraint = types(compact)[defidx]
        phi_locs = Tuple{Int, Int}[]
        while true
            if isa(def, PiNode)
                typeconstraint = typeintersect(typeconstraint, def.typ)
                if isa(def.val, SSAValue)
                    defidx = def.val.id
                    def = compact[defidx]
                else
                    def = def.val
                end
                continue
            elseif isa(def, PhiNode)
                possible_predecessors = collect(Iterators.filter(1:length(def.edges)) do n
                    isassigned(def.values, n) || return false
                    value = def.values[n]
                    edge_typ = value_typ(compact, value)
                    return (NI.:⊑)(edge_typ, typeconstraint)
                end)
                # For now, only look at unique predecessors
                if length(possible_predecessors) == 1
                    n = possible_predecessors[1]
                    pred = def.edges[n]
                    val = def.values[n]
                    if isa(val, SSAValue)
                        push!(phi_locs, (pred, defidx))
                        defidx = val.id
                        def = compact[defidx]
                    else
                        def = val
                    end
                    continue
                end
            end
            break
        end
        if is_call(def, :tuple) && isa(field, Int) && 1 <= field < length(def.args)
            forwarded = def.args[1+field]
        elseif isexpr(def, :new)
            typ = def.typ
            if isa(typ, UnionAll)
                typ = Base.unwrap_unionall(typ)
            end
            isa(typ, DataType) || continue
            !typ.mutable || continue
            if isa(field, Symbol)
                field = Base.fieldindex(typ, field, false)
                field == 0 && continue
            end
            forwarded = def.args[1+field]
        else
            continue
        end
        if !isempty(phi_locs) && isa(forwarded, SSAValue)
            # TODO: We have have to use BB ids for phi_locs
            # to avoid index invalidation.
            push!(insertions, (idx, phi_locs))
        end
        compact[idx] = forwarded
    end
    ir = finish(compact)
    for (idx, phi_locs) in insertions
        # For non-dominating load-store forward, we may have to insert extra phi nodes
        # TODO: Can use the domtree to eliminate unnecessary phis, but ok for now
        forwarded = ir.stmts[idx]
        forwarded_typ = ir.types[forwarded.id]
        for (pred, pos) in reverse!(phi_locs)
            node = PhiNode()
            push!(node.edges, pred)
            push!(node.values, forwarded)
            forwarded = insert_node!(ir, pos, forwarded_typ, node)
        end
        ir.stmts[idx] = forwarded
    end
    ir
end

function get_val_if_type_cmp(def)
    is_call(def, :(===)) || return nothing
    (val, cmp) = def.args[2:3]
    if !isa(val, SSAValue)
        (cmp, val) = (val, cmp)
    end
    isa(val, SSAValue) || return nothing
    isa(cmp, SSAValue) && return nothing
    if isa(cmp, GlobalRef)
        cmp = getfield(cmp.mod, cmp.name)
    end
    isdefined(typeof(cmp), :instance) || return nothing
    return (val, cmp)
end

function type_lift_pass!(ir::IRCode)
    type_ctx_uses = Vector{Vector{Int}}[]
    has_non_type_ctx_uses = Set{Int}()
    lifted_undef = Dict{Int, SSAValue}()
    for (idx, stmt) in pairs(ir.stmts)
        if isexpr(stmt, :isdefined) || isexpr(stmt, :undefcheck)
            val = isexpr(stmt, :isdefined) ? stmt.args[1] : stmt.args[2]
            # undef can only show up by being introduced in a phi
            # node, so lift all phi nodes that have maybe undef values
            processed = Dict{Int, SSAValue}()
            worklist = Tuple{Int, SSAValue, Int}[(val.id, SSAValue(0), 0)]
            stmt_id = val.id
            while isa(ir.stmts[stmt_id], PiNode)
                stmt_id = ir.stmts[stmt_id].val.id
            end
            def = ir.stmts[stmt_id]
            if !isa(def, PhiNode)
                if isexpr(stmt, :isdefiend)
                    ir.stmts[idx] = true
                else
                    ir.stmts[idx] = nothing
                end
                continue
            end
            if !haskey(lifted_undef, stmt_id)
                first = true
                while !isempty(worklist)
                    item, which, use = pop!(worklist)
                    edges = copy(def.edges)
                    values = Vector{Any}(uninitialized, length(edges))
                    new_phi = insert_node!(ir, item, Bool, PhiNode(edges, values))
                    processed[item] = new_phi
                    if first
                        lifted_undef[stmt_id] = new_phi
                        first = false
                    end
                    for i = 1:length(edges)
                        if !isassigned(def.values, i)
                            val = false
                        elseif !isa(def.values[i], SSAValue)
                            val = true
                        else
                            id = def.values[i].id
                            if !isa(ir.types[id], NI.MaybeUndef)
                                val = true
                            else
                                while isa(ir.stmts[id], PiNode)
                                    id = ir.stmts[id].val.id
                                end
                                if isa(ir.stmts[id], PhiNode)
                                    if haskey(processed, id)
                                        val = processed[id]
                                    else
                                        push!(worklist, (id, new_phi, i))
                                        continue
                                    end
                                else
                                    val = true
                                end
                            end
                        end
                        values[i] = val
                    end
                    if which !== SSAValue(0)
                        ir[which].values[use] = new_phi
                    end
                end
            end
            if isexpr(stmt, :isdefiend)
                ir.stmts[idx] = lifted_undef[stmt_id]
            else
                ir.stmts[idx] = Expr(:throw_undef_if_not, stmt.args[1], lifted_undef[stmt_id])
            end
        else
            x = get_val_if_type_cmp(stmt)
            x === nothing && continue
            val, cmp = x
            cmptyp = typeof(cmp)
            def = ir.stmts[val.id]
            if isa(def, PhiNode)
                # See if this is only true on one branch
                branches = collect(Iterators.filter(1:length(def.edges)) do n
                    isassigned(def.values, n) || return false
                    value = def.values[n]
                    (NI.:⊑)(cmptyp, value_typ(ir, value))
                end)
                length(branches) == 1 || continue
                value_typ(ir, def.values[branches[1]]) == cmptyp || continue
                # Ok, merge the compare into the phi node
                node = PhiNode()
                for edge in def.edges
                    push!(node.edges, edge)
                    push!(node.values, edge == def.edges[branches[1]])
                end
                ir.stmts[idx] = insert_node!(ir, val.id, Bool, node)
            end
        end
    end
    ir
end


text:
using AbstractTrees, Test

@testset "Builtins" begin include("builtins.jl") end
@testset "Custom tree types" begin include("trees.jl") end
@testset "Printing" begin include("printing.jl") end


text:
# Random Sequence Generator
# =========================
#
# Random sequence generator.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Abstract sequence generator type.
"""
abstract SequenceGenerator{T}

# Sequence generator of stationary distributions.
immutable StationaryGenerator{T} <: SequenceGenerator{T}
    elems::Vector{T}
    probs::Vector{Float64}

    function StationaryGenerator(elems, probs)
        @assert sum(probs) ≤ 1
        @assert length(elems) == length(probs) + 1
        return new(copy(elems), copy(probs))
    end
end

function StationaryGenerator{T,U<:AbstractFloat}(
        elems::AbstractVector{T}, probs::AbstractVector{U})
    return StationaryGenerator{T}(elems, probs)
end

alphatype(::Type{DNA}) = DNAAlphabet{4}
alphatype(::Type{RNA}) = RNAAlphabet{4}
alphatype(::Type{AminoAcid})     = AminoAcidAlphabet

function randseq{T}(generator::StationaryGenerator{T}, len::Integer)
    if len < 0
        throw(ArgumentError("length must be non-negative"))
    end
    seq = BioSequence{alphatype(T)}(len)
    probs = generator.probs
    for i in 1:len
        r = rand()
        j = 1
        cumprob = probs[j]
        while cumprob < r && j ≤ endof(probs)
            j += 1
            if j ≤ endof(probs)
                cumprob += probs[j]
            end
        end
        seq[i] = generator.elems[j]
    end
    return seq
end

const StationaryDNAGenerator = StationaryGenerator(
    collect(dna"ACGT"), ones(3) / 4)

const StationaryRNAGenerator = StationaryGenerator(
    collect(rna"ACGU"), ones(3) / 4)

const StationaryAAGenerator = StationaryGenerator(
    collect(aa"ARNDCQEGHILKMFPSTWYV"), ones(19) / 20)

"""
    randdnaseq(len::Integer)

Generate a random DNA sequence of length `len`.
"""
function randdnaseq(len::Integer)
    return randseq(StationaryDNAGenerator, len)
end

"""
    randrnaseq(len::Integer)

Generate a random RNA sequence of length `len`.
"""
function randrnaseq(len::Integer)
    return randseq(StationaryRNAGenerator, len)
end

"""
    randaaseq(len::Integer)

Generate a random amino acid sequence of length `len`.
"""
function randaaseq(len::Integer)
    return randseq(StationaryAAGenerator, len)
end


text:
using Test
using JSON
using MianalyzerBackend
const BasicStatisticsController = MianalyzerBackend.BasicStatisticsController

include("../TestSupport.jl")

@testset "BasicStatisticsController" begin
  @testset ".call" begin
    request = Dict([("records", [[100.0], [200.0], [300.0]])])

    actual = JSON.parse(IOBuffer(BasicStatisticsController.call(request)))

    @test actual["means"] == [200]
    @test actual["stds"] == [100]
    @test actual["maxs"] == [300]
    @test actual["mins"] == [100]
    @test actual["row_count"] == 3
  end
end


text:
# `backward` functions return gates storing derivatives

backward(g::XGate, ::AbstractMatrix) = g
backward(g::YGate, ::AbstractMatrix) = g
backward(g::ZGate, ::AbstractMatrix) = g

backward(g::HadamardGate, ::AbstractMatrix) = g

backward(g::SGate, ::AbstractMatrix) = g
backward(g::TGate, ::AbstractMatrix) = g

backward(g::SdagGate, ::AbstractMatrix) = g
backward(g::TdagGate, ::AbstractMatrix) = g

backward(g::MatrixGate, ::AbstractMatrix) = g

function backward(g::RxGate, Δ::AbstractMatrix)
    c = cos(g.θ[] / 2)
    s = sin(g.θ[] / 2)
    # using conjugated derivative matrix; factor 2 cancels 1/2 from θ/2
    RxGate(sum(real([-s im * c; im * c -s] .* Δ)))
end

function backward(g::RyGate, Δ::AbstractMatrix)
    c = cos(g.θ[] / 2)
    s = sin(g.θ[] / 2)
    # using conjugated derivative matrix; factor 2 cancels 1/2 from θ/2
    RyGate(sum(real([-s -c; c -s] .* Δ)))
end

function backward(g::RzGate, Δ::AbstractMatrix)
    # using conjugated derivative matrix; factor 2 cancels 1/2 from θ/2
    eθ = exp(im * g.θ[] / 2)
    RzGate(real(im * eθ * Δ[1, 1] - im * conj(eθ) * Δ[2, 2]))
end

function backward(g::RotationGate, Δ::AbstractMatrix)
    θ = norm(g.nθ)
    if θ == 0
        σ = Matrix{ComplexF64}[matrix(X), matrix(Y), matrix(Z)]
        RotationGate([real(sum(conj(-im * σ[i]) .* Δ)) for i in 1:3])
    else
        n = g.nθ / θ
        c = cos(θ / 2)
        s = sin(θ / 2)
        dRθ = -s * I - im * c * pauli_vector(n...)
        dn = (I - reshape(kron(n, n), 3, 3)) / θ
        # using conjugated derivative matrix
        RotationGate([2 * real(sum(conj(0.5 * n[i] * dRθ - im * s * pauli_vector(dn[:, i]...)) .* Δ)) for i in 1:3])
    end
end

function backward(g::PhaseShiftGate, Δ::AbstractMatrix)
    # using conjugated derivative matrix
    PhaseShiftGate(2 * real(-im * Base.exp(-im * g.ϕ[1]) * Δ[2, 2]))
end

backward(g::SwapGate, ::AbstractMatrix) = g

function backward(g::EntanglementXXGate, Δ::AbstractMatrix)
    c = cos(g.θ[] / 2)
    s = sin(g.θ[] / 2)
    # using conjugated derivative matrix; factor 2 cancels 1/2 from θ/2
    EntanglementXXGate(sum(real([-s 0 0 im * c; 0 -s im * c 0; 0 im * c -s 0; im * c 0 0 -s] .* Δ)))
end

function backward(g::EntanglementYYGate, Δ::AbstractMatrix)
    c = cos(g.θ[] / 2)
    s = sin(g.θ[] / 2)
    # using conjugated derivative matrix; factor 2 cancels 1/2 from θ/2
    EntanglementYYGate(sum(real([-s 0 0 -im * c; 0 -s im * c 0; 0 im * c -s 0; -im * c 0 0 -s] .* Δ)))
end

function backward(g::EntanglementZZGate, Δ::AbstractMatrix)
    eθ = exp(im * g.θ[] / 2)
    # using conjugated derivative matrix; factor 2 cancels 1/2 from θ/2
    EntanglementZZGate(sum(real(im * eθ * Δ[1, 1] - im * conj(eθ) * Δ[2, 2] - im * conj(eθ) * Δ[3, 3] + im * eθ * Δ[4, 4])))
end

function backward(g::ControlledGate{G}, Δ::AbstractMatrix) where {G}
    # Note: target qubits correspond to fastest varying indices
    ControlledGate(backward(g.U, Δ[end-2^target_wires(g)+1:end, end-2^target_wires(g)+1:end]), control_wires(g))
end

function backward(cg::CircuitGate{M,G}, ψ::AbstractVector, Δ::AbstractVector, N::Int) where {M,G}
    ρ = rdm(N, cg.iwire, Δ, ψ)
    CircuitGate{M,G}(cg.iwire, backward(cg.gate, ρ))
end

function backward(m::Moment, ψ::AbstractVector, Δ::AbstractVector, N::Int)
    gates = AbstractCircuitGate[]
    for cg in reverse(m.gates)
        Udag = Base.adjoint(cg)
        ψ = apply(ψ, Udag)
        # backward step of quantum state
        pushfirst!(gates, backward(cg, ψ, Δ, N))
        Δ = apply(Δ, Udag)
    end
    return Moment(gates), ψ, Δ
end

function backward(moments::Vector{Moment}, ψ::AbstractVector, Δ::AbstractVector, N::Int)
    dmoments = Moment[]
    for moment in reverse(moments)
        # backward step of quantum state
        (m, ψ, Δ) = backward(moment, ψ, Δ, N)
        pushfirst!(dmoments, m)
    end
    return dmoments, Δ
end

"""
    gradients(c::Circuit{N}, ψ::AbstractVector, Δ::AbstractVector{<:Real}) where {N}

Perform a backward pass to compute gradients of a (fictitious) cost function with
respect to the circuit parameters of `c` and input wavefunction `ψ`. `Δ` contains
the cost function gradients with respect to the circuit outputs (measurement operator averages).
The gradients with respect to the circuit parameters are returned in a duplicate circuit;
the overall return value is the tuple (dc::Circuit{N}, dψ::AbstractVector).
"""
function gradients(c::Circuit{N}, ψ::AbstractVector, Δ::AbstractVector{<:Real}) where {N}
    # length of circuit output vector must match gradient vector
    @assert length(Δ) == length(c.meas)
    # forward pass through unitary gates
    ψ = apply(ψ, c.moments)
    # gradient (conjugated Wirtinger derivatives) of cost function with respect to ψ
    ψbar = sum([Δ[i] * (sparse_matrix(c.meas[i]) * ψ) for i in 1:length(Δ)])
    # backward pass through unitary gates
    dcgc, ψbar = backward(c.moments, ψ, ψbar, N)
    # TODO: efficiently represent Kronecker product without explicitly storing matrix entries
    dmeas = MeasurementOperator.([Δ[i] * reshape(kron(conj(ψ), ψ), length(ψ), length(ψ)) for i in 1:length(Δ)], (Tuple(1:N),))
    return Circuit{N}(dcgc, dmeas), ψbar
end


text:
const GAEVF = JuMP.GenericAffExpr{Float64, VariableRef}
const GQEVF = JuMP.GenericQuadExpr{Float64, VariableRef}

function moi_tests(op_problem::OperationsProblem,
                   params::Bool,
                   vars::Int64,
                   interval::Int64,
                   lessthan::Int64,
                   greaterthan::Int64,
                   equalto::Int64,
                   binary::Bool)

    JuMPmodel = op_problem.canonical.JuMPmodel
    @test (:params in keys(JuMPmodel.ext)) == params
    @test JuMP.num_variables(JuMPmodel) == vars
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.Interval{Float64}) == interval
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.LessThan{Float64}) == lessthan
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.GreaterThan{Float64}) == greaterthan
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.EqualTo{Float64}) == equalto
    @test ((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(JuMPmodel)) == binary

    return

end

function psi_constraint_test(op_problem::OperationsProblem, constraint_names::Vector{Symbol})

    for con in constraint_names
        @test !isnothing(get(op_problem.canonical.constraints, con, nothing))
    end

    return

end

function psi_checkbinvar_test(op_problem::OperationsProblem, bin_variable_names::Vector{Symbol})

    for variable in bin_variable_names
        for v in op_problem.canonical.variables[variable]
            @test JuMP.is_binary(v)
        end
    end

    return

end

function psi_checkobjfun_test(op_problem::OperationsProblem, exp_type)

    @test JuMP.objective_function_type(op_problem.canonical.JuMPmodel) == exp_type

    return

end

function moi_lbvalue_test(op_problem::OperationsProblem, con_name::Symbol, value::Number)

    for con in op_problem.canonical.constraints[con_name]
        @test JuMP.constraint_object(con).set.lower == value
    end

    return

end

function psi_checksolve_test(op_problem::OperationsProblem, status)
    JuMP.optimize!(op_problem.canonical.JuMPmodel)
    @test termination_status(op_problem.canonical.JuMPmodel) in status
end

function psi_checksolve_test(op_problem::OperationsProblem, status, expected_result, tol = 0.0)
    res = solve_op_problem!(op_problem)
    @test termination_status(op_problem.canonical.JuMPmodel) in status
    @test isapprox(res.total_cost[:OBJECTIVE_FUNCTION], expected_result, atol = tol)
end


text:
function test2(alpha::String)
  println(alpha)
end


text:
using BudgetConstraints

using Test

include("mini_tb_sys.jl")

# @test size( ps ) == 4
function test_getnet( data :: Dict, gross :: Float64 ) :: Float64
     person = data[:person]
     person.income = gross
     person.hours = gross/person.wage
     rc = calculate( person, data[:params] )
     return rc[:netincome]
end

function make_one_bc( person :: Person, params :: MiniParams ) :: BudgetConstraint
    data = Dict(
        :person=>person,
        :params=>params )
    bc = BudgetConstraints.makebc( data, test_getnet )
    return bc

end

@testset "reverse lookup thing" begin
    x1 = Point2D(0,10) # 100%
    x2 = Point2D(10,10) # 100%
    x3 = Point2D(20,20)  # 0%
    x4 = Point2D(40,30)  # 50%
    x5 = Point2D(100,45) # 75%
    bc = [x1,x2,x3,x4,x5]
    yn = get_x_from_y(bc,0.0) # some abitrary number between 0 and 10
    @test 0.0 >= yn <= 10.0

    yn = get_x_from_y(bc,11.0) # 0 mr so gross=net
    @test yn ≈ 11.0

    yn = get_x_from_y(bc,21.0) # 50% mr so 20 + 2
    @test yn ≈ 22.0

    yn = get_x_from_y(bc,31.0) # 75% mr so 40 + 4
    @test yn ≈ 44.0

end

@testset "BudgetConstraints.jl" begin

    using Test

    # @testset begin
     p1 = Point2D( 1.011,2.011)
     p2 = Point2D( 1.011, 2.011 )
     p3 = Point2D( 1.011, 3.011 )
     p4 = Point2D( -1.011, -3.011 )
     p5 = Point2D( 10.011, 13.011 )

     bc = BudgetConstraint( [p1,p2,p3,p4,p5])

    sort!( bc )

    @test size( bc )[1] == 5

    ps = Set([p1,p2,p3,p4,p5] )
    push!( ps, p1 )
    push!( ps, p2 )
    push!( ps, p3 )
    push!( ps, p4 )
    push!( ps, p5 )
    println( "ps = $ps ")
    push!( ps, Point2D( 1.0011, 1.999999 ))

    bc2 = censor( ps )

    println( "bc2 = $bc2 ")
    @test size( bc2 )[1] == 5
    pers = Person()
    pars = deepcopy(MiniParams())
    res = calculate( pers, pars )
    println( "res=$res" )

    bc = make_one_bc( Person(), MiniParams() )
    println( "\nDefault Case\n")
    println(BudgetConstraints.pointstoarray( bc))
    annotations = annotate_bc( bc )
    println( annotations )

    bc = make_one_bc( Person(), ZERO_PARAMS )
    println( "\nZero Params\n")
    println( BudgetConstraints.pointstoarray( bc))
    annotations = annotate_bc( bc )
    println( annotations )

    person = deepcopy(Person())
    # rc = calculate( Person(), pars )
    pars999 = MiniParams()
    pars999.it_allow = 999.0

    println( "\n999 allowance case\n")
    bc = make_one_bc( person, pars999 )
    println( BudgetConstraints.pointstoarray( bc ))
    annotations = annotate_bc( bc )
    println( annotations )

end


text:
export Bernoulli

import Distributions: Bernoulli, succprob, failprob

vague(::Type{ <: Bernoulli }) = Bernoulli(0.5)

probvec(dist::Bernoulli) = (succprob(dist), failprob(dist))

prod_analytical_rule(::Type{ <: Bernoulli }, ::Type{ <: Bernoulli }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::Bernoulli, right::Bernoulli)
    left_p  = succprob(left)
    right_p = succprob(right)

    pprod = left_p * right_p
    norm  = pprod + (one(left_p) - left_p) * (one(right_p) - right_p)
    @assert norm > 0 "Product of $(left) and $(right) results in non-normalizable distribution"
    return Bernoulli(pprod / norm)
end

prod_analytical_rule(::Type{ <: Bernoulli }, ::Type{ <: Categorical }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::Bernoulli, right::Categorical)
    @assert length(probvec(right)) === 2 "Improper Bernoulli x Categorical product"
    return prod(ProdPreserveType(Bernoulli), left, Bernoulli(first(probvec(right))))
end


text:
# compile the protobuf files

using ProtoBuf
ENV["JULIA_PROTOBUF_MODULE_POSTFIX"] = 1
run(ProtoBuf.protoc(`-I=proto --julia_out=src/proto proto/fileformat.proto proto/osmformat.proto`))


text:
"""
```
freak_params = FREAK([pattern_scale = 22.0])
```

| Argument | Type | Description |
|----------|------|-------------|
| **pattern_scale** | Float64 | Scaling factor for the sampling window |
"""
mutable struct FREAK <: Params
    pattern_scale::Float64
    pattern_table::Vector{Vector{SamplePair}}
    smoothing_table::Vector{Vector{Float16}}
    orientation_weights::Vector{OrientationWeights}
end

function FREAK(; pattern_scale::Float64 = 22.0)
    pattern_table, smoothing_table = _freak_tables(pattern_scale)
    orientation_weights = OrientationWeights[]
    for o in freak_orientation_sampling_pattern
        offset_1 = pattern_table[1][o[1]]
        offset_2 = pattern_table[1][o[2]]
        dy = offset_1[1] - offset_2[1]
        dx = offset_1[2] - offset_2[2]
        norm = (dx ^ 2 + dy ^ 2)
        push!(orientation_weights, OrientationWeights((dy / norm, dx / norm)))
    end
    FREAK(pattern_scale, pattern_table, smoothing_table, orientation_weights)
end

function _freak_mean_intensity(int_img::AbstractArray{T, 2}, keypoint::Keypoint, offset::SamplePair, sigma::Float16) where T<:Gray
    y = keypoint[1] + offset[1]
    x = keypoint[2] + offset[2]
    if sigma < 0.5
        sigma = 1.0
    end
    xs = round(Int, x - sigma)
    ys = round(Int, y - sigma)
    xst = round(Int, x + sigma)
    yst = round(Int, y + sigma)
    intensity = boxdiff(int_img, ys:yst, xs:xst)
    intensity / ((xst - xs + 1) * (yst - ys + 1))
end

function _freak_orientation(int_img::AbstractArray{T, 2}, keypoint::Keypoint, pattern::Array{SamplePair},
                               orientation_weights::Array{OrientationWeights}, sigmas::Array{Float16}) where T<:Gray
    direction_sum_y = 0.0
    direction_sum_x = 0.0
    for (i, o) in enumerate(freak_orientation_sampling_pattern)
        offset_1 = pattern[o[1]]
        offset_2 = pattern[o[2]]
        intensity_diff = _freak_mean_intensity(int_img, keypoint, offset_1, sigmas[o[1]]) - _freak_mean_intensity(int_img, keypoint, offset_2, sigmas[o[2]])
        direction_sum_y += orientation_weights[i][1] * intensity_diff
        direction_sum_x += orientation_weights[i][2] * intensity_diff
    end
    angle = atan(direction_sum_y, direction_sum_x)
    scaled_angle = ceil(Int, (angle + pi) * freak_orientation_steps / (2 * pi))
    scaled_angle
end

function _freak_tables(pattern_scale::Float64)
    pattern_table = Vector{SamplePair}[]
    smoothing_table = Vector{Float16}[]
    for ori in 0:freak_orientation_steps - 1
        theta = ori * 2 * pi / freak_orientation_steps
        pattern = SamplePair[]
        sigmas = Float16[]
        for (i, n) in enumerate(freak_num_circular_pattern)
            for circle_number in 0:n - 1
                alt_offset = (pi / n) * ((i - 1) % 2)
                angle = (circle_number * 2 * pi / n) + alt_offset + theta

                push!(pattern, SamplePair((freak_radii[i] * sin(angle) * pattern_scale,
                                            freak_radii[i] * cos(angle) * pattern_scale)))
                push!(sigmas, freak_sigma[i] * pattern_scale)
            end
        end
        push!(pattern_table, pattern)
        push!(smoothing_table, sigmas)
    end
    pattern_table, smoothing_table
end

function create_descriptor(img::AbstractArray{T, 2}, keypoints::Keypoints, params::FREAK) where T<:Gray
    int_img = integral_image(img)
    descriptors = BitArray{1}[]
    ret_keypoints = Keypoint[]
    window_size = ceil(Int, (freak_radii[1] + freak_sigma[1]) * params.pattern_scale) + 1
    lim = CartesianIndex(window_size, window_size)
    for k in keypoints
        checkbounds(Bool, img, k - lim) && checkbounds(Bool, img, k + lim) || continue
        orientation = _freak_orientation(int_img, k, params.pattern_table[1], params.orientation_weights, params.smoothing_table[1])
        sampled_intensities = T[]
        for (i, p) in enumerate(params.pattern_table[orientation])
            push!(sampled_intensities,  _freak_mean_intensity(int_img, k, p, params.smoothing_table[orientation][i]))
        end
        descriptor = falses(512)
        for (i, f) in enumerate(freak_sampling_pattern)
            point_1 = sampled_intensities[f[1]]
            point_2 = sampled_intensities[f[2]]
            descriptor[i] = point_1 < point_2
        end
        push!(descriptors, descriptor)
        push!(ret_keypoints, k)
    end
    descriptors, ret_keypoints
end


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

"""
A list of data file sections in the order that they appear in a PTI v33 file
"""
const _pti_sections = [
    "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",
]

const _transaction_dtypes = [
    ("IC", Int64),
    ("SBASE", Float64),
    ("REV", Int64),
    ("XFRRAT", Float64),
    ("NXFRAT", Float64),
    ("BASFRQ", Float64),
]

const _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),
]

const _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),
]

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

const _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),
]

const _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),
]

const _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),
]

const _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),
]

const _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),
]

const _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),
]

const _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),
]

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

const _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),
]

const _transformer_2_3_dtypes = [("WINDV2", Float64), ("NOMV2", Float64)]

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

const _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),
]

const _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),
]

const _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),
]

const _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),
]

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

const _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),
]

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

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

const _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),
]

const _zone_dtypes = [("I", Int64), ("ZONAME", String)]

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

const _owner_dtypes = [("I", Int64), ("OWNAME", String)]

const _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),
]

const _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
const _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),
]

const _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),
]

"""
lookup array of data types for PTI file sections given by
`field_name`, as enumerated by PSS/E Program Operation Manual.
"""
const _pti_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,
)

const _default_case_identification = Dict(
    "IC" => 0,
    "SBASE" => 100.0,
    "REV" => 33,
    "XFRRAT" => 0,
    "NXFRAT" => 0,
    "BASFRQ" => 60,
)

const _default_bus = Dict(
    "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" => "            ",
)

const _default_load = Dict(
    "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" => nothing,
    "ZONE" => nothing,
    "OWNER" => nothing,
)

const _default_fixed_shunt = Dict("ID" => "1", "STATUS" => 1, "GL" => 0.0, "BL" => 0.0)

const _default_generator = Dict(
    "ID" => "1",
    "PG" => 0.0,
    "QG" => 0.0,
    "QT" => 9999.0,
    "QB" => -9999.0,
    "VS" => 1.0,
    "IREG" => 0,
    "MBASE" => nothing,
    "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" => nothing,
    "O2" => 0,
    "O3" => 0,
    "O4" => 0,
    "F1" => 1.0,
    "F2" => 1.0,
    "F3" => 1.0,
    "F4" => 1.0,
    "WMOD" => 0,
    "WPF" => 1.0,
)

const _default_branch = Dict(
    "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" => nothing,
    "O2" => 0,
    "O3" => 0,
    "O4" => 0,
    "F1" => 1.0,
    "F2" => 1.0,
    "F3" => 1.0,
    "F4" => 1.0,
)

const _default_transformer = Dict(
    "K" => 0,
    "CKT" => "1",
    "CW" => 1,
    "CZ" => 1,
    "CM" => 1,
    "MAG1" => 0.0,
    "MAG2" => 0.0,
    "NMETR" => 2,
    "NAME" => "            ",
    "STAT" => 1,
    "O1" => nothing,
    "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" => nothing,
    "R2-3" => 0.0,
    "SBASE2-3" => nothing,
    "R3-1" => 0.0,
    "SBASE3-1" => nothing,
    "VMSTAR" => 1.0,
    "ANSTAR" => 0.0,
    "WINDV1" => nothing,
    "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" => nothing,
    "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" => nothing,
    "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,
)

const _default_area_interchange =
    Dict("ISW" => 0, "PDES" => 0.0, "PTOL" => 10.0, "ARNAME" => "            ")

const _default_two_terminal_dc = Dict(
    "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,
)

const _default_vsc_dc = Dict(
    "MDC" => 1,
    "O1" => nothing,
    "O2" => 0,
    "O3" => 0,
    "O4" => 0,
    "F1" => 1.0,
    "F2" => 1.0,
    "F3" => 1.0,
    "F4" => 1.0,
    "CONVERTER BUSES" => Dict(
        "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,
    ),
)

const _default_impedance_correction = Dict(
    "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,
)

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

const _default_multi_section = Dict("ID" => "&1", "MET" => 1)

const _default_zone = Dict("ZONAME" => "            ")

const _default_interarea = Dict("TRID" => 1, "PTRAN" => 0.0)

const _default_owner = Dict("OWNAME" => "            ")

const _default_facts = Dict(
    "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" => "",
)

const _default_switched_shunt = Dict(
    "MODSW" => 1,
    "ADJM" => 0,
    "STAT" => 1,
    "VSWHI" => 1.0,
    "VSWLO" => 1.0,
    "SWREM" => 0,
    "RMPCT" => 100.0,
    "RMIDNT" => "",
    "BINIT" => 0.0,
    "N1" => 0,
    "N2" => 0,
    "N3" => 0,
    "N4" => 0,
    "N5" => 0,
    "N6" => 0,
    "N7" => 0,
    "N8" => 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,
)

const _default_gne_device = Dict(
    "NTERM" => 1,
    "NREAL" => 0,
    "NINTG" => 0,
    "NCHAR" => 0,
    "STATUS" => 1,
    "OWNER" => nothing,
    "NMETR" => nothing,
    "REAL" => 0,
    "INTG" => nothing,
    "CHAR" => "1",
)

const _default_induction_machine = Dict(
    "ID" => 1,
    "STAT" => 1,
    "SCODE" => 1,
    "DCODE" => 2,
    "AREA" => nothing,
    "ZONE" => nothing,
    "OWNER" => nothing,
    "TCODE" => 1,
    "BCODE" => 1,
    "MBASE" => nothing,
    "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,
)

const _pti_defaults = Dict(
    "BUS" => _default_bus,
    "LOAD" => _default_load,
    "FIXED SHUNT" => _default_fixed_shunt,
    "GENERATOR" => _default_generator,
    "BRANCH" => _default_branch,
    "TRANSFORMER" => _default_transformer,
    "AREA INTERCHANGE" => _default_area_interchange,
    "TWO-TERMINAL DC" => _default_two_terminal_dc,
    "VOLTAGE SOURCE CONVERTER" => _default_vsc_dc,
    "IMPEDANCE CORRECTION" => _default_impedance_correction,
    "MULTI-TERMINAL DC" => _default_multi_term_dc,
    "MULTI-SECTION LINE" => _default_multi_section,
    "ZONE" => _default_zone,
    "INTER-AREA TRANSFER" => _default_interarea,
    "OWNER" => _default_owner,
    "FACTS CONTROL DEVICE" => _default_facts,
    "SWITCHED SHUNT" => _default_switched_shunt,
    "CASE IDENTIFICATION" => _default_case_identification,
    "GNE DEVICE" => _default_gne_device,
    "INDUCTION MACHINE" => _default_induction_machine,
)

function _correct_nothing_values!(data::Dict)
    if !haskey(data, "BUS")
        return
    end

    sbase = data["CASE IDENTIFICATION"][1]["SBASE"]
    bus_lookup = Dict(bus["I"] => bus for bus in data["BUS"])

    if haskey(data, "LOAD")
        for load in data["LOAD"]
            load_bus = bus_lookup[load["I"]]
            if load["AREA"] == nothing
                load["AREA"] = load_bus["AREA"]
            end
            if load["ZONE"] == nothing
                load["ZONE"] = load_bus["ZONE"]
            end
            if load["OWNER"] == nothing
                load["OWNER"] = load_bus["OWNER"]
            end
        end
    end

    if haskey(data, "GENERATOR")
        for gen in data["GENERATOR"]
            gen_bus = bus_lookup[gen["I"]]
            if haskey(gen, "OWNER") && gen["OWNER"] == nothing
                gen["OWNER"] = gen_bus["OWNER"]
            end
            if gen["MBASE"] == nothing
                gen["MBASE"] = sbase
            end
        end
    end

    if haskey(data, "BRANCH")
        for branch in data["BRANCH"]
            branch_bus = bus_lookup[branch["I"]]
            if haskey(branch, "OWNER") && branch["OWNER"] == nothing
                branch["OWNER"] = branch_bus["OWNER"]
            end
        end
    end

    if haskey(data, "TRANSFORMER")
        for transformer in data["TRANSFORMER"]
            transformer_bus = bus_lookup[transformer["I"]]
            for base_id in ["SBASE1-2", "SBASE2-3", "SBASE3-1"]
                if haskey(transformer, base_id) && transformer[base_id] == nothing
                    transformer[base_id] = sbase
                end
            end
            for winding_id in ["WINDV1", "WINDV2", "WINDV3"]
                if haskey(transformer, winding_id) && transformer[winding_id] == nothing
                    if transformer["CW"] == 2
                        transformer[winding_id] = transformer_bus["BASKV"]
                    else
                        transformer[winding_id] = 1.0
                    end
                end
            end
        end
    end

    #=
    # TODO update this default value
    if haskey(data, "VOLTAGE SOURCE CONVERTER")
        for mdc in data["VOLTAGE SOURCE CONVERTER"]
            mdc["O1"] = Expr(:call, :_get_component_property, data["BUS"], "OWNER", "I", get(get(component, "CONVERTER BUSES", [Dict()])[1], "IBUS", 0))
        end
    end
    =#

    if haskey(data, "GNE DEVICE")
        for gne in data["GNE DEVICE"]
            gne_bus = bus_lookup[gne["I"]]
            if haskey(gne, "OWNER") && gne["OWNER"] == nothing
                gne["OWNER"] = gne_bus["OWNER"]
            end
            if haskey(gne, "NMETR") && gne["NMETR"] == nothing
                gne["NMETR"] = gne_bus["NTERM"]
            end
        end
    end

    if haskey(data, "INDUCTION MACHINE")
        for indm in data["INDUCTION MACHINE"]
            indm_bus = bus_lookup[indm["I"]]
            if indm["AREA"] == nothing
                indm["AREA"] = indm_bus["AREA"]
            end
            if indm["ZONE"] == nothing
                indm["ZONE"] = indm_bus["ZONE"]
            end
            if indm["OWNER"] == nothing
                indm["OWNER"] = indm_bus["OWNER"]
            end
            if indm["MBASE"] == nothing
                indm["MBASE"] = sbase
            end
        end
    end
end

"""
This is an experimental method for parsing elements and setting defaults at the same time.
It is not currently working but would reduce memory allocations if implemented correctly.
"""
function _parse_elements(
    elements::Array,
    dtypes::Array,
    defaults::Dict,
    section::AbstractString,
)
    data = Dict{String, Any}()

    if length(elements) > length(dtypes)
        @warn(
            "ignoring $(length(elements) - length(dtypes)) extra values in section $section, only $(length(dtypes)) items are defined"
        )
        elements = elements[1:length(dtypes)]
    end

    for (i, element) in enumerate(elements)
        field, dtype = dtypes[i]

        element = strip(element)

        if dtype == String
            if startswith(element, "'") && endswith(element, "'")
                data[field] = element[2:(end - 1)]
            else
                data[field] = element
            end
        else
            if length(element) <= 0
                # this will be set to a default in the cleanup phase
                data[field] = nothing
            else
                try
                    data[field] = parse(dtype, element)
                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
        end
    end

    if length(elements) < length(dtypes)
        for (field, dtype) in dtypes[length(elements):end]
            data[field] = defaults[field]
            #=
            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"))
                    @warn("The following fields in $section are missing: $missing_str")
                end
            end
            =#
        end
    end

    return data
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 (i, (field, dtype)) in enumerate(_pti_dtypes[section])
        if i > length(elements)
            @debug "Have run out of elements in $section at $field" _group =
                IS.LOG_GROUP_PARSING
            push!(missing_fields, field)
            continue
        else
            element = strip(elements[i])
        end

        try
            if dtype != String && element != ""
                data[field] = parse(dtype, element)
            else
                if dtype == String && startswith(element, "'") && endswith(element, "'")
                    data[field] = chop(element[nextind(element, 1):end])
                else
                    data[field] = element
                end
            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

const _comment_split = r"(?!\B[\'][^\']*)[\/](?![^\']*[\']\B)"
const _split_string = r",(?=(?:[^']*'[^']*')*[^']*$)"

"""
    _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)
    if count(i -> (i == "'"), line) % 2 == 1
        throw(
            DataFormatError(
                "There are an uneven number of single-quotes in \"{line}\", the line cannot be parsed.",
            ),
        )
    end

    line_comment = split(line, _comment_split, limit = 2)
    line = strip(line_comment[1])
    comment = length(line_comment) > 1 ? strip(line_comment[2]) : ""

    elements = split(line, _split_string)

    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 = deepcopy(_pti_sections)
    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)
        (elements, comment) = _get_line_elements(line)

        first_element = strip(elements[1])
        if line_number > 3 && length(elements) != 0 && first_element == "Q"
            break
        elseif line_number > 3 && length(elements) != 0 && first_element == "0"
            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
                @debug "At line $line_number, switched to $section" _group =
                    IS.LOG_GROUP_PARSING
            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], " ") _group = IS.LOG_GROUP_PARSING
            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
                    subsection_data = Dict{String, Any}()

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

                    if haskey(section_data, "CONVERTER BUSES")
                        push!(section_data["CONVERTER BUSES"], subsection_data)
                    else
                        section_data["CONVERTER BUSES"] = [subsection_data]
                        continue
                    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

                    if haskey(section_data, "$(subsection[2:end])")
                        section_data["$(subsection[2:end])"] =
                            push!(section_data["$(subsection[2:end])"], subsection_data)
                        if skip_sublines > 0 && subsection != "NDCLN"
                            continue
                        end
                    else
                        section_data["$(subsection[2:end])"] = [subsection_data]
                        if skip_sublines > 0 && subsection != "NDCLN"
                            continue
                        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" _group = IS.LOG_GROUP_PARSING
        end

        if haskey(pti_data, section)
            push!(pti_data[section], section_data)
        else
            pti_data[section] = [section_data]
        end
    end

    _populate_defaults!(pti_data)
    _correct_nothing_values!(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)
    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 _pti_sections
        if haskey(data, section)
            component_defaults = _pti_defaults[section]
            for component in data[section]
                for (field, field_value) in component
                    if isa(field_value, Array)
                        sub_component_defaults = component_defaults[field]
                        for sub_component in field_value
                            for (sub_field, sub_field_value) in sub_component
                                if sub_field_value == ""
                                    try
                                        sub_component[sub_field] =
                                            sub_component_defaults[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
                        end
                    elseif field_value == "" &&
                           !(field in ["Comment_Line_1", "Comment_Line_2"]) &&
                           !startswith(field, "DUM")
                        try
                            component[field] = component_defaults[field]
                        catch msg
                            if isa(msg, KeyError)
                                @warn("'$field' in '$section' has no default value",)
                            else
                                rethrow(msg)
                            end
                        end
                    end
                end
            end
        end
    end
end