text:
using Hive
using Test
using DelimitedFiles
using Dates
using Random
using Statistics

function open_database(f::Function)
    session = HiveSession()
    try
        result = execute(session, "use default")
        @test round(Int, result) == 0
        f(session)
    finally
        close(session)
    end
end

nthcol(cc, n) = cc[n][2]
ncells(cc) = nrows(cc) * ncols(cc)
ncols(cc) = length(cc)
nrows(cc) = isempty(cc) ? 0 : length(nthcol(cc,1))
firstcell(cc) = (nrows(cc) > 0) ? nthcol(cc,1)[1] : nothing

function fetch_server_metadata(session)
    server_name = get_info(session, InfoType.CLI_SERVER_NAME)
    @test length(server_name) > 0
    println("Server name: $server_name")

    cats = catalogs(session)
    @test isa(cats, Tabular)
    println("Catalogs:")
    println(cats)

    sch = schemas(session)
    # at least the default schema should be present
    @test ncells(sch.data) > 1
    println("Schemas:")
    println(sch)

    tbls = tables(session)
    @test isa(tbls, Tabular)
    println("Tables:")
    println(tbls)

    ttypes = tabletypes(session)
    @test isa(ttypes, Tabular)
    println("Table types:")
    println(ttypes)

    cols = columns(session)
    @test isa(cols, Tabular)
    println("Columns:")
    println(cols)

    fns = functions(session, "%")
    # at least the default functions should be present
    @test ncells(fns.data) > 1
    println("Functions:")
    println(fns)
    nothing
end

function create_table_twitter_small(session)
    rs = execute(session, "show tables like 'twitter_small'")
    table_exists = ncells(columnchunk(rs)) > 0

    if table_exists
        println("Use existing table: twitter_small")
    else
        println("Create table: twitter_small")
        result = execute(session, "create table twitter_small (fromid int, toid int) row format delimited fields terminated by ',' lines terminated by '\n' stored as textfile")
        @test round(Int, result) == 0
    end

    rs = execute(session, "select count(*) from twitter_small")
    rowcount = firstcell(columnchunk(rs))
    if rowcount > 0
        println("Use existing data: $rowcount rows")
    else
        rowcount = 10^4
        println("Insert data: $rowcount rows")
        filename = "/tmp/twitter_small"
        open(filename, "w") do io
            writedlm(io, convert(Array{Int}, rand(UInt16, rowcount, 2)), ',')
        end
        result = execute(session, "load data local inpath '$filename' into table twitter_small")
        @test round(Int, result) == 0
    end
end

function create_table_datatype_test(session)
    rs = execute(session, "show tables like 'datatype_test'")
    table_exists = nrows(columnchunk(rs)) > 0
    cols = (
            ("tbool"        , "boolean"         , ()->rand(Bool)),
            ("tint8"        , "tinyint"         , ()->rand(Int8)),
            ("tint16"       , "smallint"        , ()->rand(Int16)),
            ("tint32"       , "int"             , ()->rand(Int32)),
            ("tint64"       , "bigint"          , ()->rand(Int64)),
            ("tfloat32"     , "float"           , ()->rand(Float32)),
            ("tfloat64"     , "double"          , ()->rand(Float64)),
            ("tstr"         , "string"          , ()->randstring()),
            ("tdatetime"    , "timestamp"       , ()->replace(string(now() - Dates.Day(rand(UInt8))), "T"=>" ")),
            ("tdecint32"    , "decimal(9,0)"    , ()->rand(UInt16)),
            ("tdecint64"    , "decimal(18,0)"   , ()->rand(UInt16)),
            ("tdecfloat32"  , "decimal(7,6)"    , ()->rand(Float32)),
            ("tdecfloat64"  , "decimal(16,15)"  , ()->rand(Float64)),
            ("tdecimal"     , "decimal(32,31)"  , ()->(string(rand()) * string(rand(UInt16)))),
            ("tdate"        , "date"            , ()->Date(now() - Dates.Day(rand(UInt8)))),
            ("tchar"        , "char(1)"         , ()->('A' + rand(1:20))),
            ("tchar2"       , "char(2)"         , ()->randstring(2))
    )

    if table_exists
        println("Use existing table: datatype_test")
    else
        println("Create table: datatype_test")
        ct = join(["$(x[1]) $(x[2])" for x in cols], ", ")
        result = execute(session, "create table datatype_test ($ct) row format delimited fields terminated by ',' lines terminated by '\n' stored as textfile")
        @test round(Int, result) == 0
    end

    rs = execute(session, "select count(*) from datatype_test")
    rowcount = firstcell(columnchunk(rs))
    if rowcount > 0
        println("Use existing data: $rowcount rows")
    else
        rowcount = 10^4
        println("Insert data: $rowcount rows")
        filename = "/tmp/datatype_test"
        colvals = hcat([Any[fn() for idx in 1:rowcount] for fn in [x[3] for x in cols]]...)
        open(filename, "w") do io
            writedlm(io, colvals, ',')
        end
        result = execute(session, "load data local inpath '$filename' into table datatype_test")
        @test round(Int, result) == 0
    end
end

function fetch_records(session)
    rs = execute(session, "select min(fromid), max(fromid) from twitter_small")
    minmax = columnchunk(rs)
    maxval = ceil(Int, mean([nthcol(minmax,1)[1], nthcol(minmax,2)[1]]))
    lim = 10000

    println("Execute, record iterator:")
    rs = execute(session, "select * from twitter_small where fromid <= $maxval limit $lim")
    cnt = 0
    for rec in records(rs)
        if rec !== nothing
            (cnt <= 100) && println(rec)
            (cnt == 100) && println("...")
            cnt += 1
        end
    end
    close(rs)
    @test cnt <= lim

    println("Execute, tabular iterator:")
    rs = execute(session, "select * from twitter_small where fromid <= $maxval limit $lim")
    cnt = 0
    for tab in tabulars(rs)
        println(tab)
        cnt += length(nthcol(tab.data,1))
    end
    close(rs)
    @test cnt <= lim

    println("Execute, column chunk iterator:")
    rs = execute(session, "select * from twitter_small where fromid <= $maxval limit $lim")
    cnt = 0
    for colframe in columnchunks(rs)
        for cols in colframe
            println("name  : ", cols[1])
            println("values: ", cols[2][1:min(length(cols[2]), 10)])
            @test typeof(cols[2]) == Vector{Int32}
        end
        cnt += length(nthcol(colframe,1))
    end
    close(rs)
    @test cnt <= lim

    println("Execute, async:")
    rs = execute(session, "select * from twitter_small where fromid <= $maxval limit $lim"; async=true)
    while !isready(rs)
        println("waiting...")
        sleep(10)
    end
    rs = result(rs)
    cc = columnchunk(rs)
    @test nrows(cc) <= lim
    println(Tabular(cc))

    rs = execute(session, "select * from twitter_small where fromid <= $maxval limit $lim"; async=true)
    while !isready(rs)
        println("waiting...")
        sleep(10)
    end
    rs = result(rs)
    cc = columnchunk(rs)
    @test nrows(cc) <= lim
    @test typeof(nthcol(cc,1)) == Vector{Int32}
    @test typeof(nthcol(cc,2)) == Vector{Int32}
    println([(n=>(length(v), typeof(v))) for (n,v) in cc])

    println("Execute, datatypes:")
    rs = execute(session, "select * from datatype_test")
    cols = columnchunk(rs, 100)
    coltypes = [Bool, Int8, Int16, Int32, Int64, Float32, Float64, String, DateTime, Int32, Int64, Float32, Float64, BigFloat, Date, Char, String]
    for ((cn,cv),ct) in zip(cols, coltypes)
        println(cn, " => ", cv[1:min(length(cv),10)])
        @test typeof(cv) == Vector{ct}
        @test length(cv) == 10^4
    end
    close(rs)

    nothing
end

open_database() do session
    create_table_twitter_small(session)
    create_table_datatype_test(session)
    fetch_server_metadata(session)
    fetch_records(session)
end


text:
module TrussIO

using FileIO: FileIO
import FileIO: skipmagic, add_format

using ..TrussTopOptProblems: TrussFEACrossSec, TrussFEAMaterial

using StaticArrays

include("parse_json.jl")
include("parse_geo.jl")

# function load(fn::File{format}) where {format}
#     open(fn) do s
#         skipmagic(s)
#         load(s)
#     end
# end

export load_truss_json, load_truss_geo

end

ext: jl


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# joint computation of gradient and backwards
function get_grad!(grada::AbstractVector{T},
                   gradA::AbstractMatrix{T},
                   grady::AbstractArray{T},
                   betat0::Vector{T}, betat1::Vector{T},
                   Nt::Integer,
                   a::AbstractVector{T},
                   A::AbstractMatrix{T},
                   y::AbstractArray{T},
                   alpha::Matrix{T},
                   c::Vector{T}) where {T <: AbstractFloat}
  Ns = size(A,1)
  for t = Nt-1:-1:1
    betat1 ./= c[t+1]
    for j = 1:Ns
      for k = 1:Ns
        ajk = A[j,k]
        if ajk != zero(T)
          betaA         = ajk * betat1[k]
          betat0[j]    += betaA * y[t+1,k]
          grady[t+1,k] += betaA * alpha[t,j]
          gradA[j,k]   +=  alpha[t,j] * betat1[k] * y[t+1,k]
        end
      end
    end
    betat0, betat1 = betat1, betat0
    fill!(betat0,zero(T))
  end
  for k in 1:Ns
    betat1n    = betat1[k]/c[1]
    grada[k]   = betat1n * y[1,k]
    grady[1,k] = betat1n * a[k]
  end
  grada ./= -Nt
  gradA ./= -Nt
  grady ./= -Nt
  return grada, gradA, grady
end

# with constrained path
function get_grad!(
                   grada,
                   gradA,
                   grady::AbstractArray{T},
                   betat0::Vector{T}, betat1::Vector{T},
                   Nt::Integer,
                   t2IJ::Vector{Pair{Vector{D},Vector{D}}},
                   A::AbstractMatrix{T},
                   y::AbstractArray{T},
                   alpha::Matrix{T},
                   c::Vector{T}) where {D<:Integer,
                                        T <: AbstractFloat}
  Ns = size(A,1)
  for t = Nt-1:-1:1
    betat1 ./= c[t+1]
    I,J = t2IJ[t]
    for k in eachindex(I)
      i,j = I[k],J[k]
      betaA         = A[i,j] * betat1[j]
      betat0[i]    += betaA * y[t+1,j]
      grady[t+1,j] += betaA * alpha[t,i]
      # gradA[i,j]   +=  alpha[t,i] * betat1[j] * y[t+1,j]
    end
    betat0, betat1 = betat1, betat0
    fill!(betat0,zero(T))
  end
  idx0 = unique(t2IJ[1][1])
  for i in idx0
    grady[1,i] += 1/c[1] * betat1[i]
  end
  # gradA ./= -Nt
  grady ./= -Nt
  return grada, gradA, grady
end

function get_grad(Nt::Integer,
                  a,
                  A::AbstractArray{T},
                  y::AbstractArray{T},
                  alpha::Matrix{T},
                  c::Vector{T}) where {T <: AbstractFloat}
  Ns = length(a)
  grada = typeof(a) <: AbstractVector{T} ? zeros(T,Ns) : NoTangent()
  gradA = typeof(a) <: AbstractVector{T} ? zeros(T,Ns,Ns) : NoTangent()
  grady = zeros(T,size(y))
  betat1 = ones(T,Ns)
  betat0 = zeros(T,Ns)
  get_grad!(grada,gradA,grady,betat0,betat1,Nt,a,A,y,alpha,c)
end

# y is a tensor Nt2 × Ns × Nb
function get_grad(Nts::Vector{D},
                  as::AbstractMatrix{T},
                  As::AbstractArray{T,3},
                  ys::AbstractArray{T,3},
                  alphas::Vector{Matrix{T}},
                  cs::Vector{Vector{T}}
                 ) where {T <: AbstractFloat, D <: Integer}
  Ns = size(As,1)
  Nb = size(ys,3)
  grada = zeros(T,Ns,Nb)
  gradA = zeros(T,Ns,Ns,Nb)
  grady = zeros(T,size(ys,1),Ns,Nb)
  betat1 = ones(T,Ns)
  betat0 = zeros(T,Ns)
  for i = 1:Nb
    get_grad!(view(grada,:,i),
              view(gradA,:,:,i),
              view(grady,:,:,i),
              betat0,betat1,
              Nts[i],
              view(as,:,i),
              view(As,:,:,i),
              view(ys,:,:,i),
              alphas[i],cs[i])
    fill!(betat1, one(T))
    fill!(betat0, zero(T))
  end
  return grada,gradA,grady
end

# y is a tensor Nt2 × Ns × Nb constrained path
function get_grad(Nts::Vector{D},
                  as::Vector{K},
                  A::AbstractMatrix{T},
                  y::AbstractArray{T,3},
                  alphas::Vector{Matrix{T}},
                  cs::Vector{Vector{T}}
                 ) where {T <: AbstractFloat, Z,
                          K<:AbstractVector{Z}, D <: Integer}
  Ns = size(A,1)
  Nb = size(y,3)
  grada = NoTangent()
  gradA = NoTangent()
  grady = zeros(T,size(y,1),Ns,Nb)
  betat1 = ones(T,Ns)
  betat0 = zeros(T,Ns)
  for i = 1:Nb
    get_grad!(grada,
              gradA,
              view(grady,:,:,i),
              betat0,betat1,
              Nts[i],as[i],
              A,
              view(y,:,:,i),
              alphas[i],cs[i])
    fill!(betat1, one(T))
    fill!(betat0, zero(T))
  end
  return grada, gradA, grady
end


text:
# Autogenerated wrapper script for rtmidi_jll for x86_64-apple-darwin
export rtmidi

using alsa_jll
JLLWrappers.@generate_wrapper_header("rtmidi")
JLLWrappers.@declare_library_product(rtmidi, "@rpath/librtmidi.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header(alsa_jll)
    JLLWrappers.@init_library_product(
        rtmidi,
        "lib/librtmidi.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
utils.jl
"""

function interp_rows(N::Int,tf::Float64,X::AbstractMatrix)::Matrix
    n,N1 = size(X)
    t1 = range(0,stop=tf,length=N1)
    t2 = collect(range(0,stop=tf,length=N))
    X2 = zeros(n,N)
    for i = 1:n
        interp_cubic = CubicSplineInterpolation(t1, X[i,:])
        X2[i,:] = interp_cubic(t2)
    end
    return X2
end

function convertInf!(A::VecOrMat{Float64},infbnd=1.1e20)
    infs = isinf.(A)
    A[infs] = sign.(A[infs])*infbnd
    return nothing
end

function set_logger()
    if !(global_logger() isa SolverLogger)
        global_logger(SolverLogging.default_logger(true))
    end
end

function benchmark_solve!(solver; samples=10, evals=10)
    Z0 = deepcopy(get_trajectory(solver))
    if is_constrained(solver)
        λ0 = deepcopy(get_duals(solver))
    else
        λ0 = nothing
    end
    v0 = solver.opts.verbose
    s0 = solver.opts.show_summary
    solver.opts.verbose = 0
    solver.opts.show_summary = false
    b = @benchmark begin
        initial_trajectory!($solver,$Z0)
        # set_duals!($solver, $λ0)
        solve!($solver)
    end samples=samples evals=evals
    solver.opts.verbose = v0
    solver.opts.show_summary = s0
    return b
end

function benchmark_solve!(solver, data::Dict; samples=10, evals=10)
	b = benchmark_solve!(solver, samples=samples, evals=evals)

   # Run stats
   push!(data[:time], time(median(b))*1e-6)  # ms
   push!(data[:iterations], iterations(solver))
   push!(data[:cost], cost(solver))
   return b
end

function shift_fill!(A::Vector, n=1, mode=:copy)
    N = length(A)
    N <= 1 && return nothing
	@inbounds for k = n+1:N
		A[k-n] = A[k]
    end
    if mode == :copy
        a_last = A[N-n]
        @inbounds for k = N-n:N
            A[k] = copy(a_last)
        end
    else mode == :zero
        @inbounds for k = N-n:N
            A[k] .= zero(A[k])
        end
    end
	return nothing
end

function ispossemidef(A)
	eigs = eigvals(A)
	if any(real(eigs) .< 0)
		return false
	else
		return true
	end
end

struct NotImplemented <: Exception
	fun::Symbol
	type::Symbol
end

@inline NotImplemented(fun::Symbol, type::DataType) = NotImplemented(fun, Symbol(DataType))
@inline NotImplemented(fun::Symbol, type) = NotImplemented(fun, Symbol(typeof(type)))

Base.showerror(io::IO, ex::NotImplemented) =
	print(io, "Not Implemented Error: ", ex.fun, " not implemented for type ", ex.type)

function gen_zinds(n::Int, m::Int, N::Int, isequal::Bool=false)
	Nu = isequal ? N : N-1
	zinds = [(k-1)*(n+m) .+ (1:n+m) for k = 1:Nu]
	if !isequal
		push!(zinds, (N-1)*(n+m) .+ (1:n))
	end
	return zinds
end


text:
using Test
using Oceananigans
using Oceananigans.AbstractOperations: AbstractOperation
using Oceananigans.Fields: compute_at!
using Oceanostics
using Oceanostics.TKEBudgetTerms
using Oceanostics.TKEBudgetTerms: turbulent_kinetic_energy_ccc
using Oceanostics.FlowDiagnostics
using Oceanostics: SimpleProgressMessenger, SingleLineProgressMessenger

function create_model(; kwargs...)
    model = NonhydrostaticModel(grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1)); kwargs...)
    return model
end

function viscosity(model)
    try
        return model.diffusivity_fields.νₑ
    catch e
        return model.closure.ν
    end
end

include("test_progress_messengers.jl")

function computed_operation(op)
    op_cf = ComputedField(op)
    compute!(op_cf)
    return op_cf
end

function test_vel_only_diagnostics(; model_kwargs...)
    model = create_model(; model_kwargs...)
    test_vel_only_diagnostics(model)
end

function test_vel_only_diagnostics(model)
    u, v, w = model.velocities
    U = AveragedField(u, dims=(1, 2))
    V = AveragedField(v, dims=(1, 2))
    W = AveragedField(w, dims=(1, 2))

    @test begin
        tke_op = KineticEnergy(model)
        ke_c = ComputedField(tke_op)
        compute!(ke_c)
        tke_op isa AbstractOperation
    end

    @test begin
        op = TurbulentKineticEnergy(model, U=U, V=V, W=W)
        tke_c = ComputedField(op)
        compute!(tke_c)
        op isa AbstractOperation
    end

    @test begin
        op = XShearProduction(model, u, v, w, U, V, W)
        XSP = ComputedField(op)
        compute!(XSP)
        op isa AbstractOperation
    end

    @test begin
        op = YShearProduction(model, u, v, w, U, V, W)
        YSP = ComputedField(op)
        compute!(YSP)
        op isa AbstractOperation
    end

    @test begin
        op = ZShearProduction(model, u, v, w, U, V, W)
        ZSP = ComputedField(op)
        compute!(ZSP)
        op isa AbstractOperation
    end

    @test begin
        op = RossbyNumber(model;)
        Ro = ComputedField(op)
        compute!(Ro)
        op isa AbstractOperation
    end

    @test begin
        op = RossbyNumber(model; dUdy_bg=1, dVdx_bg=1, f=1e-4)
        Ro = ComputedField(op)
        compute!(Ro)
        op isa AbstractOperation
    end

    return nothing
end

function test_buoyancy_diagnostics(; model_kwargs...)
    model = create_model(; model_kwargs...)
    test_buoyancy_diagnostics(model)
end

function test_buoyancy_diagnostics(model)
    u, v, w = model.velocities
    b = model.tracers.b
    κ = model.closure.κ.b
    N²₀ = 1e-6

    Ri = RichardsonNumber(model)
    @test Ri isa AbstractOperation

    Ri = RichardsonNumber(model; N²_bg=1, dUdz_bg=1, dVdz_bg=1)
    @test Ri isa AbstractOperation

    PVe = ErtelPotentialVorticity(model)
    @test PVe isa AbstractOperation

    PVtw = ThermalWindPotentialVorticity(model)
    @test PVtw isa AbstractOperation

    PVtw = ThermalWindPotentialVorticity(model, f=1e-4)
    @test PVtw isa AbstractOperation

    return nothing
end

function test_buoyancy_diagnostics(; model_kwargs...)
    model = create_model(; model_kwargs...)
    test_buoyancy_diagnostics(model)
end

function test_pressure_terms(model)
    ∂x_up = XPressureRedistribution(model)
    @test ∂x_up isa AbstractOperation

    ∂y_vp = XPressureRedistribution(model)
    @test ∂y_vp isa AbstractOperation

    ∂z_wp = XPressureRedistribution(model)
    @test ∂z_wp isa AbstractOperation

    return nothing
end

function test_ke_dissipation_rate_terms(; model_kwargs...)
    model = create_model(; model_kwargs...)
    test_ke_dissipation_rate_terms(model)
end

function test_ke_dissipation_rate_terms(model)
    u, v, w = model.velocities
    b = model.tracers.b
    #ν = viscosity(model)

    if model.closure isa AnisotropicDiffusivity
        ε_ani = AnisotropicPseudoViscousDissipationRate(model; U=0, V=0, W=0)
        @test ε_ani isa AbstractOperation

    else
        ε_iso = IsotropicViscousDissipationRate(model; U=0, V=0, W=0)
        @test ε_iso isa AbstractOperation

        ε_iso = IsotropicPseudoViscousDissipationRate(model; U=0, V=0, W=0)
        @test ε_iso isa AbstractOperation
    end

    return nothing
end

function test_tracer_diagnostics(; model_kwargs...)
    model = create_model(; model_kwargs...)
    test_tracer_diagnostics(model)
end

function test_tracer_diagnostics(model)
    u, v, w = model.velocities
    b = model.tracers.b

    if model.closure isa AnisotropicDiffusivity
        κx = model.closure.κx.b
        κy = model.closure.κy.b
        κz = model.closure.κz.b
        χani = AnisotropicTracerVarianceDissipationRate(model, b, κx, κy, κz,)
        @test χani isa AbstractOperation

    else
        κ = model.closure.κ.b
        χiso = IsotropicTracerVarianceDissipationRate(model, b, κ)
        @test χiso isa AbstractOperation
    end

    return nothing
end

@testset "Oceanostics" begin
    model = create_model(; buoyancy=Buoyancy(model=BuoyancyTracer()),
                         coriolis=FPlane(1e-4),
                         tracers=:b,
                         closure=IsotropicDiffusivity(ν=1e-6, κ=1e-7),
                        )

    @info "Testing velocity-only diagnostics"
    test_vel_only_diagnostics(model)

    @info "Testing buoyancy diagnostics"
    test_buoyancy_diagnostics(model)

    @info "Testing pressure terms"
    test_pressure_terms(model)

    closures = (IsotropicDiffusivity(ν=1e-6, κ=1e-7),
                AnisotropicDiffusivity(νz=1e-5, νh=3e-4, κz=1e-6, κh=1e-5),
                SmagorinskyLilly(ν=1e-6, κ=1e-7))
    LESs = (false, false, true)
    messengers = (SingleLineProgressMessenger, TimedProgressMessenger)

    for (LES, closure) in zip(LESs, closures)
        model = create_model(; buoyancy=Buoyancy(model=BuoyancyTracer()),
                             coriolis=FPlane(1e-4),
                             tracers=:b,
                             closure=closure,
                            )

        @info "Testing energy dissipation rate terms with closure" closure
        test_ke_dissipation_rate_terms(model)

        @info "Testing tracer variance terms wth closure" closure
        test_tracer_diagnostics(model)

        @info "Testing SimpleProgressMessenger with closure" closure
        model.clock.iteration = 0
        time_now = time_ns()*1e-9
        test_progress_messenger(model, SimpleProgressMessenger(initial_wall_time_seconds=1e-9*time_ns(), LES=LES))

        @info "Testing SingleLineProgressMessenger with closure" closure
        model.clock.iteration = 0
        time_now = time_ns()*1e-9
        test_progress_messenger(model, SingleLineProgressMessenger(initial_wall_time_seconds=1e-9*time_ns(), LES=LES))

        @info "Testing TimedProgressMessenger with closure" closure
        model.clock.iteration = 0
        time_now = time_ns()*1e-9
        test_progress_messenger(model, TimedProgressMessenger(; LES=LES))
    end

end


text:
# loadElser
@testset "loadElser" begin
    sqrtI, nAtoms, suppSize = FPRS.loadElser("./testData/testData123X")
    @test typeof(sqrtI) <: Array{Float64,2}
    @test nAtoms == 123
end

max_stars_repo_name: tsipenyuk/FPRS

max_issues_repo_name: tsipenyuk/FPRS

max_forks_repo_name: tsipenyuk/FPRS


text:
screen_size = (w = 512, h = 512)

cam = Camera(
    Vec3(0.5f0, 0.2f0, -0.5f0),
    Vec3(0.0f0, 0.1f0,  0.0f0),
    Vec3(0.0f0, 1.0f0,  0.0f0),
    45.0f0,
    0.5f0,
    screen_size...
)

light = PointLight(
    Vec3(1.0f0), 2.0f8,
    Vec3(10.0f0, 10.0f0, 5.0f0)
)

origin, direction = get_primary_rays(cam)

@testset "Image Texture" begin
    scene = load_obj("./meshes/sign_yield.obj")

    @test_nowarn raytrace(origin, direction, scene, light, origin, 0)
end

@testset "Plain Color Texture" begin
    scene = load_obj("./meshes/tree.obj")

    @test_nowarn raytrace(origin, direction, scene, light, origin, 0)

    @inferred raytrace(origin, direction, scene, light, origin, 0)
end


text:
using Test: @test, @testset, @inferred
import Flux
import RestrictedBoltzmannMachines as RBMs
using Random: bitrand
using RestrictedBoltzmannMachines: default_optimizer, minibatch_count, minibatches
using RestrictedBoltzmannMachines: RBM, BinaryRBM, ∂free_energy
using RestrictedBoltzmannMachines: Binary, Spin, Potts, Gaussian, ReLU, dReLU, pReLU, xReLU

@testset "default_optimizer" begin
    nsamples = 60000
    batchsize = 128
    epochs = 100
    decay_final = 0.01
    decay_after = 0.5
    steps_per_epoch = minibatch_count(nsamples; batchsize)
    nsteps = steps_per_epoch * epochs
    startepoch = round(Int, epochs * decay_after)

    o = default_optimizer(
        nsamples, batchsize, epochs;
        optim=Flux.Descent(1), decay_final, decay_after, clip=Inf,
    )
    p = [0.0]
    lr = [only(Flux.Optimise.apply!(o, p, [1.0])) for b in minibatches(nsamples; batchsize) for epoch in 1:epochs]

    decay_gamma = decay_final^(1 / (epochs - startepoch))
    lr_expected = [decay_gamma^max(n - startepoch, 0) for n in 1:epochs]
    lr_expected = repeat(lr_expected; inner=steps_per_epoch)
    # due to the way Flux.ExpDecay is defined, we must shift one step
    @test lr_expected[end] ≈ decay_final
    lr_expected = [lr_expected[2:end]; max(decay_final, decay_gamma * decay_final)]
    @test lr ≈ lr_expected

    @testset "no decay" begin
        decay_after = 2
        startepoch = round(Int, epochs * decay_after)

        o = default_optimizer(
            nsamples, batchsize, epochs;
            optim=Flux.Descent(1), decay_final, decay_after, clip=Inf,
        )
        p = [0.0]
        lr = [only(Flux.Optimise.apply!(o, p, [1.0])) for b in minibatches(nsamples; batchsize) for epoch in 1:epochs]
        @test all(lr .≈ 1)
    end
end

@testset "update!" begin
    opt = Flux.Descent(rand())
    rbm = RBM(Binary(randn(5)), Spin(randn(3)), randn(5,3))
    ∂ = ∂free_energy(rbm, bitrand(5, 10))
    Δ = deepcopy(∂)
    @test RBMs.update!(Δ, rbm, opt) == Δ
    rbm0 = deepcopy(rbm)
    rbm = RBMs.update!(rbm, Δ)
    @test rbm0.visible.θ - rbm.visible.θ ≈ Δ.visible.θ ≈ ∂.visible.θ * opt.eta
    @test rbm0.hidden.θ - rbm.hidden.θ ≈ Δ.hidden.θ ≈ ∂.hidden.θ * opt.eta
    @test rbm0.w - rbm.w ≈ Δ.w ≈ ∂.w * opt.eta

    opt = Flux.Descent(rand())
    layer = Gaussian(randn(5), rand(5))
    ∂ = ∂free_energy(layer, randn(5, 10))
    Δ = deepcopy(∂)
    @test RBMs.update!(Δ, layer, opt) == Δ
    layer0 = deepcopy(layer)
    layer = RBMs.update!(layer, Δ)
    @test layer0.θ - layer.θ ≈ Δ.θ ≈ ∂.θ * opt.eta
    @test layer0.γ - layer.γ ≈ Δ.γ ≈ ∂.γ * opt.eta

    opt = Flux.Descent(rand())
    layer = dReLU(randn(5), randn(5), rand(5), rand(5))
    ∂ = ∂free_energy(layer, randn(5, 10))
    Δ = deepcopy(∂)
    @test RBMs.update!(Δ, layer, opt) == Δ
    layer0 = deepcopy(layer)
    layer = RBMs.update!(layer, Δ)
    @test layer0.θp - layer.θp ≈ Δ.θp ≈ ∂.θp * opt.eta
    @test layer0.θn - layer.θn ≈ Δ.θn ≈ ∂.θn * opt.eta
    @test layer0.γp - layer.γp ≈ Δ.γp ≈ ∂.γp * opt.eta
    @test layer0.γn - layer.γn ≈ Δ.γn ≈ ∂.γn * opt.eta

    opt = Flux.Descent(rand())
    layer = pReLU(randn(5), randn(5), rand(5), rand(5))
    ∂ = ∂free_energy(layer, randn(5, 10))
    Δ = deepcopy(∂)
    @test RBMs.update!(Δ, layer, opt) == Δ
    layer0 = deepcopy(layer)
    layer = RBMs.update!(layer, Δ)
    @test layer0.θ - layer.θ ≈ Δ.θ ≈ ∂.θ * opt.eta
    @test layer0.γ - layer.γ ≈ Δ.γ ≈ ∂.γ * opt.eta
    @test layer0.Δ - layer.Δ ≈ Δ.Δ ≈ ∂.Δ * opt.eta
    @test layer0.η - layer.η ≈ Δ.η ≈ ∂.η * opt.eta

    opt = Flux.Descent(rand())
    layer = xReLU(randn(5), randn(5), rand(5), rand(5))
    ∂ = ∂free_energy(layer, randn(5, 10))
    Δ = deepcopy(∂)
    @test RBMs.update!(Δ, layer, opt) == Δ
    layer0 = deepcopy(layer)
    layer = RBMs.update!(layer, Δ)
    @test layer0.θ - layer.θ ≈ Δ.θ ≈ ∂.θ * opt.eta
    @test layer0.γ - layer.γ ≈ Δ.γ ≈ ∂.γ * opt.eta
    @test layer0.Δ - layer.Δ ≈ Δ.Δ ≈ ∂.Δ * opt.eta
    @test layer0.ξ - layer.ξ ≈ Δ.ξ ≈ ∂.ξ * opt.eta
end

@testset "subtract_gradients" begin
    nt1 = (x = [2], y = [3])
    nt2 = (x = [1], y = [-1])
    @test @inferred(RBMs.subtract_gradients(nt1, nt2)) == (x = [1], y = [4])

    nt1 = (x = [2], y = [3], t = (a = [1], b = [2]))
    nt2 = (x = [1], y = [-1], t = (a = [2], b = [0]))
    @test @inferred(RBMs.subtract_gradients(nt1, nt2)) == (
        x = [1], y = [4], t = (a = [-1], b = [2])
    )
end

@testset "add_gradients" begin
    nt1 = (x = [2], y = [3])
    nt2 = (x = [1], y = [-1])
    @test @inferred(RBMs.add_gradients(nt1, nt2)) == (x = [3], y = [2])

    nt1 = (x = [2], y = [3], t = (a = [1], b = [2]))
    nt2 = (x = [-1], y = [1], t = (a = [-2], b = [0]))
    @test @inferred(RBMs.add_gradients(nt1, nt2)) == (
        x = [1], y = [4], t = (a = [-1], b = [2])
    )
end

@testset "accum_gradients!" begin
    nt1 = (x = [2], y = [3])
    nt2 = (x = [1], y = [-1])
    nt = deepcopy(nt1)
    @test RBMs.add_gradients(nt1, nt2) == @inferred RBMs.accum_gradients!(nt, nt2)
    @test nt == RBMs.add_gradients(nt1, nt2)

    nt1 = (x = [2], y = [3], t = (a = [1], b = [2]))
    nt2 = (x = [-1], y = [1], t = (a = [-2], b = [0]))
    nt = deepcopy(nt1)
    @test RBMs.add_gradients(nt1, nt2) == @inferred RBMs.accum_gradients!(nt, nt2)
    @test nt == RBMs.add_gradients(nt1, nt2)
end

@testset "combine_gradients" begin
    nt1 = (x = [2], y = [3])
    nt2 = (x = [1], y = [-1])
    @test RBMs.combine_gradients(-, nt1, nt2) == (x = [1], y = [4])
    @test RBMs.combine_gradients(-, nt1, nt2) == RBMs.subtract_gradients(nt1, nt2)

    nt1 = (x = [2], y = [3])
    nt2 = (x = [1], y = [-1])
    @test RBMs.combine_gradients(+, nt1, nt2) == (x = [3], y = [2])
    @test RBMs.combine_gradients(+, nt1, nt2) == RBMs.add_gradients(nt1, nt2)

    nt1 = (x = [2], y = [3], t = (a = [1], b = [2]))
    nt2 = (x = [1], y = [-1], t = (a = [2], b = [0]))
    @test RBMs.combine_gradients(-, nt1, nt2) == (
        x = [1], y = [4], t = (a = [-1], b = [2])
    )
    @test RBMs.combine_gradients(+, nt1, nt2) == RBMs.add_gradients(nt1, nt2)
    @test RBMs.combine_gradients(-, nt1, nt2) == RBMs.subtract_gradients(nt1, nt2)
end


text:
#=
day15:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-15
=#

using AoC, Test

neighbourOffsets(dims) = CartesianIndex{dims}[
    [CartesianIndex(ntuple(i->i==j ? 1 : 0, dims)) for j in 1:dims];
    [CartesianIndex(ntuple(i->i==j ? -1 : 0, dims)) for j in 1:dims]
]

function findBestMove(cave)
    risks = fill(length(cave)*max(cave...), size(cave))
    risks[1] = 0

    moves = neighbourOffsets(ndims(cave))

    finished = false
    while !finished
        finished = true
        for I in CartesianIndices(cave)
            for move in moves
                J = I + move
                if checkbounds(Bool, cave, J)
                    risk = risks[I] + cave[J]
                    if risk < risks[J]
                        risks[J] = risk
                        finished = false
                    end
                end
            end
        end
    end

    return risks[end]
end
@test findBestMove([[1,2] [3,4]]) == 6
@test findBestMove([[1,3,1] [1,1,1] [2,1,1]]) == 4

function embiggen(cave)
    realCave = similar(cave, size(cave) .* 5)

    for I in CartesianIndices(cave)
        baseValue = cave[I]
        for j in I[2]:size(cave,2):size(realCave,2)
            value = baseValue
            for i in I[1]:size(cave,1):size(realCave,1)
                realCave[i,j] = value
                value = value < 9 ? value + 1 : 1
            end
            baseValue = baseValue < 9 ? baseValue + 1 : 1
        end
    end
    return realCave
end
@test embiggen(exampleIntMap(15,1))[1:10,1:10] == exampleIntMap(15,1)
@test embiggen(exampleIntMap(15,1)) == exampleIntMap(15,2)

part1(cave) = findBestMove(cave)

part2(cave) = findBestMove(embiggen(cave))

@test part1(exampleIntMap(15,1)) == 40
@test part2(exampleIntMap(15,1)) == 315

@show exampleIntMap(15,1) |> ll -> @time part1(ll)
@show exampleIntMap(15,1) |> ll -> @time part2(ll)

@show intMap(15) |> ll -> @time part1(ll)
@show intMap(15) |> ll -> @time part2(ll)


text:
#=
admm.jl

    Alternating direction method of multipliers (ADMM), also known as
    Douglas-Rachford splitting, to optimize an objective function that can
    be split into two components, where both components can be nonsmooth.

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/01/20
=#

"""
    update_state!(state, λ, α, prox_f!, prox_g!)

Update states using proximal updates with scale `λ` and relaxation determined by
`α`, storing the result in `x`, `z`, and `u`.

#### Arguments
  - `state::ADMMState`  : state variables
  - `λ::Real`           : proximal scaling parameter
  - `α::Real`           : relaxation parameter
  - `prox_f!::Function` : proximal operator of ``f(x)``
  - `prox_g!::Function` : proximal operator of ``g(x)``
"""
function update_state!(state::ADMMState, λ::Real, α::Real, prox_f!::Function, prox_g!::Function)
    # Proximal update x
    state.x.= state.z .- state.u
    prox_f!(state.x, λ)

    # Store previous z
    copyto!(state.z_prev, state.z)

    # Proximal update z
    state.z.= α .* state.x .+ (one(α) - α) .* state.z_prev .+ state.u
    prox_g!(state.z, λ)

    # Update u
    state.u.= state.u .+ α .* state.x .+ (one(α) - α) .* state.z_prev .- state.z

     return nothing
end

"""
    admm(x0, prox_f!, prox_g!; λ=1., α=1., ϵ_abs=1e-7, ϵ_rel=1e-4, max_iter=1000)

Minimize an objective function ``f(x) + g(x)``, where ``f(x)`` and ``g(x)`` can
both be nonsmooth, using alternating direction method of multipliers, also known
as Douglas-Rachford splitting.
The alternating direction method of multipliers method has two hyperparameters,
`λ` and `α`. `λ` controls the scaling of the update steps, i.e. a pseudo step
size, and is equal to the inverse of the augmented Lagrangian parameter. ``α ∈
[0,2]`` is the relaxation parameter, where ``α < 1`` denotes under-relaxation
and ``α > 1`` over-relaxation.

#### Arguments
  - `x0::AbstractVector`: initial parameter values (n x 1)
  - `prox_f!::Function` : proximal operator of ``f(x)``
  - `prox_g!::Function` : proximal operator of ``g(x)``
  - `λ::Real`           : proximal scaling parameter
  - `α::Real`           : relaxation parameter
  - `ϵ_abs::Real`       : absolute tolerance
  - `ϵ_rel::Real`       : relative tolerance
  - `max_iter::Integer` : max number of iterations

#### Returns
  - `x::AbstractVector` : minimizer ``∈ dom f`` (n x 1)
  - `z::AbstractVector` : minimizer ``∈ dom g`` (n x 1)
"""
function admm(x0::AbstractVector, prox_f!::Function, prox_g!::Function;
                λ::Real=1., α::Real=1., ϵ_abs::Real=1e-7, ϵ_rel::Real=1e-4,
                max_iter::Integer=1000)
    # Dimensions
    n= length(x0)

    # Initialize state
    state= ADMMState(similar(x0), copy(x0), zero(x0), similar(x0), similar(x0), similar(x0))

    # Initialize stopping parameters
    iter= 1
    ℓ₂_pri= one(ϵ_abs)
    ℓ₂_dual= one(ϵ_abs)
    ϵ_pri= zero(ϵ_abs)
    ϵ_dual= zero(ϵ_abs)
    # ADMM
    while ℓ₂_pri > ϵ_pri && ℓ₂_dual > ϵ_dual && iter < max_iter
        # Update states
        update_state!(state, λ, α, prox_f!, prox_g!)

        # Primal residual
        state.r.= state.x .+ state.z
        ℓ₂_pri= norm(state.r)

        # Dual residual
        state.s.= inv(λ) .* (state.z .- state.z_prev)
        ℓ₂_dual= norm(state.s)

        # Tolerance
        ϵ_pri= sqrt(n)*ϵ_abs + ϵ_rel*max(norm(state.x), norm(state.z))
        ϵ_dual= sqrt(n)*ϵ_abs + ϵ_rel*norm(state.u)

        # Update iteration counter
        iter+=1
    end

    return (state.x, state.z)
end

"""
    admm!(x, prox_f!, prox_g!; λ=1., α=1., ϵ_abs=1e-7, ϵ_rel=1e-4, max_iter=1000)

Minimize an objective function ``f(x) + g(x)``, where ``f(x)`` and ``g(x)`` can
both be nonsmooth, using alternating direction method of multipliers, also known
as Douglas-Rachford splitting, overwriting `x`. See also `admm`.
"""
function admm!(x::AbstractVector, prox_f!::Function, prox_g!::Function;
                λ::Real=1., α::Real=1., ϵ_abs::Real=1e-7, ϵ_rel::Real=1e-4,
                max_iter::Integer=1000)
    # Dimensions
    n= length(x)

    # Initialize state
    state= ADMMState(x, copy(x), zero(x), similar(x), similar(x), similar(x))

    # Initialize stopping parameters
    iter= 1
    ℓ₂_pri= one(ϵ_abs)
    ℓ₂_dual= one(ϵ_abs)
    ϵ_pri= zero(ϵ_abs)
    ϵ_dual= zero(ϵ_abs)
    # ADMM
    while ℓ₂_pri > ϵ_pri && ℓ₂_dual > ϵ_dual && iter < max_iter
        # Update states
        update_state!(state, λ, α, prox_f!, prox_g!)

        # Primal residual
        state.r.= state.x .+ state.z
        ℓ₂_pri= norm(state.r)

        # Dual residual
        state.s.= inv(λ) .* (state.z .- state.z_prev)
        ℓ₂_dual= norm(state.s)

        # Tolerance
        ϵ_pri= sqrt(n)*ϵ_abs + ϵ_rel*max(norm(state.x), norm(state.z))
        ϵ_dual= sqrt(n)*ϵ_abs + ϵ_rel*norm(state.u)

        # Update iteration counter
        iter+=1
    end

    return state.z
end


text:
# Julia wrapper for header: /usr/local/include/libavutil/aes_ctr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_aes_ctr_alloc,
    av_aes_ctr_init,
    av_aes_ctr_free,
    av_aes_ctr_crypt,
    av_aes_ctr_get_iv,
    av_aes_ctr_set_random_iv,
    av_aes_ctr_set_iv,
    av_aes_ctr_increment_iv

function av_aes_ctr_alloc()
    ccall((:av_aes_ctr_alloc,libavutil),Ptr{AVAESCTR},())
end

function av_aes_ctr_init(a,key)
    ccall((:av_aes_ctr_init,libavutil),Cint,(Ptr{AVAESCTR},Ptr{UInt8}),a,key)
end

function av_aes_ctr_free(a)
    ccall((:av_aes_ctr_free,libavutil),Void,(Ptr{AVAESCTR},),a)
end

function av_aes_ctr_crypt(a,dst,src,size::Integer)
    ccall((:av_aes_ctr_crypt,libavutil),Void,(Ptr{AVAESCTR},Ptr{UInt8},Ptr{UInt8},Cint),a,dst,src,size)
end

function av_aes_ctr_get_iv(a)
    ccall((:av_aes_ctr_get_iv,libavutil),Ptr{UInt8},(Ptr{AVAESCTR},),a)
end

function av_aes_ctr_set_random_iv(a)
    ccall((:av_aes_ctr_set_random_iv,libavutil),Void,(Ptr{AVAESCTR},),a)
end

function av_aes_ctr_set_iv(a,iv)
    ccall((:av_aes_ctr_set_iv,libavutil),Void,(Ptr{AVAESCTR},Ptr{UInt8}),a,iv)
end

function av_aes_ctr_increment_iv(a)
    ccall((:av_aes_ctr_increment_iv,libavutil),Void,(Ptr{AVAESCTR},),a)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module TestGraphvizCategories
using Test

using Catlab, Catlab.Theories, Catlab.CategoricalAlgebra, Catlab.Graphs
using Catlab.Graphs.BasicGraphs: TheoryWeightedGraph
using Catlab.Graphics.GraphvizCategories
using Catlab.Graphics: Graphviz

const stmts = Graphviz.filter_statements

# Categories
############

@present TheorySectionRetract(FreeCategory) begin
  (A, B)::Ob
  i::Hom(A, B)
  r::Hom(B, A)
  compose(i, r) == id(A)
end

gv = to_graphviz(TheorySectionRetract)
@test stmts(gv, Graphviz.Node, :label) == ["A", "B"]
@test stmts(gv, Graphviz.Edge, :label) == ["i", "r"]

# ℳ-categories.
@present TheorySubobject₀(FreeMCategory) begin
  (A, X)::Ob
  ι::Hom(A, X)
end

gv = to_graphviz(TheorySubobject₀, edge_attrs=Dict(:arrowhead => "vee"),
                 tight_attrs=Dict(:dir => "both", :arrowtail => "crow"))
@test stmts(gv, Graphviz.Node, :label) == ["A", "X"]
@test stmts(gv, Graphviz.Edge, :arrowtail) == []

@present TheorySubobject <: TheorySubobject₀ begin
  ::Tight(ι)
end

gv = to_graphviz(TheorySubobject, edge_attrs=Dict(:arrowhead => "vee"),
                 tight_attrs=Dict(:dir => "both", :arrowtail => "crow"))
@test stmts(gv, Graphviz.Edge, :arrowtail) == ["crow"]

# Categories of elements.
el = elements(path_graph(Graph, 2))
gv = to_graphviz(el)
@test stmts(gv, Graphviz.Node, :label) == ["V", "V", "E"]
@test stmts(gv, Graphviz.Edge, :label) == ["src", "tgt"]

gv = to_graphviz(el, node_labels=true, edge_labels=true)
@test stmts(gv, Graphviz.Node, :label) == ["1:V", "2:V", "3:E"]
@test stmts(gv, Graphviz.Edge, :label) == ["1:src", "2:tgt"]

# Schemas
#########

gv = to_graphviz(TheoryWeightedGraph)
@test length(stmts(gv, Graphviz.Node)) == 3
@test stmts(gv, Graphviz.Node, :label) == ["V", "E"]
@test stmts(gv, Graphviz.Node, :xlabel) == ["Weight"]
@test stmts(gv, Graphviz.Edge, :label) == ["src", "tgt", "weight"]

end

max_forks_repo_name: epatters/CompCat


text:
module Yields

import BSplineKit
import ForwardDiff
using LinearAlgebra
using UnicodePlots
using Roots
import LsqFit
import Optim

# don't export type, as the API of Yields.Zero is nicer and
# less polluting than Zero and less/equally verbose as ZeroYieldCurve or ZeroCurve
export rate, discount, accumulation, forward,
LinearSpline, QuadraticSpline, Periodic, Continuous,
Bootstrap,NelsonSiegel,NelsonSiegelSvensson,SmithWilson

include("AbstractYield.jl")
include("Rate.jl")
const DEFAULT_COMPOUNDING = Yields.Continuous()

include("utils.jl")
include("bootstrap.jl")
include("SmithWilson.jl")
include("generics.jl")
include("RateCombination.jl")
include("NelsonSiegelSvensson.jl")

include("precompiles.jl")
end


text:
# JSON Shifts are only available from 2010 onwards
module JSON_Shifts

import JSON

include("shared.jl")

function get_shifts(game_id::AbstractString)
    request = HTTP.request("GET", "http://www.nhl.com/stats/rest/shiftcharts?cayenneExp=gameId=$game_id")
    response = String(request.body)
    return JSON.parse(response)
end

function parse_shift(shift::Dict{String, Any})
    shift_dict = Dict{Symbol, Union{String, Int64}}()
    name = shift["firstName"] * " " * shift["lastName"]
    shift_dict[:player] = name
    shift_dict[:player_id] = shift["playerId"]
    shift_dict[:period] = shift["period"]
    shift_dict[:team] = shift["teamAbbrev"]

    # At the end of the json they list when all the goal events happened. They are the only one"s which have their
    # eventDescription be not null
    if shift["eventDescription"] == nothing
        shift_dict[:start] = convert_to_seconds(shift["startTime"])
        shift_dict[:end] = convert_to_seconds(shift["endTime"])
        shift_dict[:duration] = convert_to_seconds(shift["duration"])
    else
        shift_dict = Dict{Symbol, Union{String, Int64}}()
    end

    return shift_dict
end

function parse_json(shift_json::Dict{String, Any}, game_id::String)
    shifts = [parse_shift(shift) for shift in shift_json["data"]]
    filter!(s -> !isempty(s), shifts)
    df = convert_dict_to_dataframe(shifts)
    df[:game_id] = String(game_id)[6:end]
    sort!(df, (:period, :start), rev=(false, false))
    return df
end

function scrape_shifts(game_id::AbstractString)
    shift_json = get_shifts(game_id)
    return parse_json(shift_json, game_id)
end

end


text:
@testset "test_zpk" begin
# Naming convention:
# ------------------
# {type}_{dims}
# type: C: Continuous, D: Discrete
# dims: "npnuny" (np = # poles)

# CONTINUOUS
C_111 = zpk([-2],[-5],1)
C_211 = zpk([-1+sqrt(2)im,-1-sqrt(2)im], [-5,-3],1)
C_212 = zpk(vecarray(2, 1, [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(2, 1, [-5,-3], [-5,-3]), fill(1, 2, 1))
C_221 = zpk(vecarray(1, 2, [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(1, 2, [-5,-3], [-5,-3]), fill(1, 1, 2))
C_222 = [C_221; C_221]
C_022 = zpk([4 0; 0 4])
s = zpk("s")

# DISCRETE
D_111 = zpk([-2], [0.5], 1, 0.005)
D_211 = zpk([-1+sqrt(2)im,-1-sqrt(2)im], [-0.3, 0.5], 1, 0.005)
D_221 = zpk(vecarray(1, 2, [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(1, 2, [-0.3, 0.5], [-0.3, 0.5]), [1 1], 0.005)
D_222 = [D_221; D_221]
D_022 = zpk([4 0; 0 4], 0.005)
z = zpk("z", 0.005)

# TESTS
# Constructors
@test s == zpk([0], Int64[], 1)
@test z == zpk([0], Int64[], 1, 0.005)
@test C_022 == zpk(vecarray(Int64, 2, 2, Int64[], Int64[], Int64[], Int64[]), vecarray(Int64, 2, 2, Int64[], Int64[], Int64[], Int64[]), [4 0; 0 4])
@test D_022 == zpk(vecarray(2, 2, Int64[], Int64[], Int64[], Int64[]), vecarray(2, 2, Int64[], Int64[], Int64[], Int64[]), [4 0; 0 4], 0.005)
@test C_022 == [zpk(4) 0;0 4]
#TODO We might want to fix this
@test_broken D_022 == [zpk(4, 0.005) 0;0 4]

# Test constructors with empty matrices of type Any
@test zpk([-1.0], [], 1.0) == zpk([-1.0], Float64[], 1.0)
@test zpk([], [-1.0], 1.0) == zpk(Float64[], [-1.0], 1.0)
@test zpk([], [], 1.0) == zpk(Float64[], Float64[], 1.0)
@test zpk([], [1.0+im,1.0-im], 1.0) == zpk(ComplexF64[], [1.0+im,1.0-im], 1.0)

#TODO improve polynomial accuracy se these are equal
# Addition
@test C_111 + C_111 ≈ zpk([-2,-5], [-5,-5], 2)
@test C_222 + C_222 ≈
  zpk(vecarray(2, 2, [-1+sqrt(2)im,-1-sqrt(2)im], [-2], [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(2, 2, [-5,-3], [-5,-3], [-5,-3], [-5,-3]), [2 2;2 2])
z1 = [-2.5+sqrt(9-2.5^2)im, -2.5-sqrt(9-2.5^2)im]
z2 = [-4.5-sqrt(4.5^2-17), -4.5+sqrt(4.5^2-17)]

@test C_222 + 1 ≈ zpk(vecarray(2, 2, z1, z2, z1, z2), vecarray(2, 2, [-3, -5], [-3, -5], [-3, -5], [-3, -5]), [2 1; 2 1])
@test D_111 + D_111 ≈ zpk([-2],[0.5],2,0.005)

# Subtraction
@test C_111 - C_211 ≈ zpk([-1], [-5,-3.0], 3.0)
@test minreal(zpk(tf([0,3.0,18,15], [1,13,55,75])), 1e-6) ≈ zpk([-1.0], [-5,-3.0], 3.0)

@test 1 - C_222 ≈ zpk(tf(vecarray(2, 2, [0,6,12], [1,7,13], [0,6,12], [1,7,13]),
  vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15])))
@test D_111 - D_211 - zpk(tf([0,0.3,-2.55,1.2], [1,-0.7,-0.05,0.075], 0.005)) ≈
  zpk(tf([0], [1], 0.005))

# Multiplication
@test C_111 * C_221 ≈ zpk(tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]),
  vecarray(1, 2, [1,13,55,75], [1,13,55,75])))
@test C_212 * C_111 ≈ zpk(tf(vecarray(2, 1, [1,4,7,6], [0,1,4,4]),
  vecarray(2, 1, [1,13,55,75], [1,13,55,75])))
@test 4*C_222 ≈ zpk(tf(vecarray(2, 2, [4,8,12], [0,4,8], [4,8,12], [0,4,8]),
  vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15])))
@test D_111 * D_221 - zpk(tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]),
  vecarray(1, 2, [1,-0.7,-0.05,0.075], [1.0,-0.7,-0.05,0.075]), 0.005)) ≈
  zpk(tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005))

# Division
@test 1/C_111 ≈ zpk(tf([1,5], [1,2]))
@test C_212/C_111 ≈ zpk(tf(vecarray(2, 1, [1,7,13,15], [0,1,7,10]),
  vecarray(2, 1, [1,10,31,30], [1,10,31,30])))
@test 1/D_111 ≈ zpk(tf([1.0,-0.5], [1.0,2.0], 0.005))

# Indexing
@test size(C_222) == (2, 2)
@test size(C_212) == (2, 1)
@test C_222[1,1] ≈ zpk(tf([1, 2, 3], [1, 8, 15]))
@test C_222[1:1,1] ≈ zpk(tf([1, 2, 3], [1, 8, 15]))
@test C_222[1,1:2] == C_221
@test size(C_222[1,Int64[]]) == (1,0)

# Test that number of poles matter
@test !(zpk(Int64[],[1,1],1) == zpk(Int64[],[1],1))

# Test SispZpk and SisoTf operations
C_212_tf = tf(C_212)
C_111_tf = tf(C_111)
# Add
@test C_212_tf + C_212 ≈ 2*C_212
@test C_212 + C_212_tf ≈ 2*C_212
# Approx
@test C_212_tf + C_212 ≈ 2*C_212_tf
# Minus
@test 2*C_212_tf - C_212 ≈ C_212
# Multiply
@test C_212_tf*C_111 ≈ C_212*C_111
@test C_212_tf*C_111 ≈ C_212_tf*C_111_tf

# TODO test printing when it is implemented better

# Tests of minreal
@test minreal(zpk([-5.0], [-5.0, -5.0], 1.0)) == zpk(Float64[], [-5.0], 1.0)
@test minreal(zpk([-1.0, -2.0], [-2.0, -2.0], 1.0)) == zpk([-1.0], [-2.0], 1.0)

# Tests of minreal (roudning and cancelling complex roots with real roots)
@test minreal(zpk([-2.0+1e-10im,-2.0-1e-10im], [-2.0], 1.0)) == zpk([-2.0], Float64[], 1.0)
@test minreal(zpk([-2.0], [-2.0+1e-10im,-2.0-1e-10im], 1.0)) == zpk(Float64[], [-2.0], 1.0)
@test minreal(zpk([-1.0, -2.0, -2.0+1e-10im,-2.0-1e-10im], [-2.0, -2.0], 1.0)) == zpk([-1.0, -2.0], Float64[], 1.0)
@test minreal(zpk([-1.0, -2.0+1e-10im,-2.0-1e-10im,-2.0], [-2.0, -2.0], 1.0)) == zpk([-1.0, -2.0], Float64[], 1.0)
@test minreal(zpk([-2.0, -2.0], [-1.0, -2.0, -2.0+1e-10im,-2.0-1e-10im], 1.0)) ==  zpk(Float64[], [-1.0, -2.0], 1.0)
@test minreal(zpk([-2.0, -2.0], [-1.0, -2.0+1e-10im,-2.0-1e-10im,-2.0], 1.0)) == zpk(Float64[], [-1.0, -2.0], 1.0)

# Test of minreal (systems with no poles / no zeros)
@test minreal(zpk([-1.0, -2.0], Float64[], 2.5)) == zpk([-1.0, -2.0], Float64[], 2.5)
@test minreal(zpk(Float64[], [-1.0, -2.0], 2.5)) == zpk(Float64[], [-1.0, -2.0], 2.5)

# Type stability Continuous and discrete time
# QUESTION: why not use typeof instead of eltype(fill())
@test eltype(fill(zpk("s"),2)) <: TransferFunction
@test eltype(fill(zpk([1],[1,1],1),2)) <: TransferFunction
@test eltype(fill(zpk(Int64[],[1,1],1.0),2)) <: TransferFunction
@test eltype(fill(zpk([1 2; 3 4],0.001),2)) <: TransferFunction
@test eltype(fill(zpk(1)+zpk(2),2)) <: TransferFunction
@test eltype(fill(zpk(1,0.005)/zpk(2, 0.005),2)) <: TransferFunction
@test eltype(fill(zpk(1)+1,2)) <: TransferFunction

@test eltype([tf(1,1), zpk(1,1)]) <: TransferFunction

zpk(tf([1 2; 3 4])) == zpk([1 2; 3 4])

# Errors
@test_throws ErrorException C_111 + C_222             # Dimension mismatch
@test_throws ErrorException C_111 - C_222             # Dimension mismatch
@test_throws ErrorException C_111 * C_222             # Dimension mismatch
@test_throws ErrorException [s 0; 1]                  # Dimension mismatch
@test_throws ErrorException D_111 + C_111             # Sampling time mismatch
@test_throws ErrorException D_111 - C_111             # Sampling time mismatch
@test_throws ErrorException D_111 * C_111             # Sampling time mismatch
D_diffTs = zpk(tf([1], [2], 0.1))
@test_throws ErrorException D_111 + D_diffTs          # Sampling time mismatch
@test_throws ErrorException D_111 - D_diffTs          # Sampling time mismatch
@test_throws ErrorException D_111 * D_diffTs          # Sampling time mismatch
@test_throws ErrorException zpk([1], [2], 1, -0.1)        # Negative samping time
@test_throws ErrorException zpk("s", 0.01)             # s creation can't be discrete
@test_throws ErrorException zpk("z", 0)                # z creation can't be continuous
@test_throws ErrorException zpk("z")                   # z creation can't be continuous
# Remove this when inferec is implemented
@test_throws ErrorException [z 0]                     # Sampling time mismatch (inferec could be implemented)

@test typeof(zpk(tf([1], [2], 0.1))) == TransferFunction{ControlSystems.SisoZpk{Float64,Complex{Float64}}}
@test typeof(zpk([-0.5], [], 1)) == TransferFunction{ControlSystems.SisoZpk{Float64,Float64}}
end


text:
module SampleData

using Celeste: Model, DeterministicVI
import Celeste: Infer

import Synthetic

using Distributions
using StaticArrays
import WCS, FITSIO, DataFrames

export empty_model_params, dat_dir,
       sample_ce, perturb_params,
       sample_star_fluxes, sample_galaxy_fluxes,
       gen_sample_star_dataset, gen_sample_galaxy_dataset,
       gen_two_body_dataset, gen_three_body_dataset, gen_n_body_dataset,
       make_elbo_args, true_star_init

const dat_dir = joinpath(Pkg.dir("Celeste"), "test", "data")

const sample_star_fluxes = [
    4.451805E+03,1.491065E+03,2.264545E+03,2.027004E+03,1.846822E+04]
const sample_galaxy_fluxes = [
    1.377666E+01, 5.635334E+01, 1.258656E+02,
    1.884264E+02, 2.351820E+02] * 100  # 1x wasn't bright enough

# A world coordinate system where the world and pixel coordinates are the same.
const wcs_id = WCS.WCSTransform(2,
                    cd = Float64[1 0; 0 1],
                    ctype = ["none", "none"],
                    crpix = Float64[1, 1],
                    crval = Float64[1, 1]);

"""
Turn a images and vector of catalog entries into elbo arguments
that can be used with Celeste.
"""
function make_elbo_args(images::Vector{Image},
                        catalog::Vector{CatalogEntry};
                        active_source=-1,
                        patch_radius_pix::Float64=NaN)
    vp = Vector{Float64}[init_source(ce) for ce in catalog]
    patches = Infer.get_sky_patches(images,
                                    catalog,
                                    radius_override_pix=patch_radius_pix)
    S = length(catalog)
    active_sources = active_source > 0 ? [active_source] :
                                          S <= 3 ? collect(1:S) : [1,2,3]
    ElboArgs(images, vp, patches, active_sources)
end

"""
Load a stamp into a Celeste images.
"""
function load_stamp_blob(stamp_dir, stamp_id)
    function fetch_image(b)
        band_letter = band_letters[b]
        filename = "$stamp_dir/stamp-$band_letter-$stamp_id.fits"

        fits = FITSIO.FITS(filename)
        hdr = FITSIO.read_header(fits[1])
        original_pixels = read(fits[1])
        dn = original_pixels / hdr["CALIB"] + hdr["SKY"]
        nelec_f32 = round.(dn * hdr["GAIN"])
        nelec = convert(Array{Float64}, nelec_f32)

        header_str = FITSIO.read_header(fits[1], String)
        wcs = WCS.from_header(header_str)[1]
        close(fits)

        alphaBar = [hdr["PSF_P0"]; hdr["PSF_P1"]; hdr["PSF_P2"]]
        xiBar = [
            [hdr["PSF_P3"]  hdr["PSF_P4"]];
            [hdr["PSF_P5"]  hdr["PSF_P6"]];
            [hdr["PSF_P7"]  hdr["PSF_P8"]]]'

        tauBar = Array(Float64, 2, 2, 3)
        tauBar[:,:,1] = [[hdr["PSF_P9"] hdr["PSF_P11"]];
                         [hdr["PSF_P11"] hdr["PSF_P10"]]]
        tauBar[:,:,2] = [[hdr["PSF_P12"] hdr["PSF_P14"]];
                         [hdr["PSF_P14"] hdr["PSF_P13"]]]
        tauBar[:,:,3] = [[hdr["PSF_P15"] hdr["PSF_P17"]];
                         [hdr["PSF_P17"] hdr["PSF_P16"]]]

        psf = [PsfComponent(alphaBar[k], SVector{2,Float64}(xiBar[:, k]),
                            SMatrix{2,2,Float64,4}(tauBar[:, :, k])) for k in 1:3]

        H, W = size(original_pixels)
        iota = hdr["GAIN"] / hdr["CALIB"]
        epsilon = hdr["SKY"] * hdr["CALIB"]

        run_num = round(Int, hdr["RUN"])
        camcol_num = round(Int, hdr["CAMCOL"])
        field_num = round(Int, hdr["FIELD"])

        epsilon_mat = fill(epsilon, H, W)
        iota_vec = fill(iota, H)
        empty_psf_comp = RawPSF(Array(Float64, 0, 0), 0, 0,
                                 Array(Float64, 0, 0, 0))

        Image(H, W, nelec, b, wcs, psf,
              run_num, camcol_num, field_num, epsilon_mat, iota_vec,
              empty_psf_comp)
    end

    images = map(fetch_image, 1:5)
end

function load_stamp_catalog_df(cat_dir, stamp_id, images; match_blob=false)
    # These files are generated by
    # https://github.com/dstndstn/tractor/blob/master/projects/inference/testblob2.py
    cat_fits = FITSIO.FITS("$cat_dir/cat-$stamp_id.fits")
    num_cols = FITSIO.read_key(cat_fits[2], "TFIELDS")[1]
    ttypes = [FITSIO.read_key(cat_fits[2], "TTYPE$i")[1] for i in 1:num_cols]

    df = DataFrames.DataFrame()
    for i in 1:num_cols
        tmp_data = read(cat_fits[2], ttypes[i])
        df[Symbol(ttypes[i])] = tmp_data
    end

    close(cat_fits)

    if match_blob
        camcol_matches = df[:camcol] .== images[3].camcol_num
        run_matches = df[:run] .== images[3].run_num
        field_matches = df[:field] .== images[3].field_num
        df = df[camcol_matches & run_matches & field_matches, :]
    end

    df
end

"""
Load a stamp catalog.
"""
function load_stamp_catalog(cat_dir, stamp_id, images; match_blob=false)
    df = load_stamp_catalog_df(cat_dir, stamp_id, images,
                                    match_blob=match_blob)
    df[:objid] = [ string(s) for s=1:size(df)[1] ]

    function row_to_ce(row)
        x_y = [row[1, :ra], row[1, :dec]]
        star_fluxes = zeros(5)
        gal_fluxes = zeros(5)
        fracs_dev = [row[1, :frac_dev], 1 - row[1, :frac_dev]]
        for b in 1:length(band_letters)
            bl = band_letters[b]
            psf_col = Symbol("psfflux_$bl")

            # TODO: How can there be negative fluxes?
            star_fluxes[b] = max(row[1, psf_col], 1e-6)

            dev_col = Symbol("devflux_$bl")
            exp_col = Symbol("expflux_$bl")
            gal_fluxes[b] += fracs_dev[1] * max(row[1, dev_col], 1e-6) +
                             fracs_dev[2] * max(row[1, exp_col], 1e-6)
        end

        fits_ab = fracs_dev[1] > .5 ? row[1, :ab_dev] : row[1, :ab_exp]
        fits_phi = fracs_dev[1] > .5 ? row[1, :phi_dev] : row[1, :phi_exp]
        fits_theta = fracs_dev[1] > .5 ? row[1, :theta_dev] : row[1,
:theta_exp]

        # tractor defines phi as -1 * the phi catalog for some reason.
        if !match_blob
            fits_phi *= -1.
        end

        re_arcsec = max(fits_theta, 1. / 30)  # re = effective radius
        re_pixel = re_arcsec / 0.396

        phi90 = 90 - fits_phi
        phi90 -= floor(phi90 / 180) * 180
        phi90 *= (pi / 180)

        CatalogEntry(x_y, row[1, :is_star], star_fluxes,
            gal_fluxes, row[1, :frac_dev], fits_ab, phi90, re_pixel,
            row[1, :objid], 0)
    end

    CatalogEntry[row_to_ce(df[i, :]) for i in 1:size(df, 1)]
end

function empty_model_params(S::Int)
    vp = [Model.init_source([ 0., 0. ]) for s in 1:S]
    ElboArgs(Image[],
             vp,
             Array(SkyPatch, S, 0),
             collect(1:S))
end

function sample_ce(pos, is_star::Bool)
    CatalogEntry(pos, is_star, sample_star_fluxes, sample_galaxy_fluxes,
        0.1, .7, pi/4, 4., "sample", 0)
end

function perturb_params(ea) # for testing derivatives != 0
    for vs in ea.vp
        vs[ids.a[:, 1]] = [ 0.4, 0.6 ]
        vs[ids.u[1]] += .8
        vs[ids.u[2]] -= .7
        vs[ids.r1] -= log(10)
        vs[ids.r2] *= 25.
        vs[ids.e_dev] += 0.05
        vs[ids.e_axis] += 0.05
        vs[ids.e_angle] += pi/10
        vs[ids.e_scale] *= 1.2
        vs[ids.c1] += 0.5
        vs[ids.c2] =  1e-1
    end
end

function gen_sample_star_dataset(; perturb=true)
    srand(1)
    images0 = load_stamp_blob(dat_dir, "164.4311-39.0359_2kpsf")
    for b in 1:5
        images0[b].H, images0[b].W = 20, 23
        images0[b].wcs = wcs_id
    end
    one_body = [sample_ce([10.1, 12.2], true),]
    images = Synthetic.gen_blob(images0, one_body)
    ea = make_elbo_args(images, one_body)
    if perturb
        perturb_params(ea)
    end
    images, ea, one_body
end

function gen_sample_galaxy_dataset(; perturb=true)
    srand(1)
    images0 = load_stamp_blob(dat_dir, "164.4311-39.0359_2kpsf")
    for b in 1:5
        images0[b].H, images0[b].W = 20, 23
        images0[b].wcs = wcs_id
    end
    one_body = [sample_ce([8.5, 9.6], false),]
    images = Synthetic.gen_blob(images0, one_body)
    ea = make_elbo_args(images, one_body)
    if perturb
        perturb_params(ea)
    end
    images, ea, one_body
end

function gen_two_body_dataset(; perturb=true)
    # A small two-body dataset for quick unit testing.  These objects
    # will be too close to be identifiable.

    srand(1)
    images0 = load_stamp_blob(dat_dir, "164.4311-39.0359_2kpsf")
    for b in 1:5
        images0[b].H, images0[b].W = 20, 23
        images0[b].wcs = wcs_id
    end
    two_bodies = [
        sample_ce([4.5, 3.6], false),
        sample_ce([10.1, 12.1], true)
    ]
    images = Synthetic.gen_blob(images0, two_bodies)
    ea = make_elbo_args(images, two_bodies)
    if perturb
        perturb_params(ea)
    end
    images, ea, two_bodies
end

function gen_three_body_dataset(; perturb=true)
    srand(1)
    images0 = load_stamp_blob(dat_dir, "164.4311-39.0359_2kpsf")
    for b in 1:5
        images0[b].H, images0[b].W = 112, 238
        images0[b].wcs = wcs_id
    end
    three_bodies = [
        sample_ce([4.5, 3.6], false),
        sample_ce([60.1, 82.2], true),
        sample_ce([71.3, 100.4], false),
    ];
    images = Synthetic.gen_blob(images0, three_bodies);
    ea = make_elbo_args(images, three_bodies);
    if perturb
        perturb_params(ea)
    end
    images, ea, three_bodies
end

"""
Generate a large dataset with S randomly placed bodies and non-constant
background.
"""
function gen_n_body_dataset(
    S::Int; patch_pixel_radius=20., seed=NaN)

  if !isnan(seed)
    srand(seed)
  end

  images0 = load_stamp_blob(dat_dir, "164.4311-39.0359_2kpsf");
  img_size_h = 900
  img_size_w = 1000
  for b in 1:5
      images0[b].H, images0[b].W = img_size_h, img_size_w
  end

  fluxes = [4.451805E+03,1.491065E+03,2.264545E+03,2.027004E+03,1.846822E+04]

  locations = rand(2, S) .* [img_size_h, img_size_w]
  world_locations = WCS.pix_to_world(images0[3].wcs, locations)

  S_bodies = CatalogEntry[CatalogEntry(world_locations[:, s], true,
      fluxes, fluxes, 0.1, .7, pi/4, 4., string(s), s) for s in 1:S];

  images = Synthetic.gen_blob(images0, S_bodies);

  # Make non-constant background.
  for b=1:5
    images[b].iota_vec = fill(images[b].iota_vec[1], images[b].H)
    images[b].epsilon_mat = fill(images[b].epsilon_mat[1], images[b].H, images[b].W)
  end

  ea = make_elbo_args(
    images, S_bodies, patch_radius_pix=patch_pixel_radius)

  images, ea, S_bodies
end

function true_star_init()
    images, ea, body = gen_sample_star_dataset(perturb=false)

    ea.vp[1][ids.a[:, 1]] = [ 1.0 - 1e-4, 1e-4 ]
    ea.vp[1][ids.r2] = 1e-4
    ea.vp[1][ids.r1] = log(sample_star_fluxes[3]) - 0.5 * ea.vp[1][ids.r2]
    #ea.vp[1][ids.r1] = sample_star_fluxes[3] ./ ea.vp[1][ids.r2]
    ea.vp[1][ids.c2] = 1e-4

    images, ea, body
end

end # End module


text:
@testset "Vector Spherical Harmonic mode indices" begin
    lmax = 10
    N = lmax + 1
    M = 2 * N - 1

    for v in 1:2
        mode_seen = falses(lmax + 1, 2 * lmax + 1)
        for l in 1:lmax, m in (-l):l
            lm = sphv_mode(l, m, v)
            @test !mode_seen[lm]
            mode_seen[lm] = true
        end
        @test all(mode_seen .≤ sphvones(Bool, N, M))
    end
end

@testset "Vector Spherical Harmonics for $T" for T in
                                                 [Float64, Complex{Float64}]
    # lmax = 100
    lmax = 4
    atol = rtol = 4 / lmax^2

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    # Test duality
    F = randn(T, N, M)
    C = sphv_transform(F)
    @test size(C) == size(F)
    @test eltype(C) == eltype(F)
    F′ = sphv_evaluate(C)
    @test size(F′) == size(F)
    @test eltype(F′) == eltype(F)
    C′ = sphv_transform(F′)
    @test C′ ≈ C
    F″ = sphv_evaluate(C′)
    @test F″ ≈ F′

    # Test gradient field, l=1, m=0
    Fθ = T[grad_Y_1_0_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[grad_Y_1_0_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test size(Fθ) == (N, M)
    @test size(Fϕ) == (N, M)
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 1 * 2; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    @test Cθ[sphv_mode(1, 0, 1)] ≈ sqrt(2)
    @test sum(abs2.(Cθ)) ≈ 2
    @test sum(abs2.(Cϕ)) ≈ 0
    Cgrad = Cθ                  # l=0

    # Test curl field, l=1, m=0
    Fθ = T[curl_Y_1_0_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[curl_Y_1_0_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 1 * 2; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    @test sum(abs2.(Cθ)) ≈ 0
    @test Cϕ[sphv_mode(1, 0, 2)] ≈ -sqrt(2)
    @test sum(abs2.(Cϕ)) ≈ 2
    Ccurl = Cϕ                  # l=0

    # Test gradient field, l=1, m=1
    Fθ = T[grad_Y_1p1_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[grad_Y_1p1_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 1 * 2; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    @test Cθ[sphv_mode(1, 1, 1)] ≈ sqrt(1 / 2)
    @test Cϕ[sphv_mode(1, 1, 2)] ≈ sqrt(3 / 2)
    @test sum(abs2.(Cθ)) ≈ 1 / 2
    @test sum(abs2.(Cϕ)) ≈ 3 / 2

    # Test curl field, l=1, m=1
    Fθ = T[curl_Y_1p1_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[curl_Y_1p1_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 1 * 2; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    # TODO: Using converse `v` value
    @test Cθ[sphv_mode(1, 1, 2)] ≈ sqrt(3 / 2)
    @test Cϕ[sphv_mode(1, 1, 1)] ≈ -sqrt(1 / 2)
    @test sum(abs2.(Cθ)) ≈ 3 / 2
    @test sum(abs2.(Cϕ)) ≈ 1 / 2

    # Test gradient field, l=1, m=-1
    Fθ = T[grad_Y_1m1_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[grad_Y_1m1_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 1 * 2; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    @test Cθ[sphv_mode(1, -1, 1)] ≈ sqrt(1 / 2)
    @test Cϕ[sphv_mode(1, -1, 2)] ≈ sqrt(3 / 2)
    @test sum(abs2.(Cθ)) ≈ 1 / 2
    @test sum(abs2.(Cϕ)) ≈ 3 / 2

    # Test curl field, l=1, m=-1
    Fθ = T[curl_Y_1m1_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[curl_Y_1m1_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 1 * 2; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    # TODO: Using converse `v` value
    @test Cθ[sphv_mode(1, -1, 2)] ≈ sqrt(3 / 2)
    @test Cϕ[sphv_mode(1, -1, 1)] ≈ -sqrt(1 / 2)
    @test sum(abs2.(Cθ)) ≈ 3 / 2
    @test sum(abs2.(Cϕ)) ≈ 1 / 2

    # Test gradient field, l=2, m=0
    Fθ = T[grad_Y_2_0_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[grad_Y_2_0_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 2 * 3; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    @test Cθ[sphv_mode(2, 0, 1)] ≈ sqrt(6)
    @test sum(abs2.(Cθ)) ≈ 6
    @test sum(abs2.(Cϕ)) ≈ 0

    # Test gradient field, l=2, m=1
    Fθ = T[grad_Y_2p1_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[grad_Y_2p1_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 2 * 3; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    @test Cθ[sphv_mode(2, 1, 1)] ≈ sqrt(8 / 3)
    @test Cϕ[sphv_mode(2, 1, 2)] ≈ sqrt(10 / 4)
    @show chop.(Cθ)
    @show chop.(Cϕ)
    @test sum(abs2.(Cθ)) ≈ 8 / 3
    @test sum(abs2.(Cϕ)) ≈ 10 / 4
    @assert false

    # Test gradient field, l=3, m=0
    Fθ = T[grad_Y_3_0_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[grad_Y_3_0_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 3 * 4; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    @test Cθ[sphv_mode(3, 0, 1)] ≈ sqrt(12)
    @test sum(abs2.(Cθ)) ≈ 12
    @test sum(abs2.(Cϕ)) ≈ 0

    # Test gradient field, l=3, m=1
    Fθ = T[grad_Y_3p1_θ(θ, ϕ) for θ in Θ, ϕ in Φ]
    Fϕ = T[grad_Y_3p1_ϕ(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test isapprox(integrate(Fθ, Fϕ, Fθ, Fϕ), 3 * 4; atol=atol, rtol=rtol)
    Cθ = sphv_transform(Fθ)
    Cϕ = sphv_transform(Fϕ)
    @test Cθ[sphv_mode(3, 1, 1)] ≈ sqrt(27 / 4)
    @test sum(abs2.(Cθ)) ≈ 27 / 4
    @test sum(abs2.(Cϕ)) ≈ 0

    # Test orthogonality
    lmax_test = 4
    for gc in 1:2, l in 1:lmax_test, m in (-l):l
        Cθ = zeros(T, N, M)
        Cϕ = zeros(T, N, M)
        if gc == 1
            Cθ[sphv_mode(l, m, 1)] = sqrt(1 / 2)
            Cϕ[sphv_mode(l, m, 2)] = sqrt(3 / 2)
        else
            Cθ[sphv_mode(l, m, 2)] = sqrt(3 / 2)
            Cϕ[sphv_mode(l, m, 1)] = -sqrt(1 / 2)
        end
        Fθ = sphv_evaluate(Cθ)
        Fϕ = sphv_evaluate(Cϕ)
        for gc′ in 1:2, l′ in 1:lmax_test, m′ in (-l′):l′
            gc′ == gc || continue
            Cθ′ = zeros(T, N, M)
            Cϕ′ = zeros(T, N, M)
            if gc′ == 1
                Cθ′[sphv_mode(l′, m′, 1)] = sqrt(1 / 2)
                Cϕ′[sphv_mode(l′, m′, 2)] = sqrt(3 / 2)
            else
                Cθ′[sphv_mode(l′, m′, 2)] = sqrt(3 / 2)
                Cϕ′[sphv_mode(l′, m′, 1)] = -sqrt(1 / 2)
            end
            Fθ′ = sphv_evaluate(Cθ′)
            Fϕ′ = sphv_evaluate(Cϕ′)
            int = integrate(Fθ, Fϕ, Fθ′, Fϕ′)
            δ = l * (l + 1) * (gc == gc′ && l == l′ && m == m′)
            if !isapprox(int, δ; atol=atol, rtol=rtol)
                @show gc gc′ l l′ m m′ int δ
                @assert false
            end
            @test isapprox(int, δ; atol=atol, rtol=rtol)
        end
    end
end


text:
export NFSPAgentManager

## definition
"""
    NFSPAgentManager(; agents::Dict{Any, NFSPAgent})

A special MultiAgentManager in which all agents use NFSP policy to play the game.
"""
mutable struct NFSPAgentManager <: AbstractPolicy
    agents::Dict{Any, NFSPAgent}
end

## interactions between the policy and env.
function (π::NFSPAgentManager)(env::AbstractEnv)
    player = current_player(env)
    if player == chance_player(env)
        env |> legal_action_space |> rand
    else
        env |> π.agents[player]
    end
end

RLBase.prob(π::NFSPAgentManager, env::AbstractEnv, args...) = prob(π.agents[current_player(env)], env, args...)

## update NFSPAgentManager
function RLBase.update!(π::NFSPAgentManager, env::AbstractEnv)
    while current_player(env) == chance_player(env)
        env |> legal_action_space |> rand |> env
    end
    update!(π.agents[current_player(env)], env)
end

function (π::NFSPAgentManager)(stage::Union{PreEpisodeStage, PostEpisodeStage}, env::AbstractEnv)
    @sync for (player, agent) in π.agents
        @async agent(stage, env, player)
    end
end


text:
#=
# 3D heat eq. Dirichlet bc, CG
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("heat3d");

# Try making an optional log
@useLog("heat3dlog")

n = 6;
ord = 2;

# Set up the configuration (order doesn't matter)
@domain(3, SQUARE, UNIFORM_GRID)    # dimension, geometry, decomposition
@solver(CG)                         # DG, CG, etc.
@functionSpace(LEGENDRE, ord)         # function, order (or use testFunction and trialFunction)
@nodes(LOBATTO)                     # elemental node arrangement
@stepper(EULER_IMPLICIT)            # time stepper (optional second arg is CFL#)

# Specify the problem
@mesh(HEXMESH, n)                   # .msh file or generate our own

@variable(u)                        # same as @variable(u, SCALAR)

@testSymbol(v)                    # sets the symbol for a test function

T = 1;
@timeInterval(T)                    # (start, end) using this sets problem to time dependent
@initial(u, "abs(x-0.5)+abs(y-0.5)+abs(z-0.5) < 0.2 ? 1 : 0")  # initial condition needed if time dependent

@boundary(u, 1, DIRICHLET, 0)

# Write the weak form
@coefficient(f, "-0.1*sin(pi*x)*sin(pi*y)*sin(pi*z)")
@weakForm(u, "Dt(u*v) + 0.01 * dot(grad(u),grad(v)) - f*v")

solve(u);

# solution is stored in the variable's "values"
#using Plots
#pyplot();
#N = n*ord+1;
#half = Int(round(N/2));
#range = (N*N*half+1):(N*N*(half+1));
#display(plot(Finch.grid_data.allnodes[1,range], Finch.grid_data.allnodes[2,range], u.values[range], st=:surface))

# check
log_dump_config(Finch.config);
log_dump_prob(Finch.prob);

@finalize()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 04-2022
# Functions for spherical symmetry
#--------------------------------------------------------------------

"""
 Compute initial data for η, given a and ψ
"""
function computeη(U::NTuple{3,Field{S}})::Field{S} where {S<:Space{Tag}} where {Tag}
    PS = first(U).space
    D = derivative(PS)
    I = identity(PS)
    (a, η, ϕ) = U
    B = incomingboundary(PS) + outgoingboundary(PS)
    A = D*D - (2/a)*(D*a)*D + 4pi*(D*ϕ)^2*I
    η = linsolve((I - B)*A + B, B*η)
    return η
end

"""
    Compute the Ricci scalar
"""
function R(U::NTuple{3, Field{S}})::Field{S} where {S}
    (a, η, ϕ) = U
    DU, DV = derivative(first(U).space)
    return - (8 * (DU * ϕ) * (DV * ϕ)) / a^2
end

"""
    Compute the Hawking mass
"""
function M(U::NTuple{3, Field{S}})::Field{S} where {S}
    (a, η, ϕ) = U
    DU, DV = derivative(first(U).space)
    return (η / 2) * (1 + (4 * (DU * η ) * (DV * η)) / η^2)
end

"""
    Compute the expansion
    TODO: Check if the sign is right.
"""
function expansion(U::NTuple{3, Field{S}})::Field{S} where {S}
    (a, η, ϕ) = U
    DU, DV = derivative(first(U).space)
    return (DV * η)
end

function constraints(U::NTuple{3, Field})::NTuple{1, Field}
    PS = first(U).space
    DU, DV = derivative(PS)
    (a, η, ϕ) = U
    C1 = DU*(DU*η) - (2/a)*(DU*a)*(DU*η) + (4pi*η)*(DU*ϕ)^2
    C2 = DV*(DV*η) - (2/a)*(DV*a)*(DV*η) + (4pi*η)*(DV*ϕ)^2
    return (sqrt(C1*C1 + C2*C2), )
end

function lineconstraints(U::NTuple{3, Field})::Field
    PS = first(U).space
    DU = derivative(PS)
    (a, η, ϕ) = U
    return DU*(DU*η) - (2/a)*(DU*a)*(DU*η) + (4pi*η)*(DU*ϕ)^2
end

function residuals(U::NTuple{3, Field})::NTuple{1, Field}
    PS = first(U).space
    DU, DV = derivative(PS)
    (a, η, ϕ) = U
    F1 = DU*(DV*a) - (1/a)*(DU*a)*(DV*a) + (a/η)*(DU*(DV*η)) + (4pi*a)*(DU*ϕ)*(DV*ϕ)
    F2 = DU*(DV*η) + (1/η)*(DU*η)*(DV*η) + (1/4)*(1/η)*(a^2)
    F3 = DU*(DV*ϕ) + (1/η)*(DU*η)*(DV*ϕ) + (1/η)*(DV*η)*(DU*ϕ)
    return (sqrt(F1*F1 + F2*F2 + F3*F3), )
end


text:
module Missile_Types

using Aerospace

export
    Controls,
    SoftwareMsg,
    Aerodynamics,
    Missile

type Controls
    throttle::Float64
    alpha::Float64
    alphadot::Float64
    phi::Float64
    phidot::Float64
end

type SoftwareMsg
    time::Float64
    dt::Float64
    SWmode::UInt
    controls::Controls
    navstates::NavStates
    wypt_num::UInt
    axcmd::Float64
    aycmd::Float64
    azcmd::Float64
    eta_dive_cmd::Float64
    eta_turn_cmd::Float64
    tgo::Float64
    miss::Float64
    missvec_ned::Array{Float64,1}
    vel_err_int::Float64;
end

type Missile
    time::Float64
    dt::Float64
    airframe::Aero5DOF
    control::Controls
    software::SoftwareMsg
end

type Aerodynamics
    sref::Float64
    cl::Float64
    cd::Float64
    faero_body::Array{Float64,1}
end

end ## End of Module


text:
# R: number of regions
# G: number of generators (either thermal / VG / storage)
# T: number of operations timesteps per period
# P: number of operations periods (representative)

mutable struct InvestmentProblem{R,G1,G2,G3,I,T,P} <: AbstractProblem{R,G1,G2,G3,I,T,P}

    model::Model
    regionnames::Vector{String}
    periodnames::Vector{String}
    technologies::Technologies{G1,G2,G3,I,R}
    initialconditions::InitialConditions{R,G1,G2,G3}
    discountrate::Float64
    rootscenario::Scenario{R,G1,G2,G3,I,T,P}

    function InvestmentProblem{T}(
        model::Model,
        regionnames::Vector{String},
        periodnames::Vector{String},
        techs::Technologies{G1,G2,G3,I,R},
        initconds::InitialConditions{R,G1,G2,G3},
        discountrate::Float64
    ) where {R,G1,G2,G3,I,T}

        @assert R == length(regionnames)
        P = length(periodnames)

        new{R,G1,G2,G3,I,T,P}(model, regionnames, periodnames,
                              techs, initconds, discountrate)

    end

end

function InvestmentProblem(
    regionnames::Vector{String},
    periodnames::Vector{String},
    techs::Technologies{G1,G2,G3,I,R},
    initconds::InitialConditions{R,G1,G2,G3},
    discountrate::Float64,
    rootname::String,
    investments::Investments{R,G1,G2,G3},
    operations::Operations{R,G1,G2,G3,I,T,P},
    markets::Markets{R,T,P},
    periodweights::Vector{Float64},
    optimizer::Optimizer
) where {R,G1,G2,G3,I,T,P}

    invprob = InvestmentProblem{T}(
        Model(optimizer), regionnames, periodnames,
        techs, initconds, discountrate)

    root = Scenario(invprob, rootname,
                    investments, operations, markets, periodweights)
    invprob.rootscenario = root

    return invprob

end

# For now, loading in single-scenario only
function InvestmentProblem(
    folder::String, optimizer::Optimizer
)

    # Load representative periods
    perioddata = DataFrame!(CSV.File(joinpath(folder, "periods.csv"),
                                     types=period_types))
    n_periods = length(perioddata.name)
    periodlookup = Dict(zip(perioddata.name, 1:n_periods))

    technicalpath = joinpath(folder, "technical")

    n_timesteps = 24 # TODO: Load these from a file
    discountrate = 0.95

    # Load regions

    regiondata = DataFrame!(CSV.File(joinpath(technicalpath, "regions.csv"),
                                     type=String))
    n_regions = length(regiondata.name)
    regionlookup = Dict(zip(regiondata.name, 1:n_regions))

    # Load technologies

    thermaldata = DataFrame!(CSV.File(joinpath(technicalpath, "thermal.csv")))
    thermaltechs = ThermalGenerators(
        thermaldata.name, thermaldata.owner,
        thermaldata.mingen, thermaldata.maxgen,
        thermaldata.minuptime, thermaldata.mindowntime,
        thermaldata.maxrampup, thermaldata.maxrampdown,
        thermaldata.capacitycredit)
    n_thermals = length(thermaltechs.name)
    thermallookup = Dict(zip(thermaldata.name, 1:n_thermals))

    variabledata = DataFrame!(CSV.File(joinpath(technicalpath, "variable.csv")))
    variabletechs = VariableGenerators(
        variabledata.name, variabledata.owner, variabledata.maxgen,
        variabledata.capacitycredit)
    n_variables = length(variabletechs.name)
    variablelookup = Dict(zip(variabledata.name, 1:n_variables))

    storagedata = DataFrame!(CSV.File(joinpath(technicalpath, "storage.csv")))
    storagetechs = StorageDevices(
        storagedata.name, storagedata.owner,
        storagedata.maxcharge, storagedata.maxdischarge, storagedata.maxenergy,
        storagedata.chargeefficiency, storagedata.dischargeefficiency,
        storagedata.carryoverefficiency, storagedata.capacitycredit)
    n_storages = length(storagetechs.name)
    storagelookup = Dict(zip(storagedata.name, 1:n_storages))

    interfacedata = DataFrame!(CSV.File(
        joinpath(technicalpath, "interface.csv"), type=String))

    from_idx = getindex.(Ref(regionlookup), interfacedata.from)
    to_idx = getindex.(Ref(regionlookup), interfacedata.to)
    interfaces = Interfaces{n_regions}(
        interfacedata.name, tuple.(from_idx, to_idx))
    n_interfaces = length(interfacedata.name)

    techs = Technologies(thermaltechs, variabletechs, storagetechs, interfaces)

    # Load starting builds / options

    initialdata =
        DataFrame!(CSV.File(joinpath(folder, "initialconditions.csv"),
                            types=initialcondition_types))

    thermalstarts = InitialInvestments(
        zeros(Int, n_regions, n_thermals), zeros(Int, n_regions, n_thermals))

    variablestarts = InitialInvestments(
        zeros(Int, n_regions, n_variables), zeros(Int, n_regions, n_variables))

    storagestarts = InitialInvestments(
        zeros(Int, n_regions, n_storages), zeros(Int, n_regions, n_storages))

    for r in eachrow(initialdata)

        region_idx = regionlookup[r.region]

        thermal_idx = get(thermallookup, r.class, nothing)
        if !isnothing(thermal_idx)
            thermalstarts.options[region_idx, thermal_idx] = r.optioned
            thermalstarts.builds[region_idx, thermal_idx] = r.built
            continue
        end

        variable_idx = get(variablelookup, r.class, nothing)
        if !isnothing(variable_idx)
            variablestarts.options[region_idx, variable_idx] = r.optioned
            variablestarts.builds[region_idx, variable_idx] = r.built
            continue
        end

        storage_idx = get(storagelookup, r.class, nothing)
        if !isnothing(storage_idx)
            storagestarts.options[region_idx, storage_idx] = r.optioned
            storagestarts.builds[region_idx, storage_idx] = r.built
            continue
        end

        error("Unknown class $(r.class) listed in initial conditions")

    end

    initconds = InitialConditions(thermalstarts, variablestarts, storagestarts)

    # Load scenario data

    scenariospath = joinpath(folder, "scenarios")
    scenariosdata = DataFrame!(CSV.File(joinpath(scenariospath, "tree.csv"),
                                        types=scenarios_tree_types))

    root_idx = findfirst(ismissing, scenariosdata.parent)
    rootname = scenariosdata.name[root_idx]
    rootpath = joinpath(scenariospath, rootname)
    investments, operations, markets =
        loadscenario(regionlookup, techs, periodlookup, n_timesteps, rootpath)

    result = InvestmentProblem(
        regiondata.name, perioddata.name, techs, initconds, discountrate,
        rootname, investments, operations, markets,
        perioddata.weight, optimizer)

    function makechildren!(parent::Scenario)

        for (childname, prob) in getchildren(parent.name, scenariosdata)

            investments, operations, markets = loadscenario(
                regionlookup, techs, periodlookup, n_timesteps,
                joinpath(scenariospath, childname))

            child = Scenario(parent, childname, prob,
                             investments, operations, markets,
                             perioddata.weight)

            makechildren!(child)

        end

    end

    makechildren!(result.rootscenario)

    return result

end

function solve!(invprob::InvestmentProblem; debug::Bool=false)
    @objective(invprob.model, Max, welfare(invprob.rootscenario))
    debug && open(f -> println(f, invprob.model), "problem.txt", "w")
    optimize!(invprob.model)
    return
end

getchildren(scenario::String, data::DataFrame) = [(row.name, row.probability)
    for row in eachrow(data) if !ismissing(row.parent) && row.parent == scenario]

scenarios(invprob::InvestmentProblem) = scenarios(invprob.rootscenario)

function scenarios(root::Scenario)
    iszero(length(root.children)) && return [(root, true)]
    return vcat((root, false), scenarios.(root.children)...)
end


text:
"""
    readLOC(fname::AbstractString)

Import a leadfield location matrix in the BESA loc format

#### Arguments

* `fname` = location and name of file to import
* `verbose` = should the function output user feedback

#### Returns

locations = readLOC(verbose=true)
"""
function readLOC(fname::AbstractString = Pkg.dir("Private", "test", "data", "FEM_Default50MRI.loc");
    verbose::Bool=false)

    debug("Location file: $(fname)")

    # Open file
    fid = open(fname, "r")

    # Read header information from file
    version         = read(fid, Int32, 1)[1]
    numNodes        = read(fid, Int32, 1)[1]

    # Read in positions
    numPositions = 3*numNodes
    positions = read(fid, Float64, numPositions)

    # Read in number of neighbours
    numNeighbours = read(fid, Int32, numNodes)

    maxNeighbours = maximum(numNeighbours)
    neighbours = zeros(Int32, (convert(Int64, numNodes), convert(Int64,maxNeighbours)))

    i = 1
    while i <= numNodes
        j = 1
        while j <= numNeighbours[i]
            neighbours[i,j]= read(fid, Int32, 1)[1]
            j = j + 1
        end
        i = i + 1
    end

    # Check that everything has been read
    if ~eof(fid)
        println("!! Something went wrong. File is longer than expected")
    end

    close(fid)

    # Reshape position data
    positions = reshape(vec(positions), (Int(3), Int(numNodes)))
    positions = positions*1e3           # Convert to mm
    positions = positions .- 127.5
    positions = reshape(positions, (3, Int(numNodes)))

    if verbose
        println("LFT version:            $(version)")
        println("Number of nodes:        $(numNodes)")
        println("Number of positions:    $(numPositions)")
        println("Maximum # neighbours:   $(maxNeighbours)")

    end

    return positions, neighbours
end


text:
# Autogenerated wrapper script for tectonic_jll for x86_64-w64-mingw32-cxx11
export tectonic

using Fontconfig_jll
using FreeType2_jll
using Graphite2_jll
using HarfBuzz_jll
using ICU_jll
using OpenSSL_jll
using Zlib_jll
using libpng_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "PATH"
LIBPATH_default = ""

# Relative path to `tectonic`
const tectonic_splitpath = ["bin", "tectonic.exe"]

# This will be filled out by __init__() for all products, as it must be done at runtime
tectonic_path = ""

# tectonic-specific global declaration
function tectonic(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ';', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ';', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(tectonic_path)
    end
end

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    # Initialize PATH and LIBPATH environment variable listings
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Fontconfig_jll.PATH_list, FreeType2_jll.PATH_list, Graphite2_jll.PATH_list, HarfBuzz_jll.PATH_list, ICU_jll.PATH_list, OpenSSL_jll.PATH_list, Zlib_jll.PATH_list, libpng_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Fontconfig_jll.LIBPATH_list, FreeType2_jll.LIBPATH_list, Graphite2_jll.LIBPATH_list, HarfBuzz_jll.LIBPATH_list, ICU_jll.LIBPATH_list, OpenSSL_jll.LIBPATH_list, Zlib_jll.LIBPATH_list, libpng_jll.LIBPATH_list,))

    global tectonic_path = normpath(joinpath(artifact_dir, tectonic_splitpath...))

    push!(PATH_list, dirname(tectonic_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ';')
    global LIBPATH = join(vcat(LIBPATH_list, [Sys.BINDIR, joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ';')

end  # __init__()


text:
function (==)(x::TimeDate, y::TimeDate)
     Date(x) === Date(y) && Time(x) === Time(y)
end
function (!=)(x::TimeDate, y::TimeDate)
     Date(x) !== Date(y) || (Date(x) == Date(y) && Time(x) != Time(y))
end
function (<)(x::TimeDate, y::TimeDate)
    dtx = DateTime(x)
    dty = DateTime(y)
    (dtx < dty) || ((dtx === dty) && (Time(x) < Time(y)))
end
function (>)(x::TimeDate, y::TimeDate)
    dtx = DateTime(x)
    dty = DateTime(y)
    (dtx > dty) || ((dtx === dty) && (Time(x) > Time(y)))
end
function (<=)(x::TimeDate, y::TimeDate)
    dtx = DateTime(x)
    dty = DateTime(y)
    (dtx < dty) || ((dtx === dty) && (Time(x) <= Time(y)))
end
function (>=)(x::TimeDate, y::TimeDate)
    dtx = DateTime(x)
    dty = DateTime(y)
    (dtx > dty) || ((dtx === dty) && (Time(x) >= Time(y)))
end

isless(x::TimeDate, y::DateTime) = isless(x, TimeDate(y))
isless(x::DateTime, y::TimeDate) = isless(TimeDate(x), y)

function isequal(x::TimeDateZone, y::TimeDateZone)
    xx = Microsecond(x) + Nanosecond(x)
    yy = Microsecond(y) + Nanosecond(y)
    zx = ZonedDateTime(x)
    zy = ZonedDateTime(y)
    isequal(zx, zy) && isequal(xx, yy)
end

function isless(x::TimeDateZone, y::TimeDateZone)
    xx = Microsecond(x) + Nanosecond(x)
    yy = Microsecond(y) + Nanosecond(y)
    zx = ZonedDateTime(x)
    zy = ZonedDateTime(y)
    isless(zx, zy) || (isequal(zx,zy) && isless(xx, yy))
end

(<)(x::TimeDateZone, y::TimeDateZone) = isless(x,y)
(>)(x::TimeDateZone, y::TimeDateZone) = isless(y,x)
(<=)(x::TimeDateZone, y::TimeDateZone) = isless(x,y) || isequal(x,y)
(>=)(x::TimeDateZone, y::TimeDateZone) = isless(y,x) || isequal(x,y)
(==)(x::TimeDateZone, y::TimeDateZone) = isequal(x,y)
(!=)(x::TimeDateZone, y::TimeDateZone) = !isequal(x,y)


text:
using ModelingToolkit, LinearAlgebra, SparseArrays

# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const _DD = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 8
const X = reshape([i for i in 1:N for j in 1:N],N,N)
const Y = reshape([j for i in 1:N for j in 1:N],N,N)
const α₁ = 1.0.*(X.>=4*N/5)

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

# Define the initial condition as normal arrays
@variables du[1:N,1:N,1:3] u[1:N,1:N,1:3] MyA[1:N,1:N] AMx[1:N,1:N] DA[1:N,1:N]

# Define the discretized PDE as an ODE function
function f(du,u,p,t)
   A = @view  u[:,:,1]
   B = @view  u[:,:,2]
   C = @view  u[:,:,3]
  dA = @view du[:,:,1]
  dB = @view du[:,:,2]
  dC = @view du[:,:,3]
  mul!(MyA,My,A)
  mul!(AMx,A,Mx)
  @. DA = _DD*(MyA + AMx)
  @. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
  @. dB = α₂ - β₂*B - r₁*A*B + r₂*C
  @. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

f(du,u,nothing,0.0)

multithreadedf = eval(ModelingToolkit.build_function(du,u,multithread=true)[2])
_du = rand(N,N,3)
_u = rand(N,N,3)
multithreadedf(_du,_u)

jac = sparse(ModelingToolkit.jacobian(vec(du),vec(u),simplify=false))
multithreadedjac = eval(ModelingToolkit.build_function(vec(jac),u,multithread=true)[2])

#_jac = similar(jac,Float64)
#multithreadedjac(_jac,_u)


text:
names = [
    "test_fock.jl"
    "test_nlevel.jl"
    "test_meanfield.jl"
    "test_parameters.jl"
    "test_average.jl"
    "test_diffeq.jl"
    "test_v-level.jl"
    "test_mixed-order.jl"
    "test_correlation.jl"
    "test_two-level-laser.jl"
    "test_cluster.jl"
    "test_scaling.jl"
    "test_higher-order.jl"
]

detected_tests = filter(
    name->startswith(name, "test_") && endswith(name, ".jl"),
    readdir("."))

unused_tests = setdiff(detected_tests, names)
if length(unused_tests) != 0
    @warn string("The following tests are not used:\n", join(unused_tests, "\n"))
end

unavailable_tests = setdiff(names, detected_tests)
if length(unavailable_tests) != 0
    error("The following tests could not be found:\n", join(unavailable_tests, "\n"))
end

for name=names
    if startswith(name, "test_") && endswith(name, ".jl")
        include(name)
    end
end


text:
# [[file:../fresa.org::testsimple][testsimple]]
# load in the Fresa optimization package
using Fresa
# specify the dimension of the search space
nx = 2
# create an initial point in the search space
x0 = 0.5*ones(nx)
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x)))
# the actual objective function
f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0.0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,f)]
# now invoke Fresa to solve the problem
@time best, pop = Fresa.solve(f, p0, domain)
# output the results
println("Population at end: $pop")
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")
# testsimple ends here

# [[file:../fresa.org::*simple objective function][simple objective function:2]]
#+plot: ind:1 deps:(2) with:"linespoints lt 3 pt 7 ps 0.25" set:nokey set:"yrange [0:1]" set:"xrange [0:*]"
println("#+plot: ind:1 deps:(2) with:\"linespoints pt 7 ps 0.25\" set:nokey set:\"yrange [0:1]\" set:\"xrange [0:*]\" set:\"xlabel 'Generation'\" set:\"ylabel 'fitness'\"")
Fresa.printHistoryTrace(best)
# simple objective function:2 ends here


text:
# These will get encoded into each dictionary when a struct is serialized.
const METADATA_KEY = "__metadata__"
const TYPE_KEY = "type"
const MODULE_KEY = "module"
const PARAMETERS_KEY = "parameters"
const CONSTRUCT_WITH_PARAMETERS_KEY = "construct_with_parameters"
const FUNCTION_KEY = "function"

"""
Serializes a InfrastructureSystemsType to a JSON file.
"""
function to_json(
    obj::T,
    filename::AbstractString;
    force = false,
) where {T <: InfrastructureSystemsType}
    if !force && isfile(filename)
        error("$file already exists. Set force=true to overwrite.")
    end
    result = open(filename, "w") do io
        return to_json(io, obj)
    end

    @info "Serialized $T to $filename"
    return result
end

"""
Serializes a InfrastructureSystemsType to a JSON string.
"""
function to_json(obj::T)::String where {T <: InfrastructureSystemsType}
    return JSON3.write(serialize(obj))
end

"""
Deserializes a InfrastructureSystemsType from a JSON filename.
"""
function from_json(::Type{T}, filename::String) where {T <: InfrastructureSystemsType}
    return open(filename) do io
        from_json(io, T)
    end
end

"""
Deserializes a InfrastructureSystemsType from String or IO.
"""
function from_json(io::Union{IO, String}, ::Type{T}) where {T <: InfrastructureSystemsType}
    return deserialize(T, JSON3.read(io, Dict))
end

"""
Serialize the Julia value into standard types that can be converted to non-Julia formats,
such as JSON. In cases where val is an instance of a struct, return a Dict. In cases where
val is a scalar value, return that value.
"""
function serialize(val::T) where {T <: InfrastructureSystemsType}
    @debug "serialize InfrastructureSystemsType" _group = LOG_GROUP_SERIALIZATION val T
    return serialize_struct(val)
end

function serialize(vals::Vector{T}) where {T <: InfrastructureSystemsType}
    @debug "serialize Vector{InfrastructureSystemsType}" _group = LOG_GROUP_SERIALIZATION vals T
    return serialize_struct.(vals)
end

function serialize(func::Function)
    return Dict{String, Any}(
        METADATA_KEY => Dict{String, Any}(
            FUNCTION_KEY => string(nameof(func)),
            MODULE_KEY => string(parentmodule(func)),
        ),
    )
end

function serialize_struct(val::T) where {T}
    @debug "serialize_struct" _group = LOG_GROUP_SERIALIZATION val T
    data = Dict(string(name) => serialize(getfield(val, name)) for name in fieldnames(T))
    add_serialization_metadata!(data, T)
    return data
end

"""
Add type information to the dictionary that can be used to deserialize the value.
"""
function add_serialization_metadata!(data::Dict, ::Type{T}) where {T}
    data[METADATA_KEY] = Dict{String, Any}(
        TYPE_KEY => string(nameof(T)),
        MODULE_KEY => string(parentmodule(T)),
    )
    if !isempty(T.parameters)
        data[METADATA_KEY][PARAMETERS_KEY] = string.(T.parameters)
    end

    return
end

"""
Return the type information for the serialized struct.
"""
get_serialization_metadata(data::Dict) = data[METADATA_KEY]

function get_type_from_serialization_data(data::Dict)
    return get_type_from_serialization_metadata(get_serialization_metadata(data))
end

function get_type_from_serialization_metadata(metadata::Dict)
    _module = get_module(metadata[MODULE_KEY])
    base_type = getfield(_module, Symbol(metadata[TYPE_KEY]))
    if !get(metadata, CONSTRUCT_WITH_PARAMETERS_KEY, false)
        return base_type
    end

    # This has several limitations and is only a workaround for PSY.Reserve subtypes.
    # - each parameter must be in _module
    # - does not support nested parametrics.
    # Reserves should be fixed and then we can remove this hack.
    parameters = [getfield(_module, Symbol(x)) for x in metadata[PARAMETERS_KEY]]
    return base_type{parameters...}
end

serialize(val::Base.RefValue{T}) where {T} = serialize(val[])

# The default implementation allows any scalar type (or collection of scalar types) to
# work. The JSON library must be able to encode and decode anything passed here.

serialize(val::T) where {T} = val

"""
Deserialize an object from standard types stored in non-Julia formats, such as JSON, into
Julia types.
"""
function deserialize(::Type{T}, data::Dict) where {T <: InfrastructureSystemsType}
    @debug "deserialize InfrastructureSystemsType" _group = LOG_GROUP_SERIALIZATION T data
    return deserialize_struct(T, data)
end

function deserialize_struct(::Type{TimeSeriesKey}, data::Dict)
    vals = Dict{Symbol, Any}()
    for (field_name, field_type) in
        zip(fieldnames(TimeSeriesKey), fieldtypes(TimeSeriesKey))
        val = data[string(field_name)]
        if field_name == :time_series_type
            val = getfield(InfrastructureSystems, Symbol(strip_module_name(val)))
        end
        vals[field_name] = val
    end
    return TimeSeriesKey(; vals...)
end

function deserialize_to_dict(::Type{T}, data::Dict) where {T}
    # Note: mostly duplicated in src/deterministic_metadata.jl
    vals = Dict{Symbol, Any}()
    for (field_name, field_type) in zip(fieldnames(T), fieldtypes(T))
        val = data[string(field_name)]
        if val isa Dict && haskey(val, METADATA_KEY)
            metadata = get_serialization_metadata(val)
            if haskey(metadata, FUNCTION_KEY)
                vals[field_name] = deserialize(Function, val)
            else
                vals[field_name] =
                    deserialize(get_type_from_serialization_metadata(metadata), val)
            end
        else
            vals[field_name] = deserialize(field_type, val)
        end
    end
    return vals
end

function deserialize_struct(::Type{T}, data::Dict) where {T}
    vals = deserialize_to_dict(T, data)
    return T(; vals...)
end

function deserialize(::Type{Function}, data::Dict)
    metadata = data[METADATA_KEY]
    return get_type_from_strings(metadata[MODULE_KEY], metadata[FUNCTION_KEY])
end

function deserialize(::Type{T}, data::Any) where {T}
    @debug "deserialize Any" _group = LOG_GROUP_SERIALIZATION T data
    return data
end

function deserialize(::Type{T}, data::Array) where {T <: Tuple}
    return tuple(data...)
end

function deserialize(::Type{T}, data::Any) where {T <: AbstractFloat}
    return T(data)
end

function deserialize(::Type{T}, data::Dict) where {T <: NamedTuple}
    return T(key = data[string(key)] for key in fieldnames(T))
end

function deserialize(
    ::Type{T},
    data::Union{Nothing, Dict},
) where {T <: Union{Nothing, NamedTuple}}
    return isnothing(data) ? nothing : deserialize(T.b, data)
end

function deserialize(::Type{T}, data::Array) where {T <: Vector{<:Tuple}}
    return [tuple(x...) for x in data]
end

# Enables JSON serialization of Dates.Period.
# The default implementation fails because the field is defined as abstract.
# Encode the type when serializing so that the correct value can be deserialized.
function serialize(resolution::Dates.Period)
    return Dict(
        "value" => resolution.value,
        TYPE_KEY => strip_module_name(typeof(resolution)),
    )
end

function deserialize(::Type{Dates.Period}, data::Dict)
    return getfield(Dates, Symbol(data[TYPE_KEY]))(data["value"])
end

deserialize(::Type{Dates.DateTime}, val::AbstractString) = Dates.DateTime(val)

# The next methods fix serialization of UUIDs. The underlying type of a UUID is a UInt128.
# JSON tries to encode this as a number in JSON. Encoding integers greater than can
# be stored in a signed 64-bit integer sometimes does not work - at least when using
# JSON3. The number gets converted to a float in scientific notation, and so
# the UUID is truncated and essentially lost. These functions cause JSON to encode UUIDs as
# strings and then convert them back during deserialization.

serialize(uuid::Base.UUID) = Dict("value" => string(uuid))
serialize(uuids::Vector{Base.UUID}) = serialize.(uuids)
deserialize(::Type{Base.UUID}, data::Dict) = Base.UUID(data["value"])

serialize(value::Complex) = Dict("real" => real(value), "imag" => imag(value))
deserialize(::Type{Complex}, data::Dict) = Complex(data["real"], data["imag"])
deserialize(::Type{Complex{T}}, data::Dict) where {T} =
    Complex(T(data["real"]), T(data["imag"]))

deserialize(::Type{Vector{Symbol}}, data::Vector) = Symbol.(data)

function serialize_julia_info()
    data = Dict{String, Any}("julia_version" => string(VERSION))
    io = IOBuffer()
    Pkg.status(io = io, mode = Pkg.PKGMODE_MANIFEST)
    data["package_info"] = String(take!(io))
    return data
end