text:
using GeometricBase.Config
using GeometricBase.Utils
using GeometricEquations
using GeometricIntegrators.Integrators
using GeometricIntegrators.Integrators: initial_guess!, jacobian!, update_params!
using GeometricIntegrators.Solutions
using GeometricIntegrators.Tableaus
using GeometricProblems.HarmonicOscillator
using SimpleSolvers
using Test

using GeometricProblems.HarmonicOscillator: Δt, nt, refx, refq, refp, k

ode  = harmonic_oscillator_ode()
pode = harmonic_oscillator_pode()
sode = harmonic_oscillator_sode()
iode = harmonic_oscillator_iode()
idae = harmonic_oscillator_idae()
hode = harmonic_oscillator_hode()
pdae = harmonic_oscillator_pdae()

@testset "$(rpad("Runge-Kutta integrators",80))" begin

    @test typeof(Integrator(ode, TableauExplicitMidpoint(), Δt)) <: IntegratorERK
    @test typeof(Integrator(ode, TableauCrouzeix(), Δt)) <: IntegratorDIRK
    @test typeof(Integrator(ode, TableauImplicitMidpoint(), Δt)) <: IntegratorFIRK
    @test typeof(Integrator(pode, PartitionedTableau(:eprk_midpoint, TableauExplicitMidpoint()), Δt)) <: IntegratorEPRK
    @test typeof(Integrator(pode, PartitionedTableau(:iprk_midpoint, TableauImplicitMidpoint()), Δt)) <: IntegratorIPRK
    @test typeof(Integrator(iode, PartitionedTableau(:iprk_midpoint, TableauImplicitMidpoint()), Δt)) <: IntegratorPRKimplicit

    @test_logs (:warn, r"Initializing IntegratorDIRK with explicit tableau heun.*") IntegratorDIRK(ode, TableauHeun2(), Δt)
    @test_logs (:warn, r"Initializing IntegratorFIRK with explicit tableau heun.*") IntegratorFIRK(ode, TableauHeun2(), Δt)
    @test_logs (:warn, r"Initializing IntegratorFIRK with diagonally implicit tableau crouzeix.*") IntegratorFIRK(ode, TableauCrouzeix(), Δt)

    int_erk  = Integrator(ode, TableauExplicitMidpoint(), Δt)
    int_dirk = Integrator(ode, TableauCrouzeix(), Δt)
    int_firk = Integrator(ode, TableauImplicitMidpoint(), Δt)
    int_eprk = Integrator(pode, PartitionedTableau(:eprk_midpoint, TableauExplicitMidpoint()), Δt)
    int_iprk = Integrator(pode, PartitionedTableau(:iprk_midpoint, TableauImplicitMidpoint()), Δt)

    @test ndims(ode) == ndims(int_erk)
    @test ndims(ode) == ndims(int_dirk)
    @test ndims(ode) == ndims(int_firk)
    @test ndims(pode) == ndims(int_eprk)
    @test ndims(pode) == ndims(int_iprk)

end

@testset "$(rpad("Explicit Runge-Kutta integrators",80))" begin

    sol = integrate(ode, TableauExplicitEuler(), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-2

    sol = integrate(ode, TableauExplicitMidpoint(), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-3

    sol = integrate(ode, TableauRK4(), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-7

end

function test_firk_jacobian(ode, tableau, Δt, nt; atol=eps())
    int1 = IntegratorFIRK(ode, tableau, Δt; exact_jacobian=false)
    int2 = IntegratorFIRK(ode, tableau, Δt; exact_jacobian=true)

    sol = Solution(ode, Δt, nt)

    asol1 = AtomicSolution(ode)
    asol2 = AtomicSolution(ode)

    get_initial_conditions!(sol, asol1, 1, 1)
    get_initial_conditions!(sol, asol2, 1, 1)

    update_params!(int1, asol1)
    update_params!(int2, asol2)

    initial_guess!(int1, asol1)
    initial_guess!(int2, asol2)

    computeJacobian(int1.solver.x, int1.solver.J, int1.solver.Jparams)
    computeJacobian(int2.solver.x, int2.solver.J, int2.solver.Jparams)

    @test int1.solver.J ≈ int2.solver.J atol=atol
end

@testset "$(rpad("Implicit Runge-Kutta integrators",80))" begin

    # check Jacobian of FIRK integrators

    int = Integrator(ode, TableauImplicitMidpoint(), Δt)
    jacobian!(int.solver.x, int.solver.J, IntegratorCache(int.params), int.params)
    @test int.solver.J == - [1.0  0.0; 0.0  1.0] + Δt / 2 * [0.0  1.0; -k  0.0]

    test_firk_jacobian(ode, TableauGauss(1), Δt, nt)
    test_firk_jacobian(ode, TableauGauss(2), Δt, nt)
    test_firk_jacobian(ode, TableauGauss(3), Δt, nt)
    test_firk_jacobian(ode, TableauGauss(4), Δt, nt)

    # check integrators

    sol = integrate(ode, TableauCrouzeix(), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-5

    sol = integrate(ode, TableauImplicitEuler(), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-2

    sol = integrate(ode, TableauImplicitMidpoint(), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-4

    sol = integrate(ode, TableauGauss(1), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-4

    sol = integrate(ode, TableauGauss(2), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-7

    sol = integrate(ode, TableauGauss(3), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    sol = integrate(ode, TableauGauss(4), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    sol = integrate(ode, TableauGauss(5), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    sol = integrate(ode, TableauGauss(6), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

    sol = integrate(ode, TableauGauss(7), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-15

    sol = integrate(ode, TableauSRK3(), Δt, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-7

    # Test integration with exact Jacobian

    int = IntegratorFIRK(ode, TableauGauss(1), Δt; exact_jacobian=true)
    sol = integrate(ode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 5E-4

    int = IntegratorFIRK(ode, TableauGauss(2), Δt; exact_jacobian=true)
    sol = integrate(ode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-7

    int = IntegratorFIRK(ode, TableauGauss(3), Δt; exact_jacobian=true)
    sol = integrate(ode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-11

    int = IntegratorFIRK(ode, TableauGauss(4), Δt; exact_jacobian=true)
    sol = integrate(ode, int, nt)
    @test relative_maximum_error(sol.q, refx) < 1E-15

end

@testset "$(rpad("Partitioned Runge-Kutta integrators",80))" begin

    # psol = integrate(pode, TableauSymplecticEulerA(), Δt, nt)
    # @test relative_maximum_error(psol.q, refq) < 5E-2
    # @test relative_maximum_error(psol.p, refp) < 1E-3

    # psol = integrate(pode, TableauSymplecticEulerB(), Δt, nt)
    # @test relative_maximum_error(psol.q, refq) < 5E-2
    # @test relative_maximum_error(psol.p, refp) < 1E-3

    psol = integrate(pode, PartitionedTableau(:prk4, TableauRK4()), Δt, nt)
    @test relative_maximum_error(psol.q, refq) < 5E-7
    @test relative_maximum_error(psol.p, refp) < 5E-7

    psol = integrate(pode, PartitionedTableau(:pglrk, TableauGauss(1)), Δt, nt)
    @test relative_maximum_error(psol.q, refq) < 5E-4
    @test relative_maximum_error(psol.p, refp) < 5E-4

    psol = integrate(pode, PartitionedTableau(:pglrk, TableauGauss(2)), Δt, nt)
    @test relative_maximum_error(psol.q, refq) < 1E-7
    @test relative_maximum_error(psol.p, refp) < 1E-7

    psol = integrate(pode, PartitionedTableau(:pglrk, TableauGauss(3)), Δt, nt)
    @test relative_maximum_error(psol.q, refq) < 1E-11
    @test relative_maximum_error(psol.p, refp) < 1E-11

end

@testset "$(rpad("Special Runge-Kutta integrators",80))" begin

    pgint = IntegratorPGLRK(ode, CoefficientsPGLRK(2), Δt)
    pgsol = integrate(ode, pgint, nt)
    @test relative_maximum_error(pgsol.q, refx) < 6E-6

    pgint = IntegratorPGLRK(ode, CoefficientsPGLRK(3), Δt)
    pgsol = integrate(ode, pgint, nt)
    @test relative_maximum_error(pgsol.q, refx) < 2E-12

    pgint = IntegratorPGLRK(ode, CoefficientsPGLRK(4), Δt)
    pgsol = integrate(ode, pgint, nt)
    @test relative_maximum_error(pgsol.q, refx) < 8E-16

end

@testset "$(rpad("Integrate PODE and HODE with ODE Runge-Kutta integrators",80))" begin
    for s in 1:4
        code = convert(ODE, pode)
        csol = integrate(code, TableauGauss(s), Δt, nt)
        psol = integrate(pode, PartitionedTableauGauss(s), Δt, nt)

        @test csol.q[1, end] == psol.q[1, end]
        @test csol.q[2, end] == psol.p[1, end]
    end

    for s in 1:4
        code = convert(ODE, hode)
        csol = integrate(code, TableauGauss(s), Δt, nt)
        hsol = integrate(hode, PartitionedTableauGauss(s), Δt, nt)

        @test csol.q[1, end] == hsol.q[1, end]
        @test csol.q[2, end] == hsol.p[1, end]
    end

    for s in 1:4
        code = convert(PODE, hode)
        csol = integrate(code, PartitionedTableauGauss(s), Δt, nt)
        hsol = integrate(hode, PartitionedTableauGauss(s), Δt, nt)

        @test csol.q[1, end] == hsol.q[1, end]
        @test csol.p[1, end] == hsol.p[1, end]
    end
end


text:
function spelunk(caves; revisit_small_cave = false)
    future_caves = Dict{String, Vector{String}}()
    for (cave1, cave2) in caves
        cave2 != "start" && push!(get!(valtype(future_caves), future_caves, cave1), cave2)
        cave1 != "start" && push!(get!(valtype(future_caves), future_caves, cave2), cave1)
    end
    delete!(future_caves, "end")
    n_paths = 0
    states = [(present_cave = "start", past_caves = String[], small_cave_revisited = false)]
    while !isempty(states)
        present_cave, past_caves, small_cave_revisited = pop!(states)
        if present_cave == "end"
            n_paths += 1
            continue
        end
        for future_cave in future_caves[present_cave]
            is_small_cave_revisit = is_small(future_cave) && future_cave in past_caves
            small_cave_revisit_allowed = revisit_small_cave && !small_cave_revisited
            if is_small_cave_revisit && !small_cave_revisit_allowed
                continue
            end
            state = (
                present_cave = future_cave,
                past_caves = push!(copy(past_caves), present_cave),
                small_cave_revisited = small_cave_revisited || is_small_cave_revisit
            )
            push!(states, state)
        end
    end
    n_paths
end

function is_small(cave)
    all(islowercase, collect(cave))
end

caves = split.(readlines("day12/caves.txt"), "-")

p1 = spelunk(caves)

@show p1

p2 = spelunk(caves, revisit_small_cave = true)

@show p2

ext: jl


text:
### A Pluto.jl notebook ###
# v0.14.3

using Markdown
using InteractiveUtils

# ╔═╡ 7a63ed47-abf4-4432-bdce-23b7375dc02c
using Pkg, DrWatson

# ╔═╡ fc142be9-b765-4930-8d70-00a9cfda6110
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ 44feedf7-0c4a-4f17-9384-8cdb3d30f743
begin
	a = [10. 1. -5.; -20. 3. 20.; 5. 3. 5.]
	b = [1., 2., 6.]
end;

# ╔═╡ 7813d567-6e63-4056-95e3-0b0872a0d21a
# Convert to upper triang;e
function to_upper(a::Matrix, b::Vector)
  at = deepcopy(a)
  bt = deepcopy(b)
  r, c = size(at)

  for k in 1:r-1
    if abs(at[k, k]) > 1.0e-6
    for i in k+1:r
      x = at[i, k] / at[k, k]
      at[i, k] = 0.0
      for j in k+1:r
      at[i, j] -= at[k, j] * x
      end
      bt[i] -= bt[k] * x
    end
    else
    md"Zero pivot found in row $k"
    end
  end
  (a=at, b=bt)
end

# ╔═╡ 67af792e-5577-45ed-a1b5-1e977180de34
# Direct solve upper triangular set of equations
function direct_solve(a::Matrix, b::Vector)
  at = deepcopy(a)
  bt = deepcopy(b)
  r, c = size(at)

  for i in r:-1:1
    x = bt[i]
    if i < r
    for j in i+1:r
      x -= at[i, j] * bt[j]
    end
    end
    bt[i] = x / at[i, i]
  end
  bt
end

# ╔═╡ 08b41bf9-e8d6-4773-a173-605de915441d
md" ##### Modified lhs matrix and rhs"

# ╔═╡ 0d465be8-722f-4f51-ba82-f9a41915dc4e
at, bt = to_upper(a, b)

# ╔═╡ 8ede85a1-1270-423a-86a1-58ef53f33ca6
md" ##### Solution vector:"

# ╔═╡ b817e219-d7f7-4de3-a0d6-691f3f33827b
c_nmfe = direct_solve(at, bt)

# ╔═╡ 6a048f48-eadc-4657-bd01-042f43cac5e3
c = a\b

# ╔═╡ 96afae01-56d3-4b33-8cf6-9cf9ef9490c5
@test round.(c; digits=14) == round.(c_nmfe, digits=14)

# ╔═╡ 2e4ee9a0-85e1-4933-97a5-4e15ffbcafea
@test round.(a * c; digits=14) == b

# ╔═╡ Cell order:
# ╠═7a63ed47-abf4-4432-bdce-23b7375dc02c
# ╠═fc142be9-b765-4930-8d70-00a9cfda6110
# ╠═44feedf7-0c4a-4f17-9384-8cdb3d30f743
# ╠═7813d567-6e63-4056-95e3-0b0872a0d21a
# ╠═67af792e-5577-45ed-a1b5-1e977180de34
# ╟─08b41bf9-e8d6-4773-a173-605de915441d
# ╠═0d465be8-722f-4f51-ba82-f9a41915dc4e
# ╟─8ede85a1-1270-423a-86a1-58ef53f33ca6
# ╠═b817e219-d7f7-4de3-a0d6-691f3f33827b
# ╠═6a048f48-eadc-4657-bd01-042f43cac5e3
# ╠═96afae01-56d3-4b33-8cf6-9cf9ef9490c5
# ╠═2e4ee9a0-85e1-4933-97a5-4e15ffbcafea


text:
@testset "Dummy DAQ device" begin

end


text:
using LinearAlgebra
using MAT
using STFR
using Test

@testset "STFR.jl" begin
    include("signalmodel.jl")
    include("blochsim.jl")
end


text:
#structure that contains outputting info
struct OutputParams{T <: Real, R <: Real, O}
    outputs::O         #tuple of entries defining names and quantities to be outputted
    output_freq::T     #output time
    n_iter_out::R      #number of steps per output
    output_format::String     #specify output format [mat/jld]
    prefix::String     #file prefix
    output_path::String #folder in which to save
    dump_vel::Bool     #toggle on dumping the velocity after the final timestep
    zip_format::String     #specify whether or not to zip the output, and the format
    output_start::Bool  #flag to specify whether to output the initial state or not
end

#output constructor
function OutputParams(;
    outputs = (),
    output_freq = Inf,
    output_format = "jld2",
    prefix = "outfile",
    output_path = "./",
    dump_vel = false,
    zip_format = "none",
    output_start = false)

    #default the n_iter_out to -1 (this is updated in simulation once we know timestep from timestepping_params)
    n_iter_out = -1

    #check output_freq
    ((output_freq == Inf) || (output_freq > 0)) || throw(ArgumentError("output frequency must be positive or Inf"))

    #if you don't find folder, set it to the working directory
    if ~isdir(output_path)
        @warn string("Did not find output path ", output_path, ". Any outputs will go to the working directory", pwd())
        output_path = "./"
    end

    #append a "/" to folder if it doesn't have one
    endswith(output_path, "/") || (output_path = string(output_path, "/"))

    #throw an error if we don't get an output format we know
    ((output_format == "jld2") || (output_format == "mat")) || throw(ArgumentError("Output format must be `jld2` or `mat`"))

    #revert the zip format to none if we don't recognise format
    if ~(zip_format in ["none", "nc"])
        println("detected a zip format other than none or nc...
        WAVI currently only supports zipping to nc.
        Reverting to no zipping")
        zip_format = "none"
    end

    return OutputParams(outputs, output_freq, n_iter_out, output_format, prefix, output_path, dump_vel, zip_format, output_start)
end

include("output_writing.jl")
include("zipping_output.jl")


text:
using Documenter, Weave, GitHubAPI

for file ∈ readdir(joinpath(dirname(pathof(GitHubAPI)), "..", "docs", "jmd"))
      weave(joinpath(dirname(pathof(GitHubAPI)), "..", "docs", "jmd", file),
            out_path = joinpath(dirname(pathof(GitHubAPI)), "..", "docs", "src"),
            doctype = "github")
end

makedocs(format = Documenter.HTML(assets = ["assets/custom.css"]),
         modules = [GitHubAPI],
         sitename = "GitHubAPI.jl",
         pages = ["Introduction" => "index.md",
                  "Getting Started" => "getting_started.md",
                  "Public API" => "public_api.md",
                  "Private API" => "private_api.md"]
    )

deploydocs(repo = "github.com/uva-bi-sdad/GitHubAPI.jl.git")


text:
module TestSugar

using Atomix: @atomic, @atomicreplace, @atomicswap
using Test

function test_get()
    A = [42]
    @test (@atomic A[1]) === 42
    @test (@atomic A[end]) === 42
end

function test_set()
    A = [42, 43]
    @atomic A[1] = 123
    @atomic A[end] = 124
    @test A[1] === 123
    @test A[end] === 124
end

function test_inc()
    A = [1, 1]
    @test (@atomic A[1] += 123) === 124
    @test A[1] === 124
    @test (@atomic A[end] += 123) === 124
    @test A[end] === 124
end

function test_swap()
    A = [1, 1]
    @test (@atomicswap A[1] = 123) === 1
    @test A[1] === 123
    @test (@atomicswap A[end] = 456) === 1
    @test A[end] === 456
end

function test_cas()
    A = [1, 1]
    @test (@atomicreplace A[1] 1 => 123) == (old = 1, success = true)
    @test A[1] === 123
    @test (@atomicreplace A[1] 1 => 456) == (old = 123, success = false)
    @test A[1] === 123
    @test (@atomicreplace A[end] 1 => 789) == (old = 1, success = true)
    @test A[end] === 789
end

end  # module


text:
using ..ChemistryFeaturization.AbstractType: AbstractAtoms, AbstractFeaturization

"""
    FeaturizedAtoms

Container object for an [Atoms](@ref atoms) object, a [Featurization](@ref fzn), and the resulting
`encoded_features` from applying the `featurization` to the `atoms`.

## Fields
- `atoms`: [Atoms](@ref atoms) object to be featurized
- `featurization`: [Featurization](@ref fzn) scheme meant to be used for featurizing `atoms`
- `encoded_features`: The result of featurizing `atoms` using `featurization`

!!! note
    `encoded_features` will NOT change for a given atoms-featurization pair.
"""
struct FeaturizedAtoms{A<:AbstractAtoms,F<:AbstractFeaturization}
    atoms::A
    featurization::F
    encoded_features::Any
    FeaturizedAtoms{A,F}(
        atoms,
        featurization,
    ) where {A<:AbstractAtoms,F<:AbstractFeaturization} =
        new(atoms, featurization, encode(featurization, atoms))
end

FeaturizedAtoms(
    atoms::A,
    featurization::F,
) where {A<:AbstractAtoms,F<:AbstractFeaturization} =
    FeaturizedAtoms{A,F}(atoms, featurization)

# pretty printing, short
function Base.show(io::IO, fa::FeaturizedAtoms)
    a, b = size(fa.encoded_features)
    st = string(typeof(fa), " with $a x $b encoded features")
    print(io, st)
end

# pretty printing, long
function Base.show(io::IO, ::MIME"text/plain", fa::FeaturizedAtoms)
    a, b = size(fa.encoded_features)
    st = string(typeof(fa), " with $a x $b encoded features:\n")
    st = string(st, "\tAtoms: $(fa.atoms)\n\tFeaturization: $(fa.featurization)")
    print(io, st)
end

"""
    decode(featurized_atoms::FeaturizedAtoms)

Decode a [FeaturizedAtoms](@ref) object, and return the decoded value.
"""
function decode(featurized_atoms::FeaturizedAtoms)
    decoded = decode(featurized_atoms.featurization, featurized_atoms.encoded_features)
    # is this loop universally applicable?
    # for (k, v) in decoded
    #     v["Symbol"] = sef.atoms.elements[k]
    # end
    return decoded
end

"""
    featurize(atoms::AbstractAtoms, featurization::AbstractFeaturization)

Featurize an `atoms` object using a `featurization` and return the
[FeaturizedAtoms](@ref) object created.
"""
featurize(atoms::AbstractAtoms, featurization::AbstractFeaturization) =
    FeaturizedAtoms(atoms, featurization)


text:
"""
Noop ranker, does not rank, returns the first input argument unchanged.
"""
noop_ranker(idxs, scores, args...; kwargs...) = begin
    @debug "Noop ranker, pass through..."
    idxs, scores
end


text:
# ROCm-specific runtime libraries

## GPU runtime library

# reset the runtime cache from global scope, so that any change triggers recompilation
GPUCompiler.reset_runtime()

# load or build the runtime for the most likely compilation job given a compute capability
function load_runtime(dev_isa::String)
    target = GCNCompilerTarget(; dev_isa=dev_isa)
    dummy_source = FunctionSpec(()->return, Tuple{})
    params = ROCCompilerParams()
    job = CompilerJob(target, dummy_source, params)
    GPUCompiler.load_runtime(job)
end

function signal_exception()
    flag_ptr = get_global_pointer(Val(:__global_exception_flag), Int64)
    unsafe_store!(flag_ptr, 1)
    # TODO: Can we sync here?
    endpgm() # stop this wavefront
    return
end

function report_exception(ex)
    # Add kernel ID and exception string to exception ring buffer
    ring_ptr = get_global_pointer(Val(:__global_exception_ring), DevicePtr{ExceptionEntry,AS.Global})
    ring_ptr = unsafe_load(ring_ptr)
    our_signal = _completion_signal()
    prev = UInt64(1)
    while prev != UInt64(0)
        # Try to write to this slot, and skip if we fail (because another wavefront wrote first)
        prev = atomic_cas!(convert(DevicePtr{UInt64,AS.Global}, ring_ptr), UInt64(0), our_signal)
        if prev == UInt64(0)
            # We get a ReadOnlyMemoryError without this copy because the data is pinned to the device
            # TODO: Don't use an expensive host malloc
            ex_len = string_length(ex)
            copy_buf = convert(DevicePtr{UInt8,AS.Global}, malloc(ex_len))
            memcpy!(copy_buf, DevicePtr{UInt8,AS.Global}(ex), ex_len)
            unsafe_store!(convert(DevicePtr{UInt64}, ring_ptr)+sizeof(UInt64), UInt64(copy_buf))
            break
        elseif prev == UInt64(1)
            # Tail slot, give up
            break
        end
        ring_ptr += sizeof(ExceptionEntry)
    end
    return
end

# FIXME: report_oom(sz) = @rocprintf("ERROR: Out of dynamic GPU memory (trying to allocate %i bytes)\n", sz)
report_oom(sz) = @rocprintln("ERROR: Out of dynamic GPU memory")

function report_exception_name(ex)
    #= FIXME
    @rocprintf("""
        ERROR: a %s was thrown during kernel execution.
        Stacktrace:
        """, ex)
    =#
    @rocprint("""
        ERROR: an exception was thrown during kernel execution.
        Stacktrace:
        """)
    return
end

function report_exception_frame(idx, func, file, line)
    # FIXME: @rocprintf(" [%i] %s at %s:%i\n", idx, func, file, line)
    #@rocprintln(" [%i] %s at %s:%i")
    return
end

## ROCm device library

const libcache = Dict{String, LLVM.Module}()

function load_device_libs(dev_isa)
    device_libs_path === nothing && return

    isa_short = replace(dev_isa, "gfx"=>"")
    device_libs = LLVM.Module[]
    bitcode_files = (
        "hc.amdgcn.bc",
        "hip.amdgcn.bc",
        "irif.amdgcn.bc",
        "ockl.amdgcn.bc",
        "oclc_isa_version_$isa_short.amdgcn.bc",
        "opencl.amdgcn.bc",
        "ocml.amdgcn.bc",
    )

    JuliaContext() do ctx
        for file in bitcode_files
            ispath(joinpath(device_libs_path, file)) || continue
            name, ext = splitext(file)
            lib = get!(libcache, name) do
                file_path = joinpath(device_libs_path, file)
                parse(LLVM.Module, read(file_path), ctx)
            end
            push!(device_libs, lib)
        end
    end

    @assert !isempty(device_libs) "No device libs detected!"
    return device_libs
end

function link_device_libs!(mod::LLVM.Module, dev_isa::String, undefined_fns)
    libs::Vector{LLVM.Module} = load_device_libs(dev_isa)

    ufns = undefined_fns
    # TODO: only link if used
    # TODO: make these globally/locally configurable
    ctx = LLVM.context(mod)
    link_oclc_defaults!(mod, dev_isa, ctx)
    for lib in libs
        # override libdevice's triple and datalayout to avoid warnings
        triple!(lib, triple(mod))
        datalayout!(lib, datalayout(mod))
    end

    GPUCompiler.link_library!(mod, libs)
end

function link_oclc_defaults!(mod::LLVM.Module, dev_isa::String, ctx; finite_only=false,
                             unsafe_math=false, correctly_rounded_sqrt=true, daz=false)
    # link in some defaults for OCLC knobs, to prevent undefined variable errors
    lib = LLVM.Module("OCLC")
    triple!(lib, triple(mod))
    datalayout!(lib, datalayout(mod))

    isa_short = replace(dev_isa, "gfx"=>"")
    for (name,value) in (
            "__oclc_ISA_version"=>parse(Int32, isa_short),
            "__oclc_finite_only_opt"=>Int32(finite_only),
            "__oclc_unsafe_math_opt"=>Int32(unsafe_math),
            "__oclc_correctly_rounded_sqrt32"=>Int32(correctly_rounded_sqrt),
            "__oclc_daz_opt"=>Int32(daz))
        gvtype = convert(LLVMType, typeof(value), ctx)
        gv = GlobalVariable(lib, gvtype, name, 4)
        init = ConstantInt(Int32(0), ctx)
        initializer!(gv, init)
        unnamed_addr!(gv, true)
        constant!(gv, true)
    end

    link!(mod, lib)
end


text:
function load(model::JuMP.Model, d::Data)
    d.state = Loading
    for variable in d.variables
        load(model, variable)
    end
    for (index, constraint) in d.constraints
        cref = load(model, constraint)
        if cref !== nothing
            d.transformed_constraints[index] = cref
        end
    end
    if d.objective !== nothing
        load(model, d.objective)
    end
    d.state = Modeling
end

# set as `optimize_hook` to JuMP Model in `data` so it is called at
# `JuMP.optimize!` if at least one set variable is created
function optimize_hook(model::JuMP.AbstractModel)
    d = data(model)
    set_space(d, model)
    # In case `optimize!` is called then the problem is modified and then it is
    # called again we need to clear first the space that might be wrong
    synchronize_perspective(d)
    clear_spaces(d)
    create_spaces(d)
    load(model, d)
    JuMP.optimize!(model, ignore_optimize_hook = true)
end


text:
module AAA_RNG

const cm1 = 259200
const cia1 = 7141
const cic1 = 54773
const crm1 = 3.8580247e-6

const cm2 = 134456
const cia2 = 8121
const cic2 = 28411
const crm2 = 7.4373773e-6

const cm3 = 243000
const cia3 = 4561
const cic3 = 51349

mutable struct AAARNG_State
    state::Vector{Float64}
    IX1::UInt64
    IX2::UInt64
    IX3::UInt64
end

function reseed(nSeed::UInt64, state::AAARNG_State)

    mIX1 = (cic1 - nSeed)
    mIX1 = mod(cia1*mIX1 + cic1, cm1)
    mIX2 = mod(mIX1, cm2)
    mIX1 = mod(cia1*mIX1 + cic1, cm1)
    #mIX2 = mod(cia2*mIX2 + cic2, cm2)
    mIX3 = mod(mIX1, cm3)

    for i=1:97
        mIX1 = mod(cia1*mIX1 + cic1, cm1)
        mIX2 = mod(cia2*mIX2 + cic2, cm2)
        state.state[i] = (mIX1 + (mIX2)*crm2)*crm1
    end

    return AAARNG_State(state.state, mIX1, mIX2, mIX3)

end

mutable struct AAARNG <: AbstractRNG
    seed::UInt64
    state::AAARNG_State
    function AAARNG(seed::Int)
        new(unsigned(seed), reseed(unsigned(seed), AAARNG_State(zeros(97), 0, 0, 0)))
    end
end

const GLOBAL_AAARNG = AAARNG(6499)

function rand(rng::AAARNG=GLOBAL_AAARNG)
    mIX1 = rng.state.IX1
    mIX2 = rng.state.IX2
    mIX3 = rng.state.IX3
    mIX1 = mod(cia1*(mIX1) + cic1, cm1)
    mIX2 = mod(cia2*(mIX2) + cic2, cm2)
    mIX3 = mod(cia3*(mIX3) + cic3, cm3)

    #println(typeof(mIX3))
    #println(mIX3, ", ", 1+Int(round(97*Float64(mIX3)/cm3, RoundDown)))
    j = 1+Int(floor(97*Float64(mIX3)/cm3)) #1 + Int(97*Float64(mIX3)/cm3)

    retval::Float64 = rng.state.state[j]

    u::Float64 = (mIX1 + (mIX2)*crm2)*crm1
    if u<=0.0 || u>=1.0
        u = 1.0 + floor(u) - u
    end

    rng.state.state[j] = u
    rng.state.IX1 = mIX1
    rng.state.IX2 = mIX2
    rng.state.IX3 = mIX3

    return retval
end

end # module


text:
abstract RustModel
    # beta::Float64
    # u::Array{Float64,2}                 # u[k,a]: flow utilities
    # v::Array{Float64,2}                 # v[k,a]: choice-specific discounted value
    # V::Array{Float64,1}                 # V[k]: expected value
    # pi::Array{Float64,3}                # pi[k,k',a]: condition state transitions
    # p::Array{Float64,2}                 # p[k,a]: conditional choice probabilities (CCP)

function solvedp!(model::RustModel)
    da=size(model.u)[2]

    function ff!(z,fv)
        fv[:]=-1.0
        for ia=1:da
            # fv[:]+=exp.(min(model.u[:,ia]+model.beta*model.pi[:,:,ia]*z-z,1))
            fv[:]+=exp.(model.u[:,ia]+model.beta*model.pi[:,:,ia]*z-z)
        end
    end

    function fj!(z,fm)
        fm[:]=0.0                           #inefficient reallocation of temp at each call
        for ia=1:da
            exponent=model.u[:,ia]+model.beta*model.pi[:,:,ia]*z-z
            # fm[:,:]+=diagm(1(exponent.<1).*exp.(exponent))*(model.beta*model.pi[:,:,ia]-eye(size(fm,1)))
            fm[:,:]+=diagm(exp.(exponent))*(model.beta*model.pi[:,:,ia]-eye(size(fm,1)))
        end
    end

    df = NLsolve.DifferentiableMultivariateFunction(ff!,fj!)

    res=NLsolve.nlsolve(df,model.V,autoscale=false) #will automatically start at previous value

    # println(res)
    model.V[:]=res.zero
    for ia=1:da
        model.p[:,ia]=exp(model.u[:,ia]+model.beta*model.pi[:,:,ia]*model.V-model.V)
    end
    println(checkdp(model), " DP solved?")
    # println(round(model.p,2))
end

function checkdp(model::RustModel)
    V1=\(eye(model.pi[:,:,1])-model.beta*model.pi[:,:,1],model.u[:,1]-log(model.p[:,1]))
    V2=\(eye(model.pi[:,:,2])-model.beta*model.pi[:,:,2],model.u[:,2]-log(model.p[:,2]))
    norm(V1-V2)
end

########## Identification

# "nonparametric" two-step projection is an implicit system of d_k(d_a+1) equations Phi((u,V),Pi)=A(Pi) (u,V) - B(P)
# to compute the Jacobian of eta ---> u(Pi(eta)) we use the implicit function theorem on the full X=(u,V)
# the implicit function theorem require inverting jx:=dPhi/dX. Fortunately this is closed form.
# A(Pi)=[A 0; I B] ---> A(Pi)^{-1} = [BC I-BCA;-C CA], C=(AB)^{-1}
# After three steps:
# (1) selecting the u block that we are interested in
# (2) using the fact the identifiying equation does not depend on Pi
# (3) developing dA(Pi)dPi using kronecker product formulas
# we find: du/dPi = beta (I - B(AB)^{-1}A) (I_{d_a} \otimes V' \otimes I_{d_k}) J_Pi
# where J_Pi = d (vec(Pi))/d alpha is (da dk dk x deta) is the Jacobian of Pi(eta) can be computed numerically
# the computational bottleneck is inverting (AB) which is (dk x dk)

function explicitinv(a,b)
    n=size(b,1)
    c=inv(a*b)
    [b*c eye(n)-b*c*a; -c c*a]
end

function idrank(model::RustModel,a,eta)
    dk,da=size(model.u)
    b=[1(ind2sub((dk,da),ika)[1]==jk)-model.beta*model.pi[ind2sub((dk,da),ika)[1],jk,ind2sub((dk,da),ika)[2]] for ika=1:dk*da,jk=1:dk]
    m=sparse(I-b*((a*b)\a))
    function ff(eta1)
        pi!(model,eta1)
        piv=vcat([model.pi[:,:,ia]*model.V for ia=1:da]...)
        model.beta*m*piv
    end
    jpi=Calculus.finite_difference_jacobian(ff,eta)
    rjpi=rank(jpi)
    println("rank of eta --> u is $rjpi")
    println("(unconstrained transition matrices are of dimension $(dk*(dk-1)*da))")
    println("ambient space is of dimension $(dk*(da-1))")
    println("expect dimension of partial identified set to be dtheta-$(dk*(da-1)-rjpi)")
end

########## Simulation

function rand(model::RustModel,T)
    dk,da=size(model.p)
    ka=zeros(Int,T,2)
    ka[1,1]=1
    ka[1,2]=wsample(1:da,model.p[ka[1,1],:])
    for t=1:T-1
        ka[t+1,1]=wsample(1:dk,model.pi[ka[t,1],:,ka[t,2]])
        ka[t+1,2]=wsample(1:da,model.p[ka[t+1,1],:])
    end
    ka
end


text:
"""
    TOMLConfig

Use TOML files to configure command line parsing via [ArgParse.jl](https://github.com/carlobaldassi/ArgParse.jl).
"""
module TOMLConfig

using AbstractTrees
using Dates
using Reexport

@reexport using ArgParse, TOML

export Config

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
    @eval Base.Experimental.@compiler_options compile=min optimize=0 infer=false
end

@nospecialize # use only declared type signatures, helps with compile time

"""
    Config(toml::AbstractDict{String})
    Config(; filename::String)

Basic tree structure for navigating TOML file contents.
Each node in the `Config` tree represents a single section of a TOML file.
Children of a `Config` node are the corresponding TOML subsections, if they exist.

# Examples

```jldoctest
julia> using TOMLConfig

julia> cfg = Config(TOML.parse(
       \"\"\"
       a = 1

       [sec1]
           b = 2

           [sec1.sub1]
           c = 3
       \"\"\"))
TOML Config with contents:

a = 1

[sec1]
b = 2

    [sec1.sub1]
    c = 3

julia> cfg.sec1
TOML Config with contents:

b = 2

[sub1]
c = 3

julia> cfg.sec1.sub1
TOML Config with contents:

c = 3
```
"""
struct Config <: AbstractDict{String, Any}
    "TOML section contents"
    contents::AbstractDict{String, Any}

    "Config node corresponding to parent TOML section, or `nothing` for the root node"
    parent::Union{Config, Nothing}

    "Key within the parent node which points to this node, or `nothing` for the root node"
    key::Union{String, Nothing}
end
Config(toml::AbstractDict{String}) = Config(convert(Dict{String, Any}, toml), nothing, nothing)
Config(; filename::String) = Config(TOML.parsefile(filename))

# Forward `AbstractDict`s methods to `Config.contents`
for f in [:(Base.get), :(Base.get!), :(Base.getkey), :(Base.haskey), :(Base.iterate), :(Base.keys), :(Base.length), :(Base.pop!), :(Base.values)]
    @eval $f(cfg::Config, args...; kwargs...) = $f(contents(cfg), args...; kwargs...)
end
for f in [:(Base.delete!), :(Base.pairs), :(Base.sizehint!)]
    @eval $f(cfg::Config, args...; kwargs...) = Config($f(contents(cfg), args...; kwargs...))
end
for f in [:(Base.merge), :(Base.merge!)]
    @eval $f(cfg::Config, others::Config...; kwargs...) = Config($f(contents(cfg), map(contents, others)...; kwargs...))
end

# Define getters to access struct fields, since `getproperty` is overloaded for convenience below
contents(cfg::Config) = getfield(cfg, :contents)
parent(cfg::Config) = getfield(cfg, :parent)
key(cfg::Config) = getfield(cfg, :key)

function Base.getproperty(cfg::Config, k::Symbol)
    v = contents(cfg)[String(k)]
    if v isa AbstractDict
        v = Config(v, cfg, String(k))
    end
    return v
end

Base.setproperty!(cfg::Config, k::Symbol, v) = contents(cfg)[String(k)] = v

AbstractTrees.nodetype(::Config) = Config
AbstractTrees.children(parent::Config) = [Config(contents, parent, key) for (key, contents) in contents(parent) if is_child(contents)]

function AbstractTrees.printnode(io::IO, cfg::Config)
    if key(cfg) !== nothing
        println(io, key(cfg) * ":")
    end
    print(io, join(["$k = $v" for (k,v) in contents(cfg) if !(v isa AbstractDict)], "\n"))
end

function Base.show(io::IO, ::MIME"text/plain", cfg::Config)
    println(io, "TOML Config with contents:\n")
    TOML.print(io, contents(cfg))
end

# Convenience functions for classifying nodes
is_child(v) = !is_arg(v)
is_arg(v) = is_plain_arg(v) || is_dict_arg(v)
is_plain_arg(v) = !(v isa AbstractDict)
is_dict_arg(v) = v isa AbstractDict && arg_key() ∈ keys(v)
arg_props(v::AbstractDict{String}) = recurse_convert_keytype(delete!(deepcopy(v), arg_key()), Symbol)
arg_value(v) = is_dict_arg(v) ? deepcopy(v[arg_key()]) : deepcopy(v)
arg_value(d::AbstractDict{String}, k::String) = arg_value(d[k])
arg_value!(d::AbstractDict{String}, v, k::String) = is_dict_arg(d[k]) ? (d[k][arg_key()] = deepcopy(v)) : (d[k] = deepcopy(v))

# Convenience functions for getting/setting deeply nested nodes
recurse_getindex(d::AbstractDict, keys) = foldl((dᵢ,k) -> dᵢ[k], keys; init = d)
recurse_setindex!(d::AbstractDict, v, keys) = recurse_getindex(d, keys[begin:end-1])[keys[end]] = v
recurse_convert_keytype(d::AbstractDict, ::Type{K} = Symbol) where {K} = Dict{K, Any}(K(k) => v isa AbstractDict ? recurse_convert_keytype(v, K) : v for (k,v) in d)
recurse_convert_valtype(d::AbstractDict, ::Type{V} = Symbol) where {V} = Dict{keytype(d), Any}(k => v isa AbstractDict ? recurse_convert_valtype(v, V) : V(v) for (k,v) in d)
recurse_convert_keyvaltype(d::AbstractDict, ::Type{K} = Symbol, ::Type{V} = Symbol) where {K, V} = Dict{K, Any}(K(k) => v isa AbstractDict ? recurse_convert_keyvaltype(v, K, V) : V(v) for (k,v) in d)

"""
    parse_config(toml::AbstractDict{String})
    parse_config(; filename::String)

Convenience method for parsing configuration files. Equivalent to `ArgParse.parse_args(Config(toml))`
and `ArgParse.parse_args(Config(; filename = filename))`, respectively.
"""
parse_config(toml::AbstractDict{String}) = ArgParse.parse_args(Config(toml))
parse_config(; filename::String) = ArgParse.parse_args(Config(; filename = filename))

"""
    save_config(cfg::Union{Config, AbstractDict}; filename::AbstractString, kwargs...)

Save configuration `cfg` to TOML file `filename`.
Additional keyword arguments are forwarded to `TOML.print`.
"""
function save_config(cfg::AbstractDict; filename::AbstractString, kwargs...)
    open(filename; write = true) do io
        TOML.print(io, cfg; kwargs...)
    end
end
save_config(cfg::Config; kwargs...) = save_config(deepcopy(Dict(cfg)); kwargs...)

# Parser settings
default_args_list() = ARGS
default_argparse_settings() = ArgParseSettings()
default_parser_settings() = Dict{String, String}(
    "arg_key"              => "_ARG_",
    "arg_required_value"   => "_REQUIRED_",
    "flag_delim"           => ".",
    "inherit_all_key"      => "_INHERIT_",
    "inherit_parent_value" => "_PARENT_",
)

const _parser_settings = default_parser_settings()
parser_settings()      = _parser_settings
parser_settings(k)     = _parser_settings[string(k)]
parser_settings!(k, v) = _parser_settings[string(k)] = string(v)

arg_key()              = parser_settings("arg_key")
arg_required_value()   = parser_settings("arg_required_value")
flag_delim()           = parser_settings("flag_delim")
inherit_all_key()      = parser_settings("inherit_all_key")
inherit_parent_value() = parser_settings("inherit_parent_value")

"""
    parser_settings!(;
        arg_key = "$(arg_key())",
        arg_required_value = "$(arg_required_value())",
        flag_delim = "$(flag_delim())",
        inherit_all_key = "$(inherit_all_key())",
        inherit_parent_value = "$(inherit_parent_value())",
    )

Customize global settings for parsing `Config` structs.
To restore default settings, use `TOMLConfig.parser_settings!()`.

# Keywords:
* `arg_key`: if this key is found in a TOML section, the rest of the section is interpreted as properties for the argument table entry
* `arg_required_value`: if this value is found in a TOML section, the `required = true` is passed to the argument table entry
* `inherit_all_key`: if this key is found in a TOML section, all fields from the immediate parent section (i.e., non-recursive) should be inherited
* `inherit_parent_value`: if this value is found in a TOML section, it is replaced with the value corresponding to the same key in the immediate parent section (i.e., non-recursive)
* `flag_delim`: command line flags for keys in nested TOML sections are formed by joining all parent keys together with this delimiter
"""
function parser_settings!(; kwargs...)
    for (k,v) in merge!(default_parser_settings(), recurse_convert_keyvaltype(kwargs, String, String))
        if k ∈ keys(parser_settings())
            parser_settings!(k, v)
        else
            error("Invalid parser setting: $(k). Possible settings are: $(join(sort(collect(keys(default_parser_settings()))), ", ")).")
        end
    end
end

"""
    defaults!(cfg::Config)

Populate fields of TOML config which are specified to have default values inherited from parent sections.

# Examples

```jldoctest
julia> using TOMLConfig

julia> cfg = TOMLConfig.defaults!(Config(TOML.parse(
       \"\"\"
       a = 1
       b = 2

       [sec1]
       b = \"$(inherit_parent_value())\"
       c = 3

           [sec1.sub1]
           $(inherit_all_key()) = \"$(inherit_parent_value())\"
       \"\"\")))
TOML Config with contents:

b = 2
a = 1

[sec1]
c = 3
b = 2

    [sec1.sub1]
    c = 3
    b = 2
```
"""
function defaults!(cfg::Config; replace_arg_dicts = false)
    _INHERIT_ = inherit_all_key()
    _PARENT_  = inherit_parent_value()

    # Step 0:
    #   Breadth-first search to replace arg table dictionaries with default values, which may be "_PARENT_"
    if replace_arg_dicts
        for node in StatelessBFS(cfg)
            for (k,v) in collect(contents(node))
                if is_dict_arg(v)
                    contents(node)[k] = arg_value(v)
                end
            end
        end
    end

    # Step 1:
    #   Inverted breadth-first search for "_INHERIT_" keys with value "_PARENT_".
    #   If found, copy all key-value pairs from the immediate parent (i.e. non-recursive) into the node containing "_INHERIT_".
    #   Delete the "_INHERIT_" key afterwards.
    for node in reverse(collect(StatelessBFS(cfg)))
        parent(node) === nothing && continue
        !haskey(contents(node), _INHERIT_) && continue
        contents(node)[_INHERIT_] != _PARENT_ && continue
        for (k,v) in collect(contents(parent(node)))
            if is_arg(v) && !haskey(contents(node), k)
                # If key `k` is not already present in the current section, inherit arg (possibly an arg dict) from the parent section
                contents(node)[k] = deepcopy(v)
            end
        end
        delete!(contents(node), _INHERIT_)
    end

    # Step 2:
    #   Breadth-first search for fields with value "_PARENT_".
    #   If found, copy default value from the corresponding field in the immediate parent (i.e. non-recursive).
    for node in StatelessBFS(cfg)
        parent(node) === nothing && continue
        for (k,v) in collect(contents(node))
            if arg_value(v) == _PARENT_
                # Set arg value to the arg value of the parent (both parent and/or child may be arg dicts)
                arg_value!(contents(node), arg_value(contents(parent(node)), k), k)
            end
        end
    end

    return cfg
end

"""
    arg_table_flag(cfg::Config, k::String)

Generate command flag corresponding to nested key `k` in a `Config` node.
The flag is constructed by joining the keys recursively from the parents
of the current node using the delimiter `flag_delim()` and prepending "--".

# Examples

Given a `Config` node with contents
```julia
a = 1
b = 2

[sec1]
c = 3

    [sec1.sub1]
    d = 4
```

The corresponding flags that will be generated are
```julia
--a
--b
--sec1$(flag_delim())c
--sec1$(flag_delim())sub1$(flag_delim())d
```
"""
function arg_table_flag(cfg::Config, k::String)
    flag = k
    while true
        if parent(cfg) === nothing
            flag = "--" * flag
            return flag
        else
            flag = key(cfg) * flag_delim() * flag
            cfg = parent(cfg)
        end
    end
end

"""
    ArgParse.add_arg_table!(settings::ArgParseSettings, cfg::Config)

Populate `settings` argument table using configuration `cfg`.

# Examples

```jldoctest
julia> using TOMLConfig

julia> cfg = Config(TOML.parse(
       \"\"\"
       a = 1.0
       b = 2

       [sec1]
       c = [3, 4]

           [sec1.sub1]
           d = "d"
       \"\"\"));

julia> settings = add_arg_table!(ArgParseSettings(prog = "myscript.jl"), cfg);

julia> ArgParse.show_help(settings; exit_when_done = false)
usage: myscript.jl [--b B] [--a A] [--sec1.c [SEC1.C...]]
                   [--sec1.sub1.d SEC1.SUB1.D]

optional arguments:
  --b B                 (type: Int64, default: 2)
  --a A                 (type: Float64, default: 1.0)
  --sec1.c [SEC1.C...]  (type: Int64, default: [3, 4])
  --sec1.sub1.d SEC1.SUB1.D
                        (default: "d")
```
"""
function ArgParse.add_arg_table!(settings::ArgParseSettings, cfg::Config)
    # Populate settings argument table
    for node in reverse(collect(PostOrderDFS(cfg)))
        for (k,v) in contents(node)
            if is_plain_arg(v)
                # Add to arg table using specified value as the default
                props = Dict{Symbol,Any}()
                props[:default] = deepcopy(v)
                props[:required] = false

                if v isa AbstractVector
                    props[:arg_type] = eltype(v)
                    props[:nargs] = '*'
                else
                    props[:arg_type] = typeof(v)
                end

                add_arg_table!(settings, arg_table_flag(node, k), props)

            elseif is_dict_arg(v)
                props, v = arg_props(v), arg_value(v)

                # Special-casing specific properties
                if :arg_type ∈ keys(props)
                    arg_type_dict = Dict{String, DataType}("Any" => Any, "DateTime" => DateTime, "Time" => Time, "Date" => Date, "Bool" => Bool, "Int" => Int, "Float64" => Float64, "String" => String)
                    if props[:arg_type] ∉ keys(arg_type_dict)
                        error("Invalid arg_type: $(repr(props[:arg_type])). Must be one of: $(join(sort(repr.(keys(arg_type_dict))), ", ")).")
                    end
                    props[:arg_type] = arg_type_dict[props[:arg_type]]
                end

                if :nargs ∈ keys(props)
                    nargs_dict = Dict{String, Char}("A" => 'A', "?" => '?', "*" => '*', "+" => '+', "R" => 'R')
                    if !(props[:nargs] isa Int) && !(props[:nargs] ∈ keys(nargs_dict))
                        error("Invalid nargs: $(repr(props[:nargs])). Must be a nonnegative integer, or one of: $(join(sort(repr.(keys(nargs_dict))), ", ")).")
                    end
                    if props[:nargs] ∈ keys(nargs_dict)
                        props[:nargs] = nargs_dict[props[:nargs]]
                    end
                end

                # Add to arg table using user specified props
                if v == arg_required_value()
                    get!(props, :required, true)
                else
                    get!(props, :default, deepcopy(v))
                    get!(props, :required, false)
                    if v isa AbstractVector
                        get!(props, :arg_type, eltype(v))
                        get!(props, :nargs, '*')
                    else
                        get!(props, :arg_type, typeof(v))
                    end
                end

                add_arg_table!(settings, arg_table_flag(node, k), props)
            end
        end
    end

    return settings
end

"""
    ArgParse.parse_args(
        [args_list::Vector = ARGS,]
        [settings::ArgParseSettings = ArgParseSettings(),]
        cfg::Config;
        as_dict = false,
        as_symbols = false,
        kwargs...
    )

Parse TOML configuration struct with command line arguments `args_list`.

# Arguments:
* `args_list::Vector`: (optional) vector of arguments to be parsed
* `settings::ArgParseSettings`: (optional) settings struct which will be configured according to `cfg`
* `cfg::Config`: TOML configuration settings

# Keywords:
* `as_dict`: if true, return config as a dictionary with `String` keys, otherwise return a `Config` struct
* `as_symbols`: if true and `as_dict=true`, return config dictionary with `Symbol` keys
* `kwargs...`: additional keyword arguments can be passed to customize the global parser settings; see `TOMLConfig.parser_settings!`. Note that these settings remain persistent on subsequent calls to `parse_args`; use `TOMLConfig.parser_settings!()` to restore all defaults.

# Examples

```jldoctest
julia> using TOMLConfig

julia> cfg = Config(TOML.parse(
       \"\"\"
       a = 1
       b = 2

       [sec1]
       b = \"$(inherit_parent_value())\"
       c = 3

           [sec1.sub1]
           $(inherit_all_key()) = \"$(inherit_parent_value())\"
       \"\"\"));

julia> parsed_args = parse_args(["--a", "3", "--sec1.b", "5", "--sec1.c", "10"], cfg)
TOML Config with contents:

b = 2
a = 3

[sec1]
c = 10
b = 5

    [sec1.sub1]
    c = 10
    b = 5
```
"""
function ArgParse.parse_args(args_list::Vector, settings::ArgParseSettings, cfg::Config; kwargs...)
    parse_args!(args_list, deepcopy(settings), deepcopy(cfg); kwargs...)
end
ArgParse.parse_args(cfg::Config; kwargs...) = parse_args(default_argparse_settings(), cfg; kwargs...)
ArgParse.parse_args(settings::ArgParseSettings, cfg::Config; kwargs...) = parse_args(default_args_list(), settings, cfg; kwargs...)
ArgParse.parse_args(args_list::Vector, cfg::Config; kwargs...) = parse_args(args_list, default_argparse_settings(), cfg; kwargs...)

function parse_args!(
        args_list::Vector,
        settings::ArgParseSettings,
        cfg::Config;
        as_dict::Bool = false,
        as_symbols::Bool = false,
        kwargs...
    )
    # Update global parser settings
    parser_settings!(; kwargs...)

    # Populate all "_INHERIT_" keys and/or "_PARENT_" values to establish defaults
    default_cfg = defaults!(deepcopy(cfg); replace_arg_dicts = false)

    # Add arg table entries to `settings` dynamically using default `cfg` specification
    add_arg_table!(settings, default_cfg)

    # Parse `args_list` arguments into `cfg` using the dynamically constructed arg table
    for (k,v) in parse_args(args_list, settings)
        if any(startswith("--" * k), args_list)
            # Only update `cfg` with new value if it was explicitly passed in `args_list`
            keys = String.(split(k, flag_delim()))
            recurse_setindex!(contents(cfg), deepcopy(v), keys)
        end
    end

    # Populate remaining "_INHERIT_" keys and/or "_PARENT_" values which were not overridden by `args_list`
    defaults!(cfg; replace_arg_dicts = true)

    # Return parsed config
    if as_dict
        if as_symbols
            return recurse_convert_keytype(contents(cfg), Symbol)
        else
            return contents(cfg)
        end
    else
        return cfg
    end
end

end # module TOMLConfig

max_stars_repo_name: jondeuce/TOMLConfig

max_issues_repo_name: jondeuce/TOMLConfig

max_forks_repo_name: jondeuce/TOMLConfig


text:
#!/usr/bin/env julia

import Parsers

const DATA_DIR   = "../13_data"
const INPUT_FILE = "input.txt"

# ----------------------------------------
macro assert(bool_expr)
    message = string("Assertion: ", bool_expr, " failed")
    :($(esc(bool_expr)) || error($message))
end
# ----------------------------------------

function banner(txt::String)
    println('='^10, ' ', txt, ' ', '='^10)
end

function get_earliest_ts(freq_str)
    strs = split(strip(freq_str), ",")
    str2int = s->Parsers.tryparse(Int, string(s))
    bus_freqs = [ (eid[1]-1,eid[2]) for eid in enumerate(str2int.(strs))
                  if !isnothing(eid[2])
                ]
    println("bus_freqs=$(bus_freqs)")
    bus_freqs = sort(bus_freqs, by=eid->eid[2], rev=true)
    println("bus_freqs=$(bus_freqs)")

    ts = 0
    step = 1
    for (offset, id) in bus_freqs
        while true
            if (ts + offset) % id == 0
               step *= id
               break
            end
            ts += step
        end
    end
    return ts
end

# ----------------------------------------

function test()
    test_freqs("17,13,x,x,59,x,31,19", 1068781)

    test_freqs("17,x,13,19",      3417)
    test_freqs("67,7,59,61",      754018)
    test_freqs("67,x,7,59,61",    779210)
    test_freqs("67,7,x,59,61",    1261476)
    test_freqs("1789,37,47,1889", 1202161486)
end

function test_freqs(freq_str, expected)
    result = get_earliest_ts(freq_str)
    if (result == expected)
        println("\tPASS: $(result) == $(expected)")
    else
        println("\tFAIL: $(result) != $(expected)")
    end
end

function test_maxind(starting, max_ind, expected)
    println("====================")
    result = get_maxind(starting, max_ind)
    if (result == expected)
        println("Test passed: $(starting) @ $(max_ind) == $(expected)")
    else
        println("TEST FAILED: $(starting) @ $(max_ind) != $(expected)")
    end
end

# ----------------------------------------

function part1(lines)
    banner("PART 1")
    ready_time = parse(Int, lines[1])
    strs = split(strip(lines[2]), ",")
    str2int = s->Parsers.tryparse(Int, string(s))
    bus_freqs = [k for k in str2int.(strs) if !isnothing(k)]
    wait_times = [(bf, bf - (ready_time % bf)) for bf in bus_freqs]
    bus_wait = sort(wait_times, by=((bf,wt),)->wt)[1]

    println("ready_time=$(ready_time)")
    println("bus_freqs=$(bus_freqs)")
    println("wait_times=$(wait_times)")
    # println("bus_target=$(bus_target)")
    println("solution=$(bus_wait[1] * bus_wait[2])")
end

function part2(lines)
    println()
    banner("PART 2")
    result = get_earliest_ts(lines[2])
    println("result=$(result)")
end

function main()
    filename = "$DATA_DIR/$INPUT_FILE"
    lines::Array{String} = readlines(filename)
    part1(lines)
    part2(lines)
end

main()

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
"""
main : the code you run to double-encode... the "main loop" of optimization.
"""

module main

include("utils.jl")
include("types.jl")
include("lookup.jl")
include("std_setup.jl")
include("math.jl")
include("optimize.jl")
include("mrf.jl")

using ArgParse, Logging, JLD, StatsBase, Distributions, Random, Statistics

function convert_to_saveable(cur_pop)
	new_pop = types.SaveChrome[]

	for indiv in cur_pop
		push!(new_pop, types.SaveChrome(indiv.path, indiv.full_sequence,
										indiv.deg_nuc, indiv.deg_map, indiv.deg_trns, indiv.deg_trne, indiv.deg_d, indiv.deg_skip, indiv.deg_insert,
										indiv.mark_nuc, indiv.mark_map, indiv.mark_trns, indiv.mark_trne, indiv.mark_d, indiv.mark_skip, indiv.mark_insert,
										indiv.deg_prob, indiv.deg_base_E, indiv.deg_seq,
										indiv.mark_prob, indiv.mark_base_E, indiv.mark_seq, indiv.first_weight))
	end

	return new_pop
end

"""
set_up_and_optimize is the code that sets up appropriate information about genes we want to target for double-encoding and
orchestrates the initial HMM double-encoding solution and then iteratively improves on that with greedy pseudolikelihood improvements.
Code also logs information through Logging package.
Genes are often referred to as mark/deg. This is historic from times of thinking of "marker" genes and "designated essential gene".
At this point this distinction is not meaningful but is kept to avoid equally arbitrary x_name, y_name conventions.
Arguments:
	mark_name : gene ID for 'mark' gene. Needed as a key for looking up some values associated with genes in text files.
	deg_name : gene ID for 'deg' gene.
	mark/deg_grem : path to MRF parameter files, usually a JLD/CSV file.
	mark/deg_hmm : path to HMM files, usually .hmm files.
	pop_size : size of population, i.e. number of individual HMM solutions to greedily optimize.
	frame : p1/p2.
	max_iter : maximum number of iterations.
	X_range/Y_range: positions along gene to consider for double-encoding... useful if you want to restrict where double-encoding can occur.
						Ranges are defined in terms of the end of the sequence, so if you have a 70 aa sequence you want to start at position 80,
						set the range to be something like 150:151. Generally some buffer around this value (i.e. 150:160) is useful.
						Don't set, or set to false if you don't care where the genes overlap.
"""
function set_up_and_optimize(log_io, rand_barcode, out_path, mark_name, deg_name, mark_grem, deg_grem, mark_hmm, deg_hmm, pop_size, frame, max_iter, X_range=false, Y_range=false; rand_weights = false, actually_mrf = true)

	@debug("Beginning run.")
	@debug(Libc.strftime(time()))
	@debug(Libc.strftime(time()))

	println("The random barcode on this run is: $rand_barcode")
	@debug("The random barcode on this run is: $rand_barcode")

	do_cull = false #culling reduces number of sequences we optimize over time. Just a trick to save time if you want top-performers only.
	full_even_window = 20
	half_window = div(full_even_window, 2) #window statistics to see if scores still decreasing (sometimes used).
	last_few_sf = Float64[]
	last_few_cs = Float64[]

	#Looks up mean/std. dev of family pseudolikelihoods, pre-computed.
	mu_mark, sig_mark = lookup.mu_sig(deg_name, "psls/")
	mu_deg, sig_deg = lookup.mu_sig(mark_name, "psls/")

	@debug("Stat param vars (mark / deg):\t$(mu_mark)\t$(sig_mark)\t$(mu_deg)\t$(sig_deg)")

	if true #replacing try/catch.
		gen_samples = true #this is true if starting from scratch (general case), false if we have some candidates to optimize ahead of time.
		population = types.Chromosome[]
		local mark_grem_prot, deg_grem_prot
		if !gen_samples
			population = load("trpE_population.jld")["population"][1:100] #[1:100]
			mark_gremodel, deg_gremodel = std_setup.short_set_up(mark_name, deg_name, mark_grem, deg_grem, mark_hmm, deg_hmm, length(population), rand_barcode, frame)
			deg_nNodes, mark_nNodes = deg_gremodel.nNodes, mark_gremodel.nNodes
		else
			#Generate population of sampled hmm starting points.
			if X_range == false || Y_range == false #we require both to be there...
				@debug("Doing standard full set up...")
				mark_gremodel, deg_gremodel, population, mark_grem_prot, deg_grem_prot = std_setup.full_set_up(out_path, mark_name, mark_hmm, mark_grem, deg_name, deg_hmm, deg_grem, pop_size, 1200, 1200, rand_barcode, frame)
			else
				@debug("The x range is $X_range")
				@debug("The y range is $Y_range")
				mark_gremodel, deg_gremodel, population, mark_grem_prot, deg_grem_prot = std_setup.full_sample_set_up(out_path, mark_name, mark_hmm, mark_grem, deg_name, deg_hmm, deg_grem, pop_size, 1200, 1200, rand_barcode, X_range, Y_range, frame)
			end
			deg_nNodes, mark_nNodes = deg_gremodel.nNodes, mark_gremodel.nNodes
		end

		#Look through generated candidates, check their fitness values and their correctness.
		success_chrom, fdp, fmp, fitness_values, indie_deg_maps, indie_mark_maps, deg_ull, deg_pv_w1, deg_pv_w2, deg_prot_mat, mark_ull, mark_pv_w1, mark_pv_w2, mark_prot_mat = optimize.assess_founders(population, mark_gremodel.nNodes, mark_gremodel.w1, mark_gremodel.w2, deg_gremodel.nNodes, deg_gremodel.w1, deg_gremodel.w2)

		@debug("Num success... $(length(success_chrom))")

		if actually_mrf
			founding_fitness_values = fitness_values[1:end]

			mark_base_energy = mrf.basic_energy_calc(mark_grem_prot, mark_gremodel.w1, mark_gremodel.w2)
			deg_base_energy = mrf.basic_energy_calc(deg_grem_prot, deg_gremodel.w1, deg_gremodel.w2)

			#Now we load this info into an ExChrome type object.
			cur_pop = types.ExChrome[]
			for i in 1:length(success_chrom)
				old = success_chrom[i]
				if rand_weights
					new_chrom = types.ExChrome(old.path, old.full_sequence, old.deg_nuc, indie_deg_maps[i], old.deg_trns, old.deg_trne,
															old.deg_d, old.deg_skip, old.deg_insert, old.mark_nuc, indie_mark_maps[i], old.mark_trns, old.mark_trne,
															old.mark_d, old.mark_skip, old.mark_insert, fdp[i], deg_base_energy, deg_prot_mat[i:i, 1:end], utils.aa_vec_to_seq(deg_prot_mat[i:i, 1:end]), deg_ull[i], deg_pv_w1[i], deg_pv_w2[i:i, 1:end], fmp[i], mark_base_energy,
															mark_prot_mat[i:i, 1:end], utils.aa_vec_to_seq(mark_prot_mat[i:i, 1:end]), mark_ull[i], mark_pv_w1[i], mark_pv_w2[i:i, 1:end], rand())
				else
					new_chrom = types.ExChrome(old.path, old.full_sequence, old.deg_nuc, indie_deg_maps[i], old.deg_trns, old.deg_trne,
															old.deg_d, old.deg_skip, old.deg_insert, old.mark_nuc, indie_mark_maps[i], old.mark_trns, old.mark_trne,
															old.mark_d, old.mark_skip, old.mark_insert, fdp[i], deg_base_energy, deg_prot_mat[i:i, 1:end], utils.aa_vec_to_seq(deg_prot_mat[i:i, 1:end]), deg_ull[i], deg_pv_w1[i], deg_pv_w2[i:i, 1:end], fmp[i], mark_base_energy,
															mark_prot_mat[i:i, 1:end], utils.aa_vec_to_seq(mark_prot_mat[i:i, 1:end]), mark_ull[i], mark_pv_w1[i], mark_pv_w2[i:i, 1:end], 0.5)
				end
				push!(cur_pop, new_chrom)
			end

			#There are occasionally differences between HMM/MRF models of the proteins, mostly because HMMs are more flexible with insertions.
			#This generates mappings between positions in either data type.
			mark_hmm_to_grem_map = std_setup.get_explicit_mapping(mark_grem_prot, mark_hmm)
			deg_hmm_to_grem_map = std_setup.get_explicit_mapping(deg_grem_prot, deg_hmm)

			@debug("Hey the explicit mapping for mark_hmm_to_grem is $mark_hmm_to_grem_map")
			@debug("Hey the explicit mapping for deg_hmm_to_grem is $deg_hmm_to_grem_map")
			#So now everything is in an extended chromosome object.

			changed_seq = length(cur_pop)
			no_change_iters = 0
			iter = 0;

			#Let's keep track of optimization history.
			history_matrix = zeros(Float32, round(Int64, pop_size * 1.2), max_iter)

			println("Beginning long-range optimization.")
			@debug("About to start optimization.")
			@debug(Libc.strftime(time()))

			stop_early = false #a condition to evaluate every (few) iteration(s) to see if it's worth continuing a run. When true, we clean up and move to next gene pair.
			not_worth_continuing = false

			#Okay let's set up the energy normal distributions...

			deg_energy_mu, deg_energy_sig = lookup.get_energy_params(deg_name, "energies/")
			mark_energy_mu, mark_energy_sig = lookup.get_energy_params(mark_name, "energies/")

			mark_energy_normal = Normal(mark_energy_mu, mark_energy_sig)
			deg_energy_normal = Normal(deg_energy_mu, deg_energy_sig)

			while (iter == 0 || (iter < max_iter)) #&& !(stop_early) # && minimum((fitness_values)) > 400 && no_change_iters < 100)
				sfv = sort(fitness_values)
				#best_50 = sfv[50]

				if iter % 50 == 0
					println("Step $(iter) of $(max_iter)...")
				end
				@debug("Iteration $iter out of $max_iter")
				@debug("Mean fitness is: $(mean((fitness_values)))")
				@debug("Bottom five are $(sfv[1:5])")
				@debug("... and top five are $(sfv[end-5:end]).")
				@debug("Number of sequences that were changed: $changed_seq")
				flush(log_io)

				#iterated conditional modes, multi-site keep-track = icm_multi_kt. keep track = store values of changes to sequence.
				if do_cull
					if iter < div(max_iter, 4)
						#println("$iter : no cutoff")
						cur_pop, changed_seq = optimize.icm_multi_kt(cur_pop, deg_gremodel, mark_gremodel) #, deg_energy_normal, mark_energy_normal)
					elseif iter >= div(max_iter, 4) && iter < div(max_iter, 2)
						#println("$iter : cutoff is $(sfv[div(length(sfv), 2)])")
						cur_pop, changed_seq = optimize.icm_multi_kt(cur_pop, deg_gremodel, mark_gremodel; score_cutoff = sfv[div(length(sfv), 2)])
					elseif iter >= div(max_iter, 2)
						#println("$iter : cutoff is $(sfv[div(length(sfv), 4)])")
						cur_pop, changed_seq = optimize.icm_multi_kt(cur_pop, deg_gremodel, mark_gremodel; score_cutoff = sfv[div(length(sfv), 4)])
					end
				else
					cur_pop, changed_seq = optimize.icm_multi_kt(cur_pop, deg_gremodel, mark_gremodel)
				end

				fitness_values = optimize.assess_pop(cur_pop)
				summed_fitness = sum(fitness_values) #one of the things we should be watching.

				push!(last_few_sf, summed_fitness)
				push!(last_few_cs, changed_seq)

				indiv_count = 1
				for indiv in fitness_values #history matrix stores values of function being optimized over time.
					history_matrix[indiv_count, iter + 1] = indiv
					indiv_count += 1
				end

				#This was originally just code to test the pseudolikelihoods decreasing correctly. It does.
				#Now it's still sorta somewhat useful to log a few values of sequences in both deg/mark as optimization goes.
				if (mod(iter, 50) == 0 && iter != 0)
					@debug("Second opinion, these should be ~0.")
					for i_count in 1:5
						#for the record this is not a good way to generate random numbers.
						#doesn't go from 1 to end.
						i = i_count #round(Int64, rand() * (length(cur_pop) - 1) + 1)

						bel_deg = cur_pop[i].deg_prob
						bel_mrk = cur_pop[i].mark_prob
						#Sort of assuming that end-1 is because of a "*" at stop codon being added.
						if length(strip(cur_pop[i].deg_seq, '*')) == deg_nNodes && length(strip(cur_pop[i].mark_seq, '*')) == mark_nNodes
							cal_deg = mrf.psl(strip(cur_pop[i].deg_seq, '*'), deg_gremodel.w1, deg_gremodel.w2, false)
							cal_mrk = mrf.psl(strip(cur_pop[i].mark_seq, '*'), mark_gremodel.w1, mark_gremodel.w2, false)
							@debug("DEG: $bel_deg vs. $cal_deg")
							@debug("MRK: $bel_mrk vs. $cal_mrk")
							@debug(abs(bel_deg - cal_deg) < 1e-2)
							@debug(abs(bel_mrk - cal_mrk) < 1e-2)
						else
							@debug("Proteins not correct length!")
							@debug(cur_pop[i].mark_seq)
							@debug(cur_pop[i].deg_seq)
						end
					end
				end
				if changed_seq == 0
					no_change_iters += 1
				else
					no_change_iters = 0
				end
				iter += 1
				@debug("\n")
			end

			#So we're done the iterated conditional modes step. Now we report everything...
			@debug("Done while loop...")
			@debug(Libc.strftime(time()))

			@debug("Let us save the history matrix!")
			save("$out_path/$(mark_name)_$(deg_name)_$frame/opt_his_mat_$(rand_barcode).jld", "hist_mat", history_matrix)

			deg_significance = false
			mark_significance = false #this sees if we get any hits worth keeping at all. If yes we save full population.
			all_cal_deg_scores = Float64[]
			all_cal_mark_scores = Float64[]
			#If no we save just a few (top 10) individuals. These are top-6 overall + top-2 in deg and top-2 in mark.
			out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/all_final_fitness_$(rand_barcode).txt", "w")
			write(out_file, "Ind. #\tMark Score\tDeg Score\tMark Sign\tDeg Sign\n")
			for last_indi in 1:length(cur_pop)
				rep_deg = cur_pop[last_indi].deg_prob
				rep_mrk = cur_pop[last_indi].mark_prob
				if length(strip(cur_pop[last_indi].deg_seq, '*')) == deg_nNodes && length(strip(cur_pop[last_indi].mark_seq, '*')) == mark_nNodes
					cal_deg = abs(mrf.psl(strip(cur_pop[last_indi].deg_seq, '*'), deg_gremodel.w1, deg_gremodel.w2, false))
					cal_mrk = abs(mrf.psl(strip(cur_pop[last_indi].mark_seq, '*'), mark_gremodel.w1, mark_gremodel.w2, false))
					special_deg = "-"
					if cal_deg < mu_deg + sig_deg
						special_deg = "***!"
						deg_significance = true
					elseif cal_deg < mu_deg + 2*sig_deg
						special_deg = "**."
						deg_significance = true
					elseif cal_deg < mu_deg + 3*sig_deg
						special_deg = "*?"
						deg_significance = true
					elseif cal_deg < mu_deg + 4*sig_deg #we don't mind keeping these.
						deg_significance = true
					end
					special_mark = "-"
					if cal_mrk < mu_mark + sig_mark
						special_mark = "***!"
						mark_significance = true
					elseif cal_mrk < mu_mark + 2*sig_mark
						special_mark = "**."
						mark_significance = true
					elseif cal_mrk < mu_mark + 3*sig_mark
						special_mark = "*?"
						mark_significance = true
					elseif cal_mrk < mu_mark + 4*sig_mark #we don't mind keeping these.
						mark_significance = true
					end
					push!(all_cal_mark_scores, cal_mrk)
					push!(all_cal_deg_scores, cal_deg)
					write(out_file, "$(last_indi)\t$(cal_mrk)\t$(cal_deg)\t$special_mark\t$special_deg\n")
				else
					push!(all_cal_mark_scores, Inf)
					push!(all_cal_deg_scores, Inf)
					write(out_file, "$(last_indi)\t-1.0\t-1.0\t/\t/\n")
				end
			end
			close(out_file)

			cal_deg_p = sortperm(all_cal_deg_scores)
			cal_mark_p = sortperm(all_cal_mark_scores)
			cal_both_p = sortperm(fitness_values)

			@debug(Libc.strftime(time()))
			if true #deg_significance && mark_significance
				@debug("Let's save our optimized population...")
				saved_pop = convert_to_saveable(cur_pop)
				save("$out_path/$(mark_name)_$(deg_name)_$frame/saved_pop_$(rand_barcode).jld", "variants", saved_pop)
			else #We want just a subset.
				@debug("We'll save 12 interesting members of the population.")
				selected_pop = ExChrome[]
				for ijk in 1:3
					push!(selected_pop, cur_pop[cal_deg_p[ijk]])
					push!(selected_pop, cur_pop[cal_mark_p[ijk]])
				end
				for ijk in 1:6
					push!(selected_pop, cur_pop[cal_both_p[ijk]])
				end
				save("$out_path/$(mark_name)_$(deg_name)_$frame/top_pop_$(rand_barcode).jld", "variants", selected_pop)
			end

			@debug("And we'll also output our top twelve sequences...")
			out_file = open("$out_path/$(mark_name)_$(deg_name)_$frame/top_twelve_$(rand_barcode).fa", "w")

			for ijk in 1:3
				deg_ind = cal_deg_p[ijk]
				write(out_file, ">TOP_DEG_d (ind $deg_ind ) (fit: $(fitness_values[deg_ind])) (deg: $(all_cal_deg_scores[deg_ind])) (mark: $(all_cal_mark_scores[deg_ind]))\n")
				write(out_file, "$(cur_pop[deg_ind].deg_seq)\n")
				write(out_file, ">TOP_DEG_m (ind $deg_ind ) (fit: $(fitness_values[deg_ind])) (deg: $(all_cal_deg_scores[deg_ind])) (mark: $(all_cal_mark_scores[deg_ind]))\n")
				write(out_file, "$(cur_pop[deg_ind].mark_seq)\n")
			end
			for ijk in 1:3
				deg_ind = cal_mark_p[ijk]
				write(out_file, ">TOP_MARK_d (ind $deg_ind ) (fit: $(fitness_values[deg_ind])) (deg: $(all_cal_deg_scores[deg_ind])) (mark: $(all_cal_mark_scores[deg_ind]))\n")
				write(out_file, "$(cur_pop[deg_ind].deg_seq)\n")
				write(out_file, ">TOP_MARK_m (ind $deg_ind ) (fit: $(fitness_values[deg_ind])) (deg: $(all_cal_deg_scores[deg_ind])) (mark: $(all_cal_mark_scores[deg_ind]))\n")
				write(out_file, "$(cur_pop[deg_ind].mark_seq)\n")
			end
			for ijk in 1:6
				deg_ind = cal_both_p[ijk]
				write(out_file, ">TOP_GEN_d (ind $deg_ind ) (fit: $(fitness_values[deg_ind])) (deg: $(all_cal_deg_scores[deg_ind])) (mark: $(all_cal_mark_scores[deg_ind]))\n")
				write(out_file, "$(cur_pop[deg_ind].deg_seq)\n")
				write(out_file, ">TOP_GEN_m (ind $deg_ind ) (fit: $(fitness_values[deg_ind])) (deg: $(all_cal_deg_scores[deg_ind])) (mark: $(all_cal_mark_scores[deg_ind]))\n")
				write(out_file, "$(cur_pop[deg_ind].mark_seq)\n")
			end
			close(out_file)
		end

		@debug("REALLY DONE!")

		@debug(Libc.strftime(time()))

	end
end

function parse_commandline()
	s = ArgParseSettings()
	@add_arg_table s begin
		"commands"
			help = "positional argument, path to file with tab-delimited commands to run (default: example.txt)"
			default = "example.txt"
		"--num"
			help = "optional, experimentally used for running multiple jobs at once"
			default = "0"
		"--threads"
			help = "optional, experimentally used for running multiple jobs at once"
			default = "1"
	end
	return parse_args(s)
end

function run_file()
	parsed_args = parse_commandline()

	command_file = parsed_args["commands"]
	num = parsed_args["num"]
	threads = parsed_args["threads"]

	run_i = parse(Int64, num)
	NUM_THREADS = parse(Int64, threads)

	if NUM_THREADS > 1
		sleep(rand() * 2) #randomly delay start to de-synchronize starts of runs.
	end

	if isfile(command_file)
		in_file = open(command_file)
		in_read = readlines(in_file)
		close(in_file)

		#This is an additional file that records errors independent of individual log files.
		problem_file = open("problem_runs_$(run_i).txt", "w")

		line_count = 0
		for line in in_read #
			line_count += 1
			if line[1] != '#' && (line_count % NUM_THREADS == run_i)
				try
					run_args = split(line, "\t")
					out_dir, short, long, short_jld, long_jld, short_hmm, long_hmm, pop_size, frame, num_iter = run_args

					rand_barcode = Random.randstring()

					the_out_path = "$out_dir/$(short)_$(long)_$frame/"
					if !(isdir(the_out_path))
						run(`mkdir -p $the_out_path`)
					end

					log_io = open("$out_dir/$(short)_$(long)_$frame/log_$(rand_barcode).txt", "w+")

					logger = SimpleLogger(log_io, Logging.Debug)
					global_logger(logger)

					with_logger(logger) do
						pop_size = parse(Int64, pop_size)
						num_iter = parse(Int64, num_iter)
						set_up_and_optimize(log_io, rand_barcode, out_dir, short, long, short_jld, long_jld, short_hmm, long_hmm, pop_size, frame, num_iter)
					end

					flush(log_io)
					close(log_io)
				catch y
					write(problem_file, "$line $y")
				end

			end
		end
		close(problem_file)
	else
		println("$command_file does not exist. Exiting.")
	end

end

@time run_file()

end

max_stars_repo_name: khyox/CAMEOS

max_issues_repo_name: BiosecSFA/CAMEOS

max_forks_repo_name: BiosecSFA/CAMEOS


text:
using Distributions
using LinearAlgebra
using PointwiseKDEs
using Statistics
using Test

@testset "PointwiseKDE Tests" begin
    @testset "3D Gaussian" begin
        mu = randn(3)
        Sigma = randn(3,3)
        Sigma = Sigma' * Sigma

        CSigma = cholesky(Sigma)

        pts = mu .+ CSigma.L*randn(3, 1000)

        kde = PointwiseKDE(pts)
        draw_pts = rand(kde, 1000)

        mu_draw = mean(draw_pts, dims=2)
        sigma_draw = cov(draw_pts, dims=2)

        @testset "Mean Agrees" begin
            @test all(abs.(CSigma.L \ (mu .- mu_draw)) .< 0.1)
        end

        @testset "Average density is close" begin
            dist = MvNormal(mu, Sigma)
            logp_kde = [logpdf(kde, draw_pts[:,i]) for i in 1:size(draw_pts,2)]
            logp_mvn = [logpdf(dist, draw_pts[:,i]) for i in 1:size(draw_pts,2)]
            @test abs(mean(logp_kde .- logp_mvn)) < 0.1
        end
    end
end


text:
export PumiMesh3

mutable struct PumiMesh3{T1} <: PumiMesh3CG{T1}   # 3d pumi mesh, tetrahedron only
  m_ptr::Ptr{Void}  # pointer to mesh
  mshape_ptr::Ptr{Void} # pointer to mesh's FieldShape
  f_ptr::Ptr{Void} # pointer to apf::field for storing solution during mesh adaptation

  shape_type::Int  # type of shape functions

  vert_Nptr::Ptr{Void}  # numbering of vertices
  edge_Nptr::Ptr{Void}  # numbering of edges
  face_Nptr::Ptr{Void} # numbering of faces
  el_Nptr::Ptr{Void}  # numbering of elements (faces)

  numVert::Int  # number of vertices in the mesh
  numEdge::Int # number of edges in the mesh
  numFace::Int # number of faces in the mesh
  numEl::Int  # number of elements (regions)
  order::Int # order of shape functions
  numDof::Int # number of degrees of freedom
  numNodes::Int  # number of nodes
  numDofPerNode::Int  # number of dofs per node
  numBoundaryFaces::Int # number of faces on the exterior boundary
  numInterfaces::Int # number of internal interfaces
  numNodesPerElement::Int  # number of nodes per element
  numNodesPerType::Array{Int, 1}  # number of nodes classified on each vertex, edge, face, length 4
  numEntities::Array{Int, 1}  # number of entities of each type in the mesh, length 4
  numEntitiesPerElement::Array{Int, 1}  # number of vertices, edges, and regions in each element

  # parallel info
  comm::MPI.Comm  # MPI Communicator
  myrank::Int  # MPI rank, zero based
  commsize::Int # MPI comm size

  # hold pointers to mesh entities
  verts::Array{Ptr{Void},1}  # holds pointers to mesh vertices
  edges::Array{Ptr{Void},1}  # pointers to mesh edges
  faces::Array{Ptr{Void},1}  # pointers to mesh faces
  elements::Array{Ptr{Void},1}  # pointers to regions

  dofnums_Nptr::Ptr{Void}  # pointer to Numbering of dofs (result of reordering)
  boundary_nums::Array{Int, 2}  # array of [element number, edgenumber] for each edge on the boundary

  bndryfaces::Array{Boundary, 1}  # store data on external boundary of mesh
  interfaces::Array{Interface, 1}  # store data on internal edges

  coords::Array{T1, 3}  # store coordinates of all nodes
  dxidx::Array{T1, 4}  # store scaled mapping jacobian
  jac::Array{T1,2}  # store mapping jacobian output

  dofs::Array{Int, 3}  # store dof numbers of solution array to speed assembly

 function PumiMesh3{T1}(dmg_name::AbstractString, smb_name::AbstractString, order, sbp::AbstractSBP; dofpernode=1, shape_type=0)  where {T1}
  # construct pumi mesh by loading the files named
  # dmg_name = name of .dmg (geometry) file to load (use .null to load no file)
  # smb_name = name of .smb (mesh) file to load
  # order = order of shape functions
  # dofpernode = number of dof per node, default = 1
  # shape_type = type of shape functions to use, 0 = lagrange, 1 = SBP

  mesh = new{T1}()
  mesh.comm = MPI.COMM_WORLD
  mesh.commsize = MPI.Comm_size(MPI.COMM_WORLD)
  mesh.myrank = MPI.Comm_rank(MPI.COMM_WORLD)

  @assert mesh.commsize == 1  # not paralllezed yet

  mesh.numDofPerNode = dofpernode
  mesh.order = order
  mesh.numNodesPerElement = getNumNodes(order)
  mesh.shape_type = shape_type
  @time mesh.numEntities, mesh.m_ptr, mesh.mshape_ptr, n_arr = apf.init(dmg_name, smb_name, order, shape_type=shape_type)
  mesh.f_ptr = apf.findField(mesh.m_ptr, "solution_field")
  if mesh.f_ptr == C_NULL
    @time mesh.f_ptr = apf.createPackedField(mesh, "solution_field", dofpernode)
  end

  mesh.numVert = convert(Int, mesh.numEntities[1])
  mesh.numEdge = convert(Int,  mesh.numEntities[2])
  mesh.numFace = convert(Int, mesh.numEntities[3])
  mesh.numEl = convert(Int, mesh.numEntities[4])

  # get number of nodes on each type
  mesh.numNodesPerType = zeros(Int, 4)
  mesh.numNodesPerType[1] = apf.countNodesOn(mesh.mshape_ptr, 0)  # number of nodes on a vertex
  mesh.numNodesPerType[2]= apf.countNodesOn(mesh.mshape_ptr, 1) # on edge
  mesh.numNodesPerType[3] = apf.countNodesOn(mesh.mshape_ptr, 2) # on face
  mesh.numNodesPerType[4]= apf.countNodesOn(mesh.mshape_ptr, 4) # on element interior
  mesh.numEntitiesPerElement = [4, 6, 4, 1]

  mesh.numNodesPerElement = 0
  for i=1:4
    mesh.numNodesPerElement += mesh.numNodesPerType[i]*mesh.numEntitiesPerElement[i]
  end

  # get pointers to mesh entity numberings
  mesh.vert_Nptr = n_arr[1] #getVertNumbering()
  mesh.edge_Nptr = n_arr[2] #getEdgeNumbering()
  mesh.face_Nptr = n_arr[3] #getFaceNumbering()
  mesh.el_Nptr = n_arr[4] #getElNumbering()

   mesh.verts = Array{Ptr{Void}}(mesh.numVert)
  mesh.edges = Array{Ptr{Void}}(mesh.numEdge)
  mesh.faces = Array{Ptr{Void}}(mesh.numFace)
  mesh.elements = Array{Ptr{Void}}(mesh.numEl)
  mesh.dofnums_Nptr = apf.createNumberingJ(mesh, "reordered dof numbers",
                                           dofpernode, mesh.mshape_ptr)  # 1 dof per node

  # get pointers to all MeshEntities
  # also initilize the field to zero
  for (i, e) in enumerate(apf.MeshIterator(mesh.m_ptr, 0))
    mesh.verts[i] = e
  end

  for (i, e) in enumerate(apf.MeshIterator(mesh.m_ptr, 1))
    mesh.edges[i] = e
  end

  for (i, e) in eumerate(apf.MeshIterator(mesh.m_ptr, 2))
    mesh.faces[i] = e
  end

  for (i, e) in enumerate(apf.MeshIterator(mesh.m_ptr, 3))
    mesh.elements[i] = e
  end

  # use partially constructed mesh object to populate arrays
  checkConnectivity(mesh)

#  numberDofs(mesh)

  countBoundaryFaces(mesh)
  mesh.bndryfaces = Array{Boundary}(mesh.numBoundaryFaces)
  getBoundaryArray(mesh)

  numberDofs(mesh)
  getDofNumbers(mesh)

  countBoundaryFaces(mesh)
  mesh.bndryfaces = Array{Boundary}(mesh.numBoundaryFaces)
  getBoundaryArray(mesh)

  mesh.numInterfaces = mesh.numFace - mesh.numBoundaryFaces
  mesh.interfaces = Array{Interface}(mesh.numInterfaces)
  getInterfaceArray(mesh)

  getCoordinatesAndMetrics(mesh, sbp)
#

#=
  mesh.bndryfaces = Array{Boundary}(mesh.numBoundaryFaces)
  getBoundaryArray(mesh)

  mesh.numInterfaces = mesh.numEdge - mesh.numBoundaryFaces
  mesh.interfaces = Array{Interface}(mesh.numInterfaces)
  getInterfaceArray(mesh)

  getCoordinates(mesh, sbp)  # store coordinates of all nodes into array
  getDofNumbers(mesh)  # store dof numbers

  mesh.dxidx = Array{T1}(2, 2, sbp.numnodes, mesh.numEl)
  mesh.jac = Array{T1}(sbp.numnodes, mesh.numEl)
  mappingjacobian!(sbp, mesh.coords, mesh.dxidx, mesh.jac)
=#

  printStats(mesh)

  apf.writeVtkFiles("mesh_complete", mesh.m_ptr)
  return mesh
  end

end

function getNumNodes(order::Integer)
# get the number of nodes on an element
# assumes SBP elements

  nnodes = 0
  if order == 1
    nnodes = 4
  elseif order == 2
    nnodes = 11
  elseif order == 3
    nnodes = 24
  elseif order == 4
    nnodes = 45
  else
    println("Warning, unsupported order elements are requested")
  end

  return nnodes
end

function countBoundaryFaces(mesh::PumiMesh3)
  # count boundary face
  # store array of [element number, global edge number]
#  resetFaceIt()
  bnd_faces_cnt = 0
  bnd_faces = Array{Int}(mesh.numFace, 2)
  it = apf.MeshIterator(mesh.m_ptr, 2)
  for i=1:mesh.numFace
    face_i = apf.iterate(mesh.m_ptr, it)
    numRegion = apf.countAdjacent(mesh.m_ptr, face_i, 3)  # should be count upward

    if numRegion == 1  # if an exterior edge
      elements = apf.getAdjacent(numRegion)
      elnum = apf.getNumberJ(mesh.el_Nptr, elements[1], 0, 0) + 1

#      println("face number ", i, " is part of only one element")

      bnd_faces_cnt += 1
      bnd_faces[bnd_faces_cnt, 1] = elnum
      bnd_faces[bnd_faces_cnt, 2] = i
    end
  end
  apf.free(mesh.m_ptr, it)

  mesh.boundary_nums = bnd_faces[1:bnd_faces_cnt, :] # copy, bad but unavoidable
  mesh.numBoundaryFaces = bnd_faces_cnt

return nothing

end  # end function

function getBoundaryArray(mesh::PumiMesh3)
# get an array of type Boundary for SBP
# creating an an array of a user defined type seems like a waste of memory operations
# bnd_array is a vector of type Boundary with length equal to the number of edges on the boundary of the mesh

#  mesh.bndryfaces = Array{Boundary}(mesh.numBoundaryFaces)

  for i=1:mesh.numBoundaryFaces
    elnum = mesh.boundary_nums[i,1]
    facenum_global = mesh.boundary_nums[i,2]
    facenum_local = getBoundaryEntityLocalNum(mesh, facenum_global, 2)
    mesh.bndryfaces[i] = Boundary(elnum, facenum_local)
  end

  return nothing
end

function getBoundaryEntityLocalNum(mesh::PumiMesh3, entity_num::Integer, entity_dim::Integer)
# gets the local edge number of a specified edge that is on the boundary
# of the mesh
# entity_num is an the  num from the output of getBoundaryEdgeNums() (ie. the global entity number)
# entity_dim is the dimension of the entity (0 = vertex, 1 = edge, 2 = face)
# the local edge number is the edge number within the element (1st,s 2nd, 3rd...)
  if entity_dim == 1
    entity_i = mesh.edges[entity_num]
  elseif entity_dim == 2
    entity_i = mesh.faces[entity_num]
  else
    println("unsupported entity dimension requested for getBoundaryEntityLocalNum")
    println("requested entity dimension = ", entity_dim)
  end

  # get mesh face or region associated with edge
  apf.countAdjacent(mesh.m_ptr, entity_i, entity_dim+1)
  parent = apf.getAdjacent(1)[1]  # get the single face (not an array)

  local_num = apf.getEntityLocalNumber(mesh.m_ptr, entity_i, parent, entity_dim, entity_dim+1)

  return local_num

end

function numberDofs(mesh::PumiMesh3)
# assign dof numbers to entire mesh
# calculates numNodes, numDof
# assumes mesh elements have already been reordered
# this method minimizes the assembly/disassembly time
# because blocks of entries in res go into res_vec

  # calculate total number of nodes
  # Continuous Galerkin only
  numnodes = 0
  for i=1:4
    numnodes += mesh.numNodesPerType[i]*mesh.numEntities[i]
  end

  numDof = mesh.numDofPerNode*numnodes

  if (numDof > 2^30)
    println("Warning: too many dofs, renumbering will fail")
  end

  # save numbers to mesh
  mesh.numNodes = numnodes
  mesh.numDof = numDof

  # initally number all dofs as numDof+1 to 2*numDof
  # this allows quick check to see if somthing is labelled or not

#  resetAllIts2(mesh.m_ptr)
  # mesh iterator increment, retreval functions
#  iterators_inc = [incrementVertIt, incrementEdgeIt, incrementFaceIt, incrementElIt]
#  iterators_get = [getVert, getEdge, getFace, getEl]
  num_entities = mesh.numEntities
  num_nodes_entity = mesh.numNodesPerType

  #  curr_dof = 1
  curr_dof = numDof+1
  for etype = 1:length(num_nodes_entity) # loop over entity types
    if (num_nodes_entity[etype] != 0)  # if no nodes on this type of entity, skip
      it = apf.MeshIterator(mesh.m_ptr, etype - 1)
      for entity = 1:num_entities[etype]  # loop over all entities of this type
#	println("  entity number: ", entity)
        entity_ptr = apf.iterate(mesh.m_ptr, it)
#	entity_ptr = iterators_get[etype]()  # get entity

	for node = 1:num_nodes_entity[etype]
#	  println("    node : ", node)
	  for dof = 1:mesh.numDofPerNode
	    apf.numberJ(mesh.dofnums_Nptr, entity_ptr, node-1, dof-1, curr_dof)
#	    println("      entity ", entity_ptr, " labelled ", curr_dof)
	    curr_dof += 1
	  end  # end loop over dof
	end  # end loop over node

#      iterators_inc[etype]()
      end  # end loop over entitiesaa
      apf.free(mesh.m_ptr, it)
    end  # end if
  end  # end loop over entity types

  # final dof numbering
# move all if statements out one for loop (check only first dof on each node)
  curr_dof = 1
  for i=1:mesh.numEl
#    println("element number: ", i)

    el_i = mesh.elements[i]
    # get vertices, edges for this element
    (verts_i, numVert) = apf.getDownward(mesh.m_ptr, mesh.elements[i], 0)
#    println("verts_i = ", verts_i)
#    println("mesh.verts = ", mesh.verts)
    (edges_i, numEdge) = apf.getDownward(mesh.m_ptr, mesh.elements[i], 1)

    if mesh.numEntitiesPerElement[3] > 1  # don't call get downward if el_i is a face
      (face_i, numFace) = apf.getDownward(mesh.m_ptr, mesh.elements[i], 2)
      regions_i = [mesh.elements[i]]  # make this an array for consistency
      numRegions = 1
    else  # this is a 2d mesh
      face_i = [mesh.elements[i]]
      numFace = 1
      numRegion = 0
    end

    # loop over vertices
    for j=1:numVert  # loop over vertices, edges
#      println("  vertex number: ", j)
      vert_j = verts_i[j]
#      println("  vert_j = ", vert_j)
#      println("  edge_j = ", edge_j)
      for k=1:mesh.numDofPerNode # loop over vertex dofs
#	println("    dof number: ", k)
#	println("    vert_j = ", vert_j)
        dofnum_k = apf.getNumberJ(mesh.dofnums_Nptr, vert_j, 0, k-1)
	if dofnum_k > numDof  # still has initial number
	  # give it new (final) number
	  apf.numberJ(mesh.dofnums_Nptr, vert_j, 0, k-1, curr_dof)
#	  println("numbering a vertex ", curr_dof)
	  curr_dof += 1
	end
      end  # end loop over vertex dofs
    end  # end loop over vertices

    # now loop over edges
    for j=1:numEdge
#      println("edge number: ", j)
      edge_j = edges_i[j]
      # loop over nodes on edge
      for k=1:num_nodes_entity[2]  # loop over all nodes on the edge
	for p=1:mesh.numDofPerNode  # loop over dofs
	  dofnum_p = apf.getNumberJ(mesh.dofnums_Nptr, edge_j, k-1, p-1)
	  if dofnum_p > numDof  # still has initial number
	    # give it new (final) number
	    apf.numberJ(mesh.dofnums_Nptr, edge_j, k-1, p-1, curr_dof)
	    curr_dof += 1
	  end
	end
      end  # en loop over nodes on edge
    end  # end loop over edges

    # now loop over faces
    for j=1:numFace
#      println("face number: ", j)
      face_j = face_i[j]
      for k=1:num_nodes_entity[3]  # loop over nodes on faces
	for p=1:mesh.numDofPerNode  # loop over dofs
	  dofnum_p = apf.getNumberJ(mesh.dofnums_Nptr, face_j, k-1, p-1)
	  if dofnum_p > numDof
	    apf.numberJ(mesh.dofnums_Nptr, face_j, k-1, p-1, curr_dof)
	    curr_dof += 1
	  end
	end
      end  # end loop over face nodes
    end  # end loop over faces

    # now loop over region
    for j=1:numRegions
#      println("region number: ", j)
      region_j = regions_i[j]
      for k=1:num_nodes_entity[4]
	for p=1:mesh.numDofPerNode
	  dofnum_p = apf.getNumberJ(mesh.dofnums_Nptr, region_j, k-1, p-1)
	  if dofnum_p > numDof
	    apf.numberJ(mesh.dofnums_Nptr, region_j, k-1, p-1, curr_dof)
	    curr_dof += 1
	  end
	end
      end  # end loop over nodes on the region
    end  # end loop over region

  end  # end loop over elements

  if (curr_dof -1) != numDof
    println("Warning: number of dofs assigned is not equal to teh expected number")
    println("number of dofs assigned = ", curr_dof-1, " ,expected number = ", numDof)
  end

#  resetAllIts2(mesh.m_ptr)
  return nothing

end

function getDofNumbers(mesh::PumiMesh3)
# populate array of dof numbers, in same shape as solution array u (or q)

mesh.dofs = Array{Int}(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)

for i=1:mesh.numEl
  dofnums = getGlobalNodeNumbers(mesh, i)

  for j=1:mesh.numNodesPerElement
    for k=1:mesh.numDofPerNode  # loop over dofs on the node
      mesh.dofs[k, j, i] = dofnums[k,j]
    end
  end
end

#println("mesh.dof = ", mesh.dofs)

return nothing

end

function getGlobalNodeNumbers(mesh::PumiMesh3, elnum::Integer)
# gets global node numbers of all dof on all nodes of the element
# output formap is array [numdofpernode, nnodes]  (each column contains dof numbers for a node)
el_i = mesh.elements[elnum]
#type_i = apf.getType(mesh.m_ptr, el_i)

# use apf::getElementNumbers to ensure everything is in the proper orientation

  nnodes = mesh.numNodesPerElement
  numdof = nnodes*mesh.numDofPerNode
  nums = zeros(Int32, numdof)

apf.getElementNumbers(mesh.dofnums_Nptr, el_i, numdof, nums::Array{Int32, 1})

# reshape into output format
# underlying memory is shared
dofnums_reshaped = reshape(nums, mesh.numDofPerNode, nnodes)

return dofnums_reshaped

end

function getCoordinates(mesh::PumiMesh3, sbp::AbstractSBP)
# populate the coords array of the mesh object

mesh.coords = Array{Float64}(3, sbp.numnodes, mesh.numEl)

println("entered getCoordinates")

coords_i = zeros(3, 4)
coords_it = zeros(4, 3)
for i=1:mesh.numEl  # loop over elements

  el_i = mesh.elements[i]
  (sizex, sizey) = size(coords_i)
  apf.getElCoords(mesh.m_ptr, el_i, coords_i, sizex, sizey)  # populate coords

#  println("coords_i = ", coords_i)

  coords_it[:,:] = coords_i.'
 # println("coords_it = ", coords_it)
  mesh.coords[:, :, i] = calcnodes(sbp, coords_it)
#  println("mesh.coords[:,:,i] = ", mesh.coords[:,:,i])
end

return nothing

end

function getInterfaceArray(mesh::PumiMesh3)
# get array of [elementL, elementR, faceL, faceR] for each internal face,
# where elementL and R are the elements that use the face, faceL R are the
# local number of the edge within the element
# interfaces is the array to be populated with the data, it must be
# number of internal edges by 1

  # only need internal boundaries (not external)
#  num_ext_edges = size(getBoundaryEdgeNums(mesh))[1]  # bad memory efficiency
#  num_int_edges = getNumEdges(mesh) - num_ext_edges

#  new_bndry = Boundary(2, 3)
#   println("new_bndry = ", new_bndry)

#  new_interface = Interface(1, 2, 3, 4)
#   println("new_interface = ", new_interface)

#   println("num_int_edges = ", num_int_edges)

#  interfaces = Array{typeof(new_interface}(num_int_edges)
  # get number of nodes affecting an edge
  num_face_nodes = apf.countAllNodes(mesh.mshape_ptr, 2)

  # store whether the node map has been calculated for each case
  # index = relative rotation of faces + 1
  nodemap_mask = zeros(Bool, 3)
  nodemaps = Array{Array{Uint8,1}}(3)

#  nodemap = collect(num_edge_nodes:(-1):1)

  pos = 1 # current position in interfaces
  for i=1:mesh.numFace
#     println("i = ", i)
#     println("pos = ", pos)
    # get number of elements using the edge
    adjacent_nums, num_adjacent = getAdjacentEntityNums(mesh, i, 2, 3)
#     println("num_adjacent = ", num_adjacent)
#     println("adjacent_nums = ", adjacent_nums)
    if num_adjacent > 1  # internal edge
#       println("this is an internal edge")
      element1 = adjacent_nums[1]
      element2 = adjacent_nums[2]

#      coords_1 = x[:, :, element1]
#      coords_2 = x[:, :, element2]

      coords_1 = zeros(3,4)
      coords_2 = zeros(3,4)
      apf.getElCoords(mesh.m_ptr, mesh.elements[element1], coords_1, 3, 4)
      apf.getElCoords(mesh.m_ptr, mesh.elements[element2], coords_2, 3, 4)

      # calculate centroid
      centroid1 = sum(coords_1, 2)
      centroid2 = sum(coords_2, 2)

      if abs(centroid1[1] - centroid2[2]) < 1e-10  # if big enough difference
        if centroid1[1] < centroid2[1]
    elementL = element1
    elementR = element2
        else
    elementL = element2
    elementR = element1
        end
      else  # use y coordinate to decide
        if centroid1[2] < centroid2[2]
    elementL = element1
    elementR = element2
        else
    elementL = element2
    elementR = element1
        end
      end

      faceL = apf.getEntityLocalNumber(mesh.m_ptr, mesh.faces[i], mesh.elements[elementL], 2, 3)
      faceR = apf.getEntityLocalNumber(mesh.m_ptr, mesh.faces[i], mesh.elements[elementR], 2, 3)
#      edgeL = getFaceLocalNum(mesh, i, elementL)
#      edgeR = getFaceLocalNum(mesh, i, elementR)

      # construct nodemap if needed
      rel_rotate = calcRelRotate(mesh, elementL, elementR, i)
      if !nodemap_mask[rel_rotate + 1]  # if nodemap not yet constructed
	nodemaps[rel_rotate + 1] = constructNodemap(mesh, rel_rotate)
	nodemap_mask[rel_rotate + 1] = true
      end

      mesh.interfaces[pos] = Interface(elementL, elementR, faceL, faceR, nodemaps[rel_rotate + 1])
#       println("updating pos")
      pos += 1

  #    print("\n")

    end  # end if internal edge

#     print("\n")
  end  # end loop over edges

  return nothing

end  # end function

function calcRelRotate(mesh::PumiMesh3, elementL::Integer, elementR::Integer, facenum::Integer)
# calculate the rotation of faceR relative to faceL, where faceL and faceR
# are shared between the two element
# this uses getAligmnet, which describes how to transform the face from
# its current orientation to the canonical orientation as part of its parent
# element
# the difference in the required rotations is the relative rotation
# elementL : global element index
# elementR : global element index
# faceL : local face number of shared face on elementL
# faceR : local face number of shared face on elementR
# face : global face number of the shared face

eL = mesh.elements[elementL]
eR = mesh.elements[elementR]
face = mesh.faces[facenum]

# get rotations to bring faces into canonical orientation
# for each element
# we can ignore flip because one of the faces will always be flipped
which, flipL, rotateL = apf.getAlignment(mesh.m_ptr, eL, face)
which, flipR, rotateR = apf.getAlignment(mesh.m_ptr, eR, face)

if flipL && flipR
  println("Warning both faces are flipped, nodemap may  be incorrect")
end

# sum rotations because they each rotate ccw in their parent
# element's orientation, so they rotate in opposite directions
rel_rotate = rotateR + rotateL
#println("rel_rotate before wrapping = ", rel_rotate)
rel_rotate = wrapNumber(rel_rotate, 0, 2)
#println("rel_rotate after wrapping = ", rel_rotate)

return rel_rotate

end  # end function

function constructNodemap(mesh::PumiMesh3, rel_rotate::Integer)
# calculate node map from elementL in canonical order to the canonical order
# on elementR
# both orders are canonical because the dofs array ensure everything
# is in canonical order for the element
# element must be a tetrahedron
# make nodemap UInt8?
# calculate total number of nodes on a face
num_nodes_face = 0

numEntitiesPerFace = [3, 3, 1]

for i=1:3
  num_nodes_face += numEntitiesPerFace[i]*mesh.numNodesPerType[i]
end

# get number of nodes classified on edges and face
nnodes_edge = mesh.numNodesPerType[2]
nnodes_face = mesh.numNodesPerType[3]

nodemap = zeros(Uint8, num_nodes_face)

for i=1:3  # loop over verts
  index_i = wrapNumber(i + rel_rotate, 1, 3)  # correct for relative rotation
#  println("after first wrap, index_i = ", index_i)
  index_i = 5 - index_i  # permuate because one triangle is flipped
#  println("pre wrapping, vertex index_i = ", index_i)
  index_i = wrapNumber(index_i, 1, 3)
#  println("post wrapping, vert index_i = ", index_i)
  nodemap[i] = index_i
end

for i=1:3  # loop over edges
  edgeR = 3 - i + 1 - rel_rotate

#  println("before wrap, edgeR = ", edgeR)
  edgeR = wrapNumber(edgeR, 1, 3)
#  println("after wrap, edgeR = ", edgeR)
  for j=1:nnodes_edge
    nodemap_index = 3 + nnodes_edge*(i-1) + j
    nodemap_val = 3 + nnodes_edge*(edgeR-1) + (nnodes_edge -j + 1)
#    println("nodemap_index = ", nodemap_index)
#    println("nodemap_val = ", nodemap_val)
    nodemap[nodemap_index] = nodemap_val
  end
end

#=
# threshold = location of wraparound
threshold = rel_rotate + 1
sum1 = 2 + (mesh.order - 2)*rel_rotate  # index sum below threshold
sum2 = 2 + (mesh.order - 2)*rel_rotate + nnodes_face  # index sum above threshold
=#

if ((rel_rotate % 2) == 0)
  threshold = rel_rotate + 1
  sum1 = nnodes_face - 1
  sum2 = 2*nnodes_face - 1
else  # rel_rotate is odd
  threshold = 0  # always above threshold
  sum1 = 0  # irrelevent
  sum2 = nnodes_face + 1
end

if mesh.order == 3
  threshold = 0
  sum1 = 2  # irrelevent
  sum2 = 2
end

if mesh.order == 4
  if rel_rotate == 0
    threshold = 1
    sum1 = 2
    sum2 = 5
  elseif rel_rotate == 1
    threshold == 0
    sum1 = 4  # irrelevent
    sum2 = 4
  else  # rel_rotate = 2
    threshold = 2
    sum1 = 3
    sum2 = 5
  end
end

nodemap_index = 3 + 3*nnodes_edge
for j=1:nnodes_face  # loop over interior face points
  if j <= threshold  # choose which sum to use
    sum = sum1
  else
    sum = sum2
  end
#  println("face loop j = ", j, " sum = ", sum)
  nodemap[nodemap_index + j] = nodemap_index + sum - j
end

return nodemap

end


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

__precompile__(true)

module LegendDataTypes

using ArraysOfArrays
using ElasticArrays
using EncodedArrays
using RadiationDetectorSignals
using RadiationSpectra
using RecipesBase
using StaticArrays
using StatsBase
using StructArrays
using UnsafeArrays
using Unitful

import Tables
import TypedTables
import Colors

include("daq.jl")
include("type_registry.jl")
include("array_codecs.jl")
include("abstract_io.jl")
include("output_generation.jl")
include("radware_sigcompress.jl")
include("utils.jl")

const array_codecs = TypeRegistry{AbstractArrayCodec}()

function __init__()
    array_codecs[:varlend_diff1] = VarlenDiffArrayCodec
    array_codecs[:radware_sigcompress] = RadwareSigcompress
end

end # module