text:
# println("Testing vec2state! and state2vec!...")
# srand(1)
# P = 20
# K = 10
# Sigma_X = spdiagm(rand(P) .+ 0.2)
# A = sparse(1:P, Int64[ceil(i/2) for i in 1:P], ones(P))
# Sigma_L = FactorAnalysis.randcor(K, 0.2)
# d = CFADistribution(Sigma_X, A, Sigma_L)
# N = 1000
# X = rand(d, N)
# S = X*X' / N
# Base.cov2cor!(S, sqrt(diag(S)))
# function ll_test(d, S, N)
#     Theta = inv(d.Sigma_X + d.A*d.Sigma_L*d.A')
#     logdet(Theta) - trace(S*Theta)
# end
# @test abs(loglikelihood(d, S, N) - ll_test(d, S, N)) < 1e-8

println("Testing CFADistributionSigma: vec2state! and state2vec!...")
srand(10)
P = 20
K = 10
Sigma_X = spdiagm(rand(P) .+ 0.2)
A = sparse(1:P, Int64[ceil(i/2) for i in 1:P], ones(P))
Sigma_L = FactorAnalysis.randcor(K, 0.2)

d = CFADistributionSigma(Sigma_X, A, Sigma_L)
x = zeros(P + length(A.nzval) + div((K+1)*K,2))
FactorAnalysis.state2vec!(x, d)
x_copy = copy(x)
FactorAnalysis.vec2state!(d, x)
FactorAnalysis.state2vec!(x, d)
@test maximum(abs(x .- x_copy)) <= 1e-8

g = FactorAnalysis.CFADistributionSigmaGradient(Sigma_X, A, Sigma_L)
x = zeros(P + length(A.nzval) + div((K+1)*K,2))
FactorAnalysis.state2vec!(x, g)
x_copy = copy(x)
FactorAnalysis.vec2state!(g, x)
FactorAnalysis.state2vec!(x, g)
@test maximum(abs(x .- x_copy)) <= 1e-8

println("Testing CFADistributionSigma: loglikelihood...")
N = 100000
X = rand(d, N)
S = X*X' / N
Base.cov2cor!(S, sqrt(diag(S)))
@test loglikelihood(d, S, N) < 0.0
@test loglikelihood(d, S, N) > loglikelihood(CFADistributionSigma(Sigma_X, A, Sigma_L .+ eye(K)*10), S, N)
function ll_test(d, S, N)
    Theta = inv(d.Sigma_X + d.A*d.Sigma_L*d.A')
    logdet(Theta) - trace(S*Theta)
end
@test abs(loglikelihood(d, S, N) - ll_test(d, S, N)) < 1e-8

println("Testing CFADistributionSigma: rand...")
d = CFADistributionSigma(Sigma_X, A, Sigma_L)
@test size(rand(d, N)) == (P,N)

println("Testing CFADistributionSigma: gradient (compare to finite difference)...")
g = dloglikelihood(d, S, N)
eps = 1e-8
tol = 1e-4
for i in 1:P
    tmp = Sigma_X[i,i]
    Sigma_X[i,i] = tmp + eps
    l1 = loglikelihood(CFADistributionSigma(Sigma_X, A, Sigma_L), S, N)
    Sigma_X[i,i] = tmp - eps
    l2 = loglikelihood(CFADistributionSigma(Sigma_X, A, Sigma_L), S, N)
    Sigma_X[i,i] = tmp
    @test abs(g.Sigma_X[i,i] - (l1-l2)/(2*eps)) < tol
    #@test abs(dSigma_X[i,i] - (l1-l2)/(2*eps)) < tol
end
for i in 1:K, j in i+1:K
    tmp = Sigma_L[j,i]
    Sigma_L[j,i] = Sigma_L[i,j] = tmp + eps
    l1 = loglikelihood(CFADistributionSigma(Sigma_X, A, Sigma_L), S, N)
    Sigma_L[j,i] = Sigma_L[i,j] = tmp - eps
    l2 = loglikelihood(CFADistributionSigma(Sigma_X, A, Sigma_L), S, N)
    Sigma_L[j,i] = Sigma_L[i,j] = tmp
    @test abs(g.Sigma_L[j,i] - (l1-l2)/(2*eps)) < tol
end
rows = rowvals(A)
vals = nonzeros(A)
for col = 1:K
    for j in nzrange(A, col)
        tmp = vals[j]
        vals[j] = tmp + eps
        l1 = loglikelihood(CFADistributionSigma(Sigma_X, A, Sigma_L), S, N)
        vals[j] = tmp - eps
        l2 = loglikelihood(CFADistributionSigma(Sigma_X, A, Sigma_L), S, N)
        vals[j] = tmp
        @test abs(g.A[rows[j],col] - (l1-l2)/(2*eps)) < tol
        #@test abs(dA[rows[j],col] - (l1-l2)/(2*eps)) < tol
    end
end

println("Testing CFADistributionSigma: fit_mle...")
truthLL = loglikelihood(d, S, N)
dopt = fit_mle(CFADistributionSigma, ones(K,K), spones(A), S, N, show_trace=false, iterations=1000)
@test loglikelihood(dopt, S, N) > truthLL

# println("Testing CFADistributionSigma: fit_map...")
# dopt2 = fit_map(Normal(0, 1.0), CFADistributionSigma, spones(A), S, N, show_trace=true, iterations=1000)
# @test loglikelihood(dopt2, S, N) < loglikelihood(dopt, S, N)
# @test loglikelihood(dopt2, S, N) > truthLL

println("Testing CFADistributionSigma: fit_mle with a fully connected latent variable...")
srand(10)
P = 40
K = 20
Sigma_X = spdiagm(rand(P) .+ 0.2)
A = sparse(1:P, Int64[ceil(i/2) for i in 1:P], ones(P))
Sigma_L = FactorAnalysis.randcor(K, 0.2)
d = CFADistributionSigma(Sigma_X, A, Sigma_L)
N = 100000
X = rand(d, N)

# add a confounding latent factor
confounder = randn(P)
for i in 1:N
    X[:,i] .+= randn()*confounder
end

S = X*X' / N
Base.cov2cor!(S, sqrt(diag(S)))
truthLL = loglikelihood(d, S, N)

# create an independent component to handle the confounder
maskSigma_L = ones(K+1,K+1)
for i in 1:K
    maskSigma_L[K+1,i] = maskSigma_L[i,K+1] = 0
end
A2 = sparse([1:P; 1:P], [Int64[ceil(i/2) for i in 1:P]; ones(Int64, P)*(K+1)], ones(2*P))

dopt = fit_mle(CFADistributionSigma, maskSigma_L, A2, S, N, show_trace=false, iterations=15000)
Th = inv(dopt.Sigma_X + dopt.A*dopt.Sigma_L*dopt.A')
# display(dopt.A)
# tmp = deepcopy(dopt.Sigma_L)
# dopt.Sigma_L = dopt.Sigma_L * (2*eye(K+1))
# #Base.cov2cor!(dopt.Sigma_L, 2*ones(K+1))
# FactorAnalysis.normalize_Sigma_L!(dopt::CFADistributionSigma)
# display(dopt.A)
# display(tmp)
# println()
# display(dopt.Sigma_L)
# @test maximum(abs(tmp .- dopt.Sigma_L)) < 1e-8
#
# display(dopt.Sigma_X + dopt.A*dopt.Sigma_L*dopt.A')
# println()
# println("GGG2")
# display(S)
# println()
#
# println(logdet(Th) - trace(S*Th))
# println(logdet(inv(S)) - trace(S*inv(S)))
# # display(Sigma_X)
# display(dopt.Sigma_X)
# display(Sigma_L)
# println()
# display(dopt.Sigma_L)
# display(dopt.A)
@test loglikelihood(dopt, S, N) > truthLL
@test FactorAnalysis.area_under_pr(inv(Sigma_L), inv(dopt.Sigma_L[1:K,1:K])) > 0.9

# println(truthLL)
#println(v)
# println()
# v = FactorAnalysis.area_under_pr(abs(FactorAnalysis.upper(inv(Sigma_L))) .> 0.01, abs(FactorAnalysis.upper(inv(dopt.Sigma_L + 0.01*eye(K)))))
# println(v)

#display(Sigma_L)


text:
module BoundaryConditions

export
    BCType, Flux, Gradient, Value, NoPenetration,

    BoundaryCondition, bctype, getbc, setbc!,

    PeriodicBoundaryCondition, NoPenetrationBoundaryCondition, NoFluxBoundaryCondition,
    FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition,

    CoordinateBoundaryConditions,

    FieldBoundaryConditions, UVelocityBoundaryConditions, VVelocityBoundaryConditions,
    WVelocityBoundaryConditions, TracerBoundaryConditions, PressureBoundaryConditions,

    DiffusivityBoundaryConditions,

    BoundaryFunction, TracerBoundaryCondition,
    UVelocityBoundaryCondition, VVelocityBoundaryCondition, WVelocityBoundaryCondition,

    ParameterizedBoundaryCondition, ParameterizedBoundaryConditionFunction,

    apply_z_bcs!, apply_y_bcs!,

    fill_halo_regions!, zero_halo_regions!

using CUDAnative

using Oceananigans.Architectures
using Oceananigans.Grids

include("boundary_condition_types.jl")
include("boundary_condition.jl")
include("coordinate_boundary_conditions.jl")
include("field_boundary_conditions.jl")
include("boundary_function.jl")
include("parameterized_boundary_condition.jl")
include("show_boundary_conditions.jl")

include("fill_halo_regions.jl")
include("zero_halo_regions.jl")

include("apply_flux_bcs.jl")
include("apply_value_gradient_bcs.jl")
include("apply_no_penetration_bcs.jl")

end

ext: jl


text:
@testset "Remove reactions" begin
    lp = CoreModel(
        [1.0 1 1 0; 1 1 1 0; 1 1 1 0; 0 0 0 1],
        zeros(4),
        zeros(4),
        zeros(4),
        zeros(4),
        ["r1"; "r2"; "r3"; "r4"],
        ["m1"; "m2"; "m3"; "m4"],
    )
end

@testset "Find exchange reactions and metabolites" begin
    cp = test_LP()
    @test isempty(find_exchange_reactions(cp))
    @test isempty(find_exchange_metabolites(cp))

    cp = test_simpleLP()
    @test isempty(find_exchange_reactions(cp))
    @test isempty(find_exchange_metabolites(cp))

    cp = CoreModel(
        [-1.0 -1 -2; 0 -1 0; 0 0 0],
        zeros(3),
        ones(3),
        ones(3),
        ones(3),
        ["EX_m1"; "r2"; "r3"],
        ["m1"; "m2"; "m3"],
    )
    @test find_exchange_reactions(cp) == [1]

    cp = CoreModel(
        [-1.0 0 0; 0 0 -1; 0 -1 0],
        zeros(3),
        ones(3),
        ones(3),
        ones(3),
        ["EX_m1"; "Exch_m3"; "Ex_m2"],
        ["m1"; "m2"; "m3"],
    )
    @test find_exchange_reactions(cp) == [1; 2; 3]
    @test find_exchange_metabolites(cp) == [1; 3; 2]
    @test find_exchange_reactions(cp, exc_prefs = ["Exch_"]) == [2]
    @test find_exchange_metabolites(cp, exc_prefs = ["Exch_"]) == [3]

    cp = load_model(CoreModel, joinpath("data", "toyModel1.mat"))
    @test find_exchange_reactions(cp) == [4; 5; 6]
    @test find_exchange_metabolites(cp) == [4; 5; 6]
    @test find_exchange_reactions(cp, exclude_biomass = true) == [4; 5]
    @test find_exchange_metabolites(cp, exclude_biomass = true) == [4; 5]
    @test find_exchange_reactions(cp, exclude_biomass = true, biomass_str = "biom") ==
          [4; 5]
    @test find_exchange_metabolites(cp, exclude_biomass = true, biomass_str = "biom") ==
          [4; 5]
end

@testset "Change bounds" begin
    cp = test_LP()
    change_bounds!(cp, [3; 1], xl = [-10.0; -20], xu = [10.0; 20])
    @test cp.xl == [-20; 1; -10]
    @test cp.xu == [20; 1; 10]
    change_bounds!(
        cp,
        ["gibberish1"; "r3"; "r1"; "gibberish2"],
        xl = [0; -30.0; -40; 0],
        xu = [0; 30.0; 40; 0],
    )
    @test cp.xl == [-40; 1; -30]
    @test cp.xu == [40; 1; 30]
    change_bounds!(cp, ["r1"; "r3"], xl = [-50.0; -60])
    @test cp.xl == [-50; 1; -60]
end

@testset "Verify consistency" begin
    cp = test_LP()
    @test size(cp.S) == (4, 3)
    (new_reactions, new_mets) = verify_consistency(
        cp,
        reshape(cp.S[:, end], :, 1),
        [1.0, 2.0, 3.0, 4.0],
        [2.0],
        [-1.0],
        [1.0],
        ["r4"],
        ["m1", "m2", "m3", "m6"],
        [1],
        [4],
    )
    @test new_reactions == [1]
    @test new_mets == [4]

    (new_reactions, new_mets) = verify_consistency(
        cp,
        reshape(cp.S[:, end], :, 1),
        [1.0, 2.0, 3.0, 4.0],
        [2.0],
        [-1.0],
        [1.0],
        ["r1"],
        ["m1", "m2", "m3", "m6"],
        [],
        [4],
    )
    @test new_reactions == []
    @test new_mets == [4]
end

@testset "Add reactions (checking existence and consistency)" begin
    cp = test_LP()
    @test size(cp.S) == (4, 3)
    (new_cp, new_reactions, new_mets) = add_reactions(
        cp,
        cp.S[:, end],
        [1.0, 2.0, 3.0, 4.0],
        2.0,
        -1.0,
        1.0,
        check_consistency = true,
    )
    @test n_reactions(cp) + 1 == n_reactions(new_cp)

    (new_cp, new_reactions, new_mets) = add_reactions(
        cp,
        cp.S[:, end],
        [1.0, 2.0, 3.0, 4.0],
        2.0,
        -1.0,
        1.0,
        "r1",
        ["m1", "m2", "m3", "m6"],
        check_consistency = true,
    )
    @test n_reactions(cp) == n_reactions(new_cp)
    @test n_metabolites(cp) + 1 == n_metabolites(new_cp)
end

@testset "Add reactions" begin
    cp = test_LP()
    @test size(cp.S) == (4, 3)
    cp = add_reactions(cp, 2.0 * ones(4), 3 .* ones(4), 2.0, -1.0, 1.0)
    @test size(cp.S) == (8, 4)
    cp = add_reactions(cp, 2.0 * ones(4, 1), 3 .* ones(4), 2 .* ones(1), -ones(1), ones(1))
    @test size(cp.S) == (12, 5)
    cp = add_reactions(
        cp,
        2.0 * ones(4, 10),
        3 .* ones(4),
        2 .* ones(10),
        -ones(10),
        ones(10),
    )
    @test size(cp.S) == (16, 15)

    cp = test_sparseLP()
    @test size(cp.S) == (4000, 3000)
    cp = add_reactions(cp, 2.0 * sprand(4000, 0.5), 3 .* sprand(4000, 0.5), 2.0, -1.0, 1.0)
    @test size(cp.S) == (8000, 3001)
    cp = add_reactions(
        cp,
        2.0 * sprand(4000, 1, 0.5),
        3 .* sprand(4000, 0.5),
        2 .* sprand(1, 0.5),
        -sprand(1, 0.5),
        sprand(1, 0.5),
    )
    @test size(cp.S) == (12000, 3002)
    cp = add_reactions(
        cp,
        2.0 * sprand(4000, 1000, 0.5),
        3 .* sprand(4000, 0.5),
        2 .* sprand(1000, 0.5),
        -sprand(1000, 0.5),
        sprand(1000, 0.5),
    )
    @test size(cp.S) == (16000, 4002)

    cp = test_sparseLP()
    @test size(cp.S) == (4000, 3000)
    cp = add_reactions(cp, 2.0 * ones(4000), 3 .* ones(4000), 2.0, -1.0, 1.0)
    @test size(cp.S) == (8000, 3001)
    cp = add_reactions(
        cp,
        2.0 * ones(4000, 1),
        3 .* ones(4000),
        2 .* ones(1),
        -ones(1),
        ones(1),
    )
    @test size(cp.S) == (12000, 3002)
    cp = add_reactions(
        cp,
        2.0 * ones(4000, 1000),
        3 .* ones(4000),
        2 .* ones(1000),
        -ones(1000),
        ones(1000),
    )
    @test size(cp.S) == (16000, 4002)

    # proper subset of existing metabolites
    cp = test_LP()
    new_cp = add_reactions(cp, [-1.0], zeros(1), 1.0, 0.0, 1.0, "r4", ["m1"])
    @test n_reactions(cp) + 1 == n_reactions(new_cp)

    @test_throws DimensionMismatch add_reactions(
        cp,
        2.0 * ones(4000, 1),
        3 .* ones(4000),
        2 .* ones(2),
        -ones(1),
        ones(1),
    )
end

@testset "Remove reactions" begin
    cp = test_LP()
    @test size(cp.S) == (4, 3)
    cp = remove_reactions(cp, 2)
    @test size(cp.S) == (0, 2)
    cp = remove_reactions(cp, [2, 1])
    @test size(cp.S) == (0, 0)

    cp = test_LP()
    @test size(cp.S) == (4, 3)
    cp = remove_reactions(cp, "r0")
    @test size(cp.S) == (4, 3)
    cp = remove_reactions(cp, "r1")
    @test size(cp.S) == (0, 2)
    cp = remove_reactions(cp, ["r2"])
    @test size(cp.S) == (0, 1)
end


text:
using Shapefile: Shapefile
using GeoInterface: AbstractMultiPolygon, AbstractPoint
using WeightedOnlineStats: WeightedMean
using Dates: Day

function getlabeldict(shapepath,labelsym,T,labelsleft)
  t = Shapefile.Table(shapepath)
  labels = getproperty(t,labelsym)
  labeldict = Dict(T(i)=>stripc0x(labels[i]) for i in 1:length(labels) if T(i) in labelsleft)
  return Dict("labels"=>labeldict)
end
getlabeldict(shapepath, ::Nothing,T, labelsleft)=Dict{String,Any}()

function aggregate_out(allout, highmat, labelsleft,n)
  dsort = Dict(i[2]=>i[1] for i in enumerate(labelsleft))
  for jout in 1:size(allout,2)
    for iout in 1:size(allout,1)
      v = view(highmat,(iout*n-(n-1)):iout*n,(jout*n-(n-1)):jout*n,:)
      for value in skipmissing(v)
        allout[iout,jout,dsort[value]]+=1
      end
    end
  end
  map!(i->i/(n*n),allout,allout)
end

function cubefromshape_fraction(shapepath,lonaxis,lataxis;labelsym=nothing, T=Float64, nincrease=10)
  s = (length(lonaxis), length(lataxis))
  outmat = zeros(T,s.*nincrease)
  lon1,lon2 = get_bb(lonaxis)
  lat1,lat2 = get_bb(lataxis)
  rasterize!(outmat, shapepath, bb = (left = lon1, right=lon2, top=lat1, bottom=lat2))

  outmat = replace(outmat,zero(T)=>missing)
  labelsleft = collect(skipmissing(unique(outmat)))
  pp = getlabeldict(shapepath,labelsym,T,Set(labelsleft))["labels"]

  allout = zeros(T,s...,length(labelsleft))
  aggregate_out(allout,outmat,labelsleft,nincrease)

  if labelsym===nothing
    newax = CategoricalAxis("Label", labelsleft)
  else
    newax = CategoricalAxis(labelsym,[pp[l] for l in labelsleft])
  end

  return YAXArray(CubeAxis[lonaxis, lataxis, newax], allout)

end
function cubefromshape_single(shapepath, lonaxis, lataxis; labelsym = nothing, T=Int32)
  s = (length(lonaxis), length(lataxis))
  outmat = zeros(T,s)
  lon1,lon2 = get_bb(lonaxis)
  lat1,lat2 = get_bb(lataxis)
  rasterize!(outmat, shapepath, bb = (left = lon1, right=lon2, top=lat1, bottom=lat2))

  outmat = replace(outmat,zero(T)=>missing)
  labelsleft = collect(skipmissing(unique(outmat)))
  properties = getlabeldict(shapepath,labelsym,T,labelsleft)

  return YAXArray(CubeAxis[lonaxis, lataxis], outmat,properties)
end

"""
    cubefromshape(shapepath, refcube; samplefactor=nothing, labelsym = nothing, T=Int32)

Read the shapefile at `shapepath`, which is required to contain
a list of polygons to construct a cube. The polygons will be *rasterized* to the same grid and bounding box as the reference data cube `refcube`.
If the shapefile comes with additional labels in a .dbf file, the labels can be mapped to the respective field codes by
providing the label to choose as a symbol. If a `samplefactor` is supplied, the polygons will be rasterized on an n-times higher
resolution and the fraction of area within each polygon will be returned.
"""
cubefromshape(shapepath, c; kwargs...) = cubefromshape(shapepath, getAxis("Lon",c), getAxis("Lat",c);kwargs...)
function cubefromshape(args...; samplefactor=nothing, kwargs...)
  if samplefactor===nothing
    cubefromshape_single(args...; kwargs...)
  else
    cubefromshape_fraction(args...; nincrease = samplefactor, kwargs...)
  end
end

function prune_labels!(c::YAXArray)
  if haskey(c.properties,"labels")
    labelsleft = Set(skipmissing(unique(c.data)))
    dold = c.properties["labels"]
    dnew = Dict(k=>dold[k] for k in filter(i->in(i,labelsleft),keys(dold)))
    c.properties["labels"] = dnew
  end
  c
end
stripc0x(a) = replace(a, r"[^\x20-\x7e]"=> "")

function rasterize!(outar,shapepath;bb = (left = -180.0, right=180.0, top=90.0,bottom=-90.0),label=nothing)
  t = Shapefile.Table(shapepath)
  p = Shapefile.shapes(t)
  if length(p)>1
    rasterizepoly!(outar,p,bb)
  else
    rasterizepoly!(outar,p[1],bb)
  end
end

function rasterizepoly!(outmat,poly::Vector{<:Union{Missing,AbstractMultiPolygon}},bb)
    foreach(1:length(poly)) do ipoly
        rasterizepoly!(outmat,poly[ipoly],bb,value=ipoly)
    end
    outmat
end

function rasterizepoly!(outmat,poly::Vector{T},bb;value=one(eltype(outmat)), wrap = (left = -180.0,right = 180.0)) where T<:AbstractPoint
nx,ny = size(outmat)
resx = (bb.right-bb.left)/nx
resy = (bb.top-bb.bottom)/ny
xr = range(bb.left+resx/2,bb.right-resx/2,length=nx)
yr = range(bb.top-resy/2,bb.bottom+resy/2,length=ny)
wrapwidth = wrap.right-wrap.left

for (iy,pixelY) in enumerate(yr)
    nodeX = Float64[]
    j = length(poly)
    for i = 1:length(poly)
        p1 = poly[i]
        p2 = poly[j]
        if wrap !==nothing && abs(p1.x-p2.x)>wrapwidth
            p1,p2 = p1.x < p2.x ? (T(p1.x+wrapwidth,p1.y),T(p2.x,p2.y)) : (T(p1.x,p1.y),T(p2.x+wrapwidth,p2.y))
        end
        if (p1.y < pixelY) && (p2.y >= pixelY) || (p2.y < pixelY) && (p1.y >= pixelY)
            push!(nodeX, p1.x + (pixelY-p1.y)/(p2.y-p1.y)*(p2.x-p1.x))
        end
        j = i
    end
    #Add intersect points at start and end for wrapped polygons
    if wrap!==nothing && any(i->i>wrap.right,nodeX)
        push!(nodeX,wrap.left)
        push!(nodeX,wrap.right)
        map!(i->i>wrap.right ? i-wrapwidth : i,nodeX,nodeX)
    end
    sort!(nodeX)
    @assert(iseven(length(nodeX)))
    for i = 1:2:length(nodeX)
        outmat[searchsortedfirst(xr,nodeX[i]):searchsortedlast(xr,nodeX[i+1]),iy].=value
    end
end
    outmat
end

function rasterizepoly!(outmat,pp::Shapefile.Polygon,bb;value=one(eltype(outmat)))
    points = map(1:length(pp.parts)) do ipart
        i0 = pp.parts[ipart]+1
        iend = ipart == length(pp.parts) ? length(pp.points) : pp.parts[ipart+1]
        pp.points[i0:iend]
    end
    foreach(1:length(points)) do ipoly
        rasterizepoly!(outmat,points[ipoly],bb,value=value)
    end
    outmat
end

function getboundingbox(data)
  defaultval = (CartesianIndex(size(data)),CartesianIndex{ndims(data)}())
  r1, r2 = mapreduce(
    i->ismissing(data[i]) ? defaultval : (i,i),
    (x,y)->(min(x[1],y[1]),max(x[2],y[2])),
    CartesianIndices(data),
    init = defaultval
  )
  CartesianIndices(map(Colon(),r1.I,r2.I))
end

"""
    aggregate_shapefile(shapepath,cube,refdate; npast=3, nfuture=10)

Calculates spatially aggregated statistics over a shapefile specified by `shapepath`.
Returns a time series cube for variables in `cube` that contains `npast` time steps
prior to and `nfuture` time steps after the event.
"""
function aggregate_shapefile(shapepath,cube,refdate;npast=3,nfuture=10)

  craster = cubefromshape(shapepath, cube)
  bspace = getboundingbox(craster.data)
  blon,blat = map(i->CartesianIndices((i,)),bspace.indices)
  taxall = getAxis("Time",cube)
  itime = axVal2Index(taxall,refdate)
  refdate = taxall.values[itime]
  btime = LinearIndices((itime-npast:itime+nfuture,))
  csmallbox = cube[lon=blon,lat=blat,time=btime]
  maskcube = craster[lon=blon,lat=blat]
  if findAxis("Variable",cube)===nothing
    by = ()
    incax = ("lat","time")
  else
    by = (:var,)
    incax = ("lat","var","time")
  end
  ct = CubeTable(cdata = csmallbox, cmask = maskcube, include_axes=incax)

  by = (i->iszero(i.cmask),by...,i->Int((i.time-refdate)/Day(1)))
  fitres = cubefittable(ct,WeightedMean,:cdata,by=by,weight=i->cosd(i.lat),showprog=false)

  fr = fitres[Category1=false]
  renameaxis!(fr,"Category"=>RangeAxis("Time",(-npast:nfuture)*8))
  fr
end


text:
# psf-gauss.jl

export psf_gauss

"""
    psf_gauss( ; nx, nx_psf, fwhm_start, fwhm_end, fwhm, T)

Create depth-dependent Gaussian PSFs
having specified full-width half-maximum (FHWM) values.

# Options
- 'nx::Int = 128'
- 'nx_psf::Int = 11' (should be odd)
- 'fwhm_start::Real = 1'
- 'fwhm_end::Real = 4'
- 'fwhm::AbstractVector{<:Real} = LinRange(fwhm_start, fwhm_end, nx)'
- 'T::DataType == Float32'

Returned `psf` is `[nx_psf, nx_psf, nx]` where each PSF sums to 1.
"""
function psf_gauss( ;
    nx::Int = 128,
    nx_psf::Int = 11,
    fwhm_start::Real = 1,
    fwhm_end::Real = 4,
    fwhm::AbstractVector{<:Real} = LinRange(fwhm_start, fwhm_end, nx),
    T::DataType = Float32,
)
    isodd(nx_psf) || @warn("even nx_psf = $nx_psf ?")
    psf = zeros(T, nx_psf, nx_psf, nx)

    for iy in 1:nx # depth-dependent blur
        r = (-(nx_psf-1)÷2):((nx_psf-1)÷2)
        σ = fwhm[iy] / sqrt(log(256))
        r2 = abs2.(r / σ)
        tmp = @. exp(-π * (r2 + r2'))
        psf[:,:,iy] = tmp / sum(tmp)
    end
    return psf
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the API that converts Euler angles to rotations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/angle_to_rot.jl
# =======================================

# Functions: angle_to_rot
# -----------------------

@testset "Euler angles => Rotation (Float64)" begin
    T = Float64

    for rot_seq in valid_rot_seqs
        # Sample Euler angles.
        ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rot_seq)

        # Convert to DCM and quaternion.
        q1 = angle_to_rot(Quaternion, ea)
        D1 = angle_to_rot(DCM, ea)
        D2 = angle_to_rot(ea)
        @test eltype(q1) === eltype(D1) === eltype(D2) === T

        # Expected values.
        qe = angle_to_quat(ea)
        De = angle_to_dcm(ea)

        @test q1 ≈ qe
        @test D1 ≈ De
        @test D2 ≈ De
    end
end

@testset "Euler angles => Rotation (Float32)" begin
    T = Float32

    for rot_seq in valid_rot_seqs
        # Sample Euler angles.
        ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rot_seq)

        # Convert to DCM and quaternion.
        q1 = angle_to_rot(Quaternion, ea)
        D1 = angle_to_rot(DCM, ea)
        D2 = angle_to_rot(ea)
        @test eltype(q1) === eltype(D1) === eltype(D2) === T

        # Expected values.
        qe = angle_to_quat(ea)
        De = angle_to_dcm(ea)

        @test q1 ≈ qe
        @test D1 ≈ De
        @test D2 ≈ De
    end
end

# Functions: smallangle_to_rot
# ----------------------------

@testset "Small Euler angles => Rotation (Float64)" begin
    T = Float64

    # Sample three small angles.
    θx = _rand_ang(T) * T(0.001)
    θy = _rand_ang(T) * T(0.001)
    θz = _rand_ang(T) * T(0.001)

    # Convert to DCM and quaternion.
    q1 = smallangle_to_rot(Quaternion, θx, θy, θz)
    D1 = smallangle_to_rot(DCM, θx, θy, θz; normalize = false)
    D2 = smallangle_to_rot(DCM, θx, θy, θz; normalize = true)
    D3 = smallangle_to_rot(θx, θy, θz; normalize = false)
    D4 = smallangle_to_rot(θx, θy, θz; normalize = true)
    @test eltype(q1) === eltype(D1) === eltype(D2) === T

    # Expected result.
    qe = smallangle_to_quat(θx, θy, θz)
    D1e = smallangle_to_dcm(θx, θy, θz; normalize = false)
    D2e = smallangle_to_dcm(θx, θy, θz; normalize = true)

    @test qe ≈ q1
    @test D1 ≈ D1e
    @test D2 ≈ D2e
    @test D3 ≈ D1e
    @test D4 ≈ D2e
end

@testset "Small Euler angles => Rotation (Float32)" begin
    T = Float32

    # Sample three small angles.
    θx = _rand_ang(T) * T(0.001)
    θy = _rand_ang(T) * T(0.001)
    θz = _rand_ang(T) * T(0.001)

    # Convert to DCM and quaternion.
    q1 = smallangle_to_rot(Quaternion, θx, θy, θz)
    D1 = smallangle_to_rot(DCM, θx, θy, θz; normalize = false)
    D2 = smallangle_to_rot(DCM, θx, θy, θz; normalize = true)
    D3 = smallangle_to_rot(θx, θy, θz; normalize = false)
    D4 = smallangle_to_rot(θx, θy, θz; normalize = true)
    @test eltype(q1) === eltype(D1) === eltype(D2) === T

    # Expected result.
    qe = smallangle_to_quat(θx, θy, θz)
    D1e = smallangle_to_dcm(θx, θy, θz; normalize = false)
    D2e = smallangle_to_dcm(θx, θy, θz; normalize = true)

    @test qe ≈ q1
    @test D1 ≈ D1e
    @test D2 ≈ D2e
    @test D3 ≈ D1e
    @test D4 ≈ D2e
end


text:
"""
spin15_standalone.jl - t1 optimized pulses
"""
# Usually a good idea to ship code with a Project.toml file to specify the environment
import Pkg; Pkg.activate(@__DIR__)

using Altro
using Dierckx
using HDF5
using Interpolations
using LinearAlgebra
using Plots
using Printf
using RobotDynamics
using StaticArrays
using TrajectoryOptimization
const TO = TrajectoryOptimization

# paths
WDIR = get(ENV, "ROBUST_QOC_PATH", "../../")
EXPERIMENT_META = "spin"
EXPERIMENT_NAME = "spin15"
SAVE_PATH = joinpath(WDIR, "out", EXPERIMENT_META, EXPERIMENT_NAME)

# defs
# this type is used for specifying a target transformation
@enum GateType begin
    zpiby2 = 1
    ypiby2 = 2
    xpiby2 = 3
end

# this type specifies common save formats i use
@enum SaveType begin
    jl = 1
    samplejl = 2
    py = 3
end

# this type is used for choosing an Altro solver
@enum SolverType begin
    ilqr = 1
    alilqr = 2
    altro = 3
end

"""
get an unused file name like XXXXX_<save_file_path>/<save_file_name>.<extension>
where XXXXX is a unique numeric prefix
"""
function generate_save_file_path(extension, save_file_name, save_path)
    # Ensure the path exists.
    mkpath(save_path)

    # Create a save file name based on the one given; ensure it will
    # not conflict with others in the directory.
    max_numeric_prefix = -1
    for (_, _, files) in walkdir(save_path)
        for file_name in files
            if occursin("_$(save_file_name).$(extension)", file_name)
                max_numeric_prefix = max(parse(Int, split(file_name, "_")[1]))
            end
        end
    end

    save_file_name = "_$(save_file_name).$(extension)"
    save_file_name = @sprintf("%05d%s", max_numeric_prefix + 1, save_file_name)

    return joinpath(save_path, save_file_name)
end

"""
read_save - Read all data from an h5 file into memory.
"""
function read_save(save_file_path)
    dict = h5open(save_file_path, "r") do save_file
        dict = Dict()
        for key in names(save_file)
            dict[key] = read(save_file, key)
        end
        return dict
    end

    return dict
end

"""
sample_controls - Sample controls and d2controls_dt2
on the preferred time axis using a spline.
"""
function sample_controls(save_file_path; dt=DT_PREF, dt_inv=DT_PREF_INV,
                         plot=false, plot_file_path=nothing)
    # Grab data to sample from.
    save = read_save(save_file_path)
    controls = save["astates"][1:end - 1, (save["controls_idx"])]
    d2controls_dt2 = save["acontrols"][1:end, save["d2controls_dt2_idx"]]
    (control_knot_count, control_count) = size(controls)
    if "dt_idx" in keys(save)
        dts = save["acontrols"][1:end, save["dt_idx"]]
    elseif "dt" in keys(save)
        dts = save["dt"] * ones(control_knot_count)
    end
    time_axis = [0; cumsum(dts, dims=1)[1:end - 1]]

    # Construct time axis to sample over.
    final_time_sample = sum(dts)
    knot_count_sample = Int(floor(final_time_sample * dt_inv))
    # The last control should be DT_PREF before final_time_sample.
    time_axis_sample = Array(0:1:knot_count_sample - 1) * dt

    # Sample time_axis_sample via spline.
    controls_sample = zeros(knot_count_sample, control_count)
    d2controls_dt2_sample = zeros(knot_count_sample, control_count)
    for i = 1:control_count
        controls_spline = Spline1D(time_axis, controls[:, i])
        controls_sample[:, i] = map(controls_spline, time_axis_sample)
        d2controls_dt2_spline = Spline1D(time_axis, d2controls_dt2[:, i])
        d2controls_dt2_sample[:, i] = map(d2controls_dt2_spline, time_axis_sample)
    end

    # Plot.
    if plot
        DPI = 300
        MS_SMALL = 2
        ALPHA = 0.2
        fig = Plots.plot(dpi=DPI)
        Plots.scatter!(time_axis, controls[:, 1], label="controls data", markersize=MS_SMALL, alpha=ALPHA)
        Plots.scatter!(time_axis_sample, controls_sample[:, 1], label="controls fit",
                       markersize=MS_SMALL, alpha=ALPHA)
        Plots.scatter!(time_axis, d2controls_dt2[:, 1], label="d2_controls_dt2 data")
        Plots.scatter!(time_axis_sample, d2controls_dt2_sample[:, 1], label="d2_controls_dt2 fit")
        Plots.xlabel!("Time (ns)")
        Plots.ylabel!("Amplitude (GHz)")
        Plots.savefig(fig, plot_file_path)
        println("Plotted to $(plot_file_path)")
    end
    return (controls_sample, d2controls_dt2_sample, final_time_sample)
end

# system defs
FQ = 1.4e-2 #GHz
MAX_CONTROL_NORM_0 = 5e-1 #GHz
FBFQ_A = 0.202407
FBFQ_B = 0.5
NEGI = SA_F64[0   0  1  0 ;
              0   0  0  1 ;
              -1  0  0  0 ;
              0  -1  0  0 ;]
SIGMAX_ISO = SA_F64[0   1   0   0;
                    1   0   0   0;
                    0   0   0   1;
                    0   0   1   0]
SIGMAZ_ISO = SA_F64[1   0   0   0;
                    0  -1   0   0;
                    0   0   1   0;
                    0   0   0  -1]
NEGI_H0_ISO = pi * NEGI * SIGMAZ_ISO
NEGI_H1_ISO = pi * NEGI * SIGMAX_ISO
const FQ_NEGI_H0_ISO = FQ * NEGI_H0_ISO

# raw T1 times are in units of microseconds
T1_ARRAY = [
    1597.923, 1627.93, 301.86, 269.03, 476.33, 1783.19, 2131.76, 2634.50,
    4364.68, 2587.82, 1661.915, 1794.468, 2173.88, 1188.83,
    1576.493, 965.183, 560.251, 310.88
] * 1e3
FBFQ_ARRAY = [
    0.26, 0.28, 0.32, 0.34, 0.36, 0.38, 0.4,
    0.42, 0.44, 0.46, 0.465, 0.47, 0.475,
    0.48, 0.484, 0.488, 0.492, 0.5
]
const FBFQ_T1_SPLINE_ITP = extrapolate(interpolate((FBFQ_ARRAY,), T1_ARRAY, Gridded(Linear())), Flat())
@inline amp_fbfq_lo(amplitude) = -abs(amplitude) * FBFQ_A + FBFQ_B
@inline amp_t1_spline(amplitude) = FBFQ_T1_SPLINE_ITP(amp_fbfq_lo(amplitude))

# Define the optimization.
const CONTROL_COUNT = 1   # the sizes are one of the things that really NEED to be const
DT_PREF = 1e-2
DT_PREF_INV = 1e2
# static is the time step used for non time-optimal problems
DT_STATIC = DT_PREF
DT_STATIC_INV = DT_PREF_INV
# dt_init is the initial time step at all knot points for time-optimal problems
DT_INIT = DT_PREF
DT_INIT_INV = DT_PREF_INV
DT_MIN = DT_INIT / 2
DT_MAX = DT_INIT * 2
CONSTRAINT_TOLERANCE = 1e-8
AL_KICKOUT_TOLERANCE = 1e-7
PN_STEPS = 5

# Define the problem.
INITIAL_STATE1 = SA[1., 0, 0, 0]
INITIAL_STATE2 = SA[0., 1, 0, 0]
STATE_SIZE, = size(INITIAL_STATE1)
INITIAL_ASTATE = [
    INITIAL_STATE1; # state1, this is a quantum evolving according to the schroedinger equation
    INITIAL_STATE2; # state2, also a quantum state
    @SVector zeros(CONTROL_COUNT); # int_control, this is the integral of the control
    @SVector zeros(CONTROL_COUNT); # control, this is the control
    @SVector zeros(CONTROL_COUNT); # dcontrol_dt, this is the first time derivative of the control
    @SVector zeros(1); # int_gamma, this is a thing i want to minimize subject to the constraints of the state evolution
]
const ASTATE_SIZE, = size(INITIAL_ASTATE)
# definitions of target quantum states
ZPIBY2_1 = SA[1., 0, -1, 0] / sqrt(2)
ZPIBY2_2 = SA[0., 1, 0, 1] / sqrt(2)
YPIBY2_1 = SA[1., 1, 0, 0] / sqrt(2)
YPIBY2_2 = SA[-1., 1, 0, 0] / sqrt(2)
XPIBY2_1 = SA[1., 0, 0, -1] / sqrt(2)
XPIBY2_2 = SA[0., 1, -1, 0] / sqrt(2)

# state indices
const STATE1_IDX = 1:STATE_SIZE
const STATE2_IDX = STATE1_IDX[end] + 1:STATE1_IDX[end] + STATE_SIZE
const INT_CONTROLS_IDX = STATE2_IDX[end] + 1:STATE2_IDX[end] + CONTROL_COUNT
const CONTROLS_IDX = INT_CONTROLS_IDX[end] + 1:INT_CONTROLS_IDX[end] + CONTROL_COUNT
const DCONTROLS_DT_IDX = CONTROLS_IDX[end] + 1:CONTROLS_IDX[end] + CONTROL_COUNT
const INT_GAMMA_IDX = DCONTROLS_DT_IDX[end] + 1:DCONTROLS_DT_IDX[end] + 1
# control indices
const D2CONTROLS_DT2_IDX = 1:CONTROL_COUNT
const DT_IDX = D2CONTROLS_DT2_IDX[end] + 1:D2CONTROLS_DT2_IDX[end] + 1

# Specify logging.
VERBOSE = true
SAVE = true

# Use type parameter `TO` to specify whether or not the model is Time Optimal
struct Model{TO} <: AbstractModel
    # make sure `TO` is a boolean by defining the inner constructor
    Model(TO::Bool=true) = new{TO}()
end
RobotDynamics.state_dim(::Model) = ASTATE_SIZE
RobotDynamics.control_dim(::Model{true}) = CONTROL_COUNT + 1
RobotDynamics.control_dim(::Model{false}) = CONTROL_COUNT

function RobotDynamics.discrete_dynamics(
        ::Type{RK4},
        model::Model{true},
        x::StaticVector,
        u::StaticVector,
        t::Real,
        dt::Real
    )

    println("using dt from controls")
    u = control(z)
    dt = u[end]

    k1 = dynamics(model, x,        u, t       )*dt
	k2 = dynamics(model, x + k1/2, u, t + dt/2)*dt
	k3 = dynamics(model, x + k2/2, u, t + dt/2)*dt
    k4 = dynamics(model, x + k3,   u, t + dt  )*dt
	x + (k1 + 2k2 + 2k3 + k4)/6
end

# time-optimal dynamics
function RobotDynamics.dynamics(model::Model{true}, astate, acontrols, time)
    negi_h = (
        FQ_NEGI_H0_ISO
        + astate[CONTROLS_IDX][1] * NEGI_H1_ISO
    )
    delta_state1 = negi_h * astate[STATE1_IDX]
    delta_state2 = negi_h * astate[STATE2_IDX]
    delta_int_control = astate[CONTROLS_IDX]
    delta_control = astate[DCONTROLS_DT_IDX]
    delta_dcontrol_dt = acontrols[D2CONTROLS_DT2_IDX]
    # int_gamma is a bad thing that I want to be small. the
    # value of amp_t1_spline is proportional to the value of astate[CONTROLS_IDX][1]
    # so making int_gamma small is often at odds with keeping the controls small
    delta_int_gamma = amp_t1_spline(astate[CONTROLS_IDX][1])^(-1)
    return [
        delta_state1;
        delta_state2;
        delta_int_control;
        delta_control;
        delta_dcontrol_dt;
        delta_int_gamma;
        # I have been doing time optimal problems by multiplying the dynamics
        # by the dt in the acontrols vector, and telling TO that my
        # dt is 1.
    ] .* acontrols[DT_IDX][1]^2
end

# non time-optimal dynamics
#   NOTE: you could easily reduce code duplication by defining a single utility
#         dynamics function (e.g. `_dynamics`), and call it for both.
function RobotDynamics.dynamics(model::Model, astate, acontrols, time)
    negi_h = (
        FQ_NEGI_H0_ISO
        + astate[CONTROLS_IDX][1] * NEGI_H1_ISO
    )
    delta_state1 = negi_h * astate[STATE1_IDX]
    delta_state2 = negi_h * astate[STATE2_IDX]
    delta_int_control = astate[CONTROLS_IDX]
    delta_control = astate[DCONTROLS_DT_IDX]
    delta_dcontrol_dt = acontrols[D2CONTROLS_DT2_IDX]
    delta_int_gamma = amp_t1_spline(astate[CONTROLS_IDX][1])^(-1)
    return [
        delta_state1;
        delta_state2;
        delta_int_control;
        delta_control;
        delta_dcontrol_dt;
        delta_int_gamma;
    ]
end

# Running this function with the default arguments given here reproduces
# the problem that I want to solve.
function run_traj(;evolution_time=60., gate_type=xpiby2,
                  initial_save_file_path=nothing,
                  initial_save_type=jl, time_optimal=true,
                  solver_type=altro)
    # Choose dynamics
    # It is likely the first time you run the script than an error will be thrown
    # saying that the dynamics can not be found. You should be able to reinclude
    # the file and try running it again. There is some weird behavior with world
    # contexts that I don't understand yet due to the eval statements. There
    # might be another way to define the dynamics function conditionally
    # like this but I do not know one yet.

    # RESPONSE: You really don't want to call `eval` like this, it's almost always
    #           bad practice to have `eval` within runtime code. There's several
    #           ways to get around having different dynamics functions. One would be
    #           to simply have two different models. I've implemented a slighly more
    #           "Julian" version leveraging the power of multiple dispatch.

    # Convert to trajectory optimization language.
    n = ASTATE_SIZE
    t0 = 0.
    x0 = INITIAL_ASTATE
    if time_optimal
        # add one to make dt a decision variable
        m = CONTROL_COUNT + 1
    else
        m = CONTROL_COUNT
    end
    if gate_type == xpiby2
        target_state1 = XPIBY2_1
        target_state2 = XPIBY2_2
    elseif gate_type == ypiby2
        target_state1 = YPIBY2_1
        target_state2 = YPIBY2_2
    elseif gate_type == zpiby2
        target_state1 = ZPIBY2_1
        target_state2 = ZPIBY2_2
    end
    xf = [
        target_state1;
        target_state2;
        @SVector zeros(CONTROL_COUNT); # int_control
        @SVector zeros(CONTROL_COUNT); # control
        @SVector zeros(CONTROL_COUNT); # dcontrol_dt
        @SVector zeros(1); # int_gamma
    ]

    # Bound the control amplitude.
    x_max = SVector{n}([
        fill(Inf, STATE_SIZE);
        fill(Inf, STATE_SIZE);
        fill(Inf, CONTROL_COUNT);
        fill(MAX_CONTROL_NORM_0, CONTROL_COUNT); # control
        fill(Inf, CONTROL_COUNT);
        fill(Inf, 1)
    ])
    x_min = SVector{n}([
        fill(-Inf, STATE_SIZE);
        fill(-Inf, STATE_SIZE);
        fill(-Inf, CONTROL_COUNT);
        fill(-MAX_CONTROL_NORM_0, CONTROL_COUNT); # control
        fill(-Inf, CONTROL_COUNT);
        fill(-Inf, 1)
    ])
    # Controls start and end at 0.
    x_max_boundary = SVector{n}([
        fill(Inf, STATE_SIZE);
        fill(Inf, STATE_SIZE);
        fill(Inf, CONTROL_COUNT);
        fill(0, CONTROL_COUNT); # control
        fill(Inf, CONTROL_COUNT);
        fill(Inf, 1)
    ])
    x_min_boundary = SVector{n}([
        fill(-Inf, STATE_SIZE);
        fill(-Inf, STATE_SIZE);
        fill(-Inf, CONTROL_COUNT);
        fill(0, CONTROL_COUNT); # control
        fill(-Inf, CONTROL_COUNT);
        fill(-Inf, 1)
    ])
    # Bound dt.
    if time_optimal
        u_min = SVector{m}([
            fill(-Inf, CONTROL_COUNT);
            fill(sqrt(DT_MIN), 1); # dt
        ])
        u_max = SVector{m}([
            fill(Inf, CONTROL_COUNT);
            fill(sqrt(DT_MAX), 1); # dt
        ])
    else
        u_min = SVector{m}([
            fill(-Inf, CONTROL_COUNT);
        ])
        u_max = SVector{m}([
            fill(Inf, CONTROL_COUNT);
        ])
    end

    # Generate initial trajectory.
    model = Model(time_optimal)
    U0 = nothing
    if time_optimal
        # Default initial guess w/ optimization over dt.
        dt = 1
        N = Int(floor(evolution_time * DT_INIT_INV)) + 1
        U0 = [SVector{m}([
            fill(1e-4, CONTROL_COUNT);
            fill(DT_INIT, 1);
        ]) for k = 1:N - 1]
    else
        if initial_save_file_path == nothing
            # Default initial guess.
            dt = DT_STATIC
            N = Int(floor(evolution_time * DT_STATIC_INV)) + 1
            U0 = [SVector{m}(
                fill(1e-4, CONTROL_COUNT)
            ) for k = 1:N - 1]
        else
            # Initial guess pulled from initial_save_file_path.
            (d2controls_dt2, evolution_time) = h5open(initial_save_file_path, "r") do save_file
                 if initial_save_type == jl
                     d2controls_dt2_idx = read(save_file, "d2controls_dt2_idx")
                     d2controls_dt2 = read(save_file, "acontrols")[:, d2controls_dt2_idx]
                     evolution_time = read(save_file, "evolution_time")
                 elseif initial_save_type == samplejl
                     d2controls_dt2 = read(save_file, "d2controls_dt2_sample")
                     evolution_time = read(save_file, "evolution_time_sample")
                 end
                 return (d2controls_dt2, evolution_time)
            end
            # Without variable dts, evolution time will be a multiple of DT_STATIC.
            evolution_time = Int(floor(evolution_time * DT_STATIC_INV)) * DT_STATIC
            dt = DT_STATIC
            N = Int(floor(evolution_time * DT_STATIC_INV)) + 1
            U0 = [SVector{m}(d2controls_dt2[k, 1]) for k = 1:N-1]
        end
    end
    X0 = [SVector{n}(
        fill(NaN, n)
    ) for k = 1:N]
    Z = Traj(X0, U0, dt * ones(N))

    # Define penalties.
    Q = Diagonal(SVector{n}([
        # The state1 and state2 evolution is the key thing. I want to simulate their evolution
        # for very long periods of time using the controls obtained from this optimization.
        # Therefore, I need to make sure the optimized controls get very low tolerances
        # on reaching the target state, otherwise when I put the controls into the long
        # time simulation I get numerical error that builds with each iterate.
        fill(1e1, STATE_SIZE); # state 1
        fill(1e1, STATE_SIZE); # state 2
        fill(1e2, CONTROL_COUNT); # int
        fill(1e2, CONTROL_COUNT); # control
        fill(1e-1, CONTROL_COUNT); # dcontrol_dt
        # I would like this value to be as large as possible so int_gamma is as small as possible.
        # int_gamma is typically on the order of 1e-5 so this 1e7 is not that large.
        # I can typically afford to use 1e9 for similar problems.
        # COMMENT: This is screaming "numerical ill-conditioning."
        fill(1e7, 1); # int_gamma
    ]))
    Qf = Q * N
    if time_optimal
        R = Diagonal(SVector{m}([
            fill(1e-1, CONTROL_COUNT); # d2control_dt2
            fill(5e2, 1); # dt, I have found this typically needs to be high
        ]))
    else
        R = Diagonal(SVector{m}([
            fill(1e-1, CONTROL_COUNT); # d2control_dt2
        ]))
    end
    obj = LQRObjective(Q, R, Qf, xf, N)

    # Must satisfy control amplitude bound.
    control_bnd = BoundConstraint(n, m, x_max=x_max, x_min=x_min)
    # Must satisfy controls start and stop at zero.
    control_bnd_boundary = BoundConstraint(n, m, x_max=x_max_boundary, x_min=x_min_boundary)
    # Must satisfy dt bound.
    dt_bnd = BoundConstraint(n, m, u_max=u_max, u_min=u_min)
    # States must reach target. Controls must have zero net flux.
    target_astate_constraint = GoalConstraint(xf, [STATE1_IDX; STATE2_IDX; INT_CONTROLS_IDX])
    # States must obey unit norm. This constraint can be removed if necessary.
    # This constraint is rarely the reason for tolerance violation if dt is
    # adequately small / integration is accurate.
    normalization_constraint_1 = NormConstraint(n, m, 1, TO.Equality(), STATE1_IDX)
    normalization_constraint_2 = NormConstraint(n, m, 1, TO.Equality(), STATE2_IDX)

    constraints = ConstraintList(n, m, N)
    add_constraint!(constraints, control_bnd, 2:N-2)
    add_constraint!(constraints, control_bnd_boundary, N-1:N-1)
    add_constraint!(constraints, target_astate_constraint, N:N)
    add_constraint!(constraints, normalization_constraint_1, 2:N-1)
    add_constraint!(constraints, normalization_constraint_2, 2:N-1)
    if time_optimal
        add_constraint!(constraints, dt_bnd, 1:N-1)
    end

    # Instantiate problem and solve.
    prob = Problem(model, obj, xf, evolution_time, x0=x0,
        constraints=constraints, integration=RK4)
    initial_controls!(prob, U0)
    solver = nothing
    opts = SolverOptions(verbose=VERBOSE)
    if solver_type == alilqr
        solver = AugmentedLagrangianSolver(prob, opts)
        solver.opts.constraint_tolerance = CONSTRAINT_TOLERANCE
        solver.opts.constraint_tolerance_intermediate = CONSTRAINT_TOLERANCE
    elseif solver_type == altro
        solver = ALTROSolver(prob, opts)
        solver.opts.constraint_tolerance = CONSTRAINT_TOLERANCE
        # Let ALTRO do this for you via `projected_newton_tolerance`
        # solver.solver_al.opts.constraint_tolerance = AL_KICKOUT_TOLERANCE
        # solver.solver_al.opts.constraint_tolerance_intermediate = AL_KICKOUT_TOLERANCE
        solver.opts.projected_newton_tolerance = AL_KICKOUT_TOLERANCE
        solver.solver_pn.opts.constraint_tolerance = CONSTRAINT_TOLERANCE
        solver.solver_pn.opts.n_steps = PN_STEPS
        solver.solver_al.opts.iterations = 1
        solver.solver_al.solver_uncon.opts.iterations = 1
    end
    Altro.solve!(solver)

    # Post-process.
    acontrols_raw = controls(solver)
    acontrols_arr = permutedims(reduce(hcat, map(Array, acontrols_raw)), [2, 1])
    astates_raw = TO.states(solver)
    astates_arr = permutedims(reduce(hcat, map(Array, astates_raw)), [2, 1])
    Q_raw = Array(Q)
    Q_arr = [Q_raw[i, i] for i in 1:size(Q_raw)[1]]
    Qf_raw = Array(Qf)
    Qf_arr = [Qf_raw[i, i] for i in 1:size(Qf_raw)[1]]
    R_raw = Array(R)
    R_arr = [R_raw[i, i] for i in 1:size(R_raw)[1]]
    cidx_arr = Array(CONTROLS_IDX)
    d2cdt2idx_arr = Array(D2CONTROLS_DT2_IDX)
    dtidx_arr = Array(DT_IDX)
    # Square the dts.
    if time_optimal
        acontrols_arr[:, DT_IDX] = acontrols_arr[:, DT_IDX] .^2
    end
    cmax = TO.max_violation(solver)
    cmax_info = TO.findmax_violation(get_constraints(solver))

    # Save.
    if SAVE
        save_file_path = generate_save_file_path("h5", EXPERIMENT_NAME, SAVE_PATH)
        println("Saving this optimization to $(save_file_path)")
        h5open(save_file_path, "cw") do save_file
            write(save_file, "acontrols", acontrols_arr)
            write(save_file, "controls_idx", cidx_arr)
            write(save_file, "d2controls_dt2_idx", d2cdt2idx_arr)
            write(save_file, "dt_idx", dtidx_arr)
            write(save_file, "evolution_time", evolution_time)
            write(save_file, "astates", astates_arr)
            write(save_file, "Q", Q_arr)
            write(save_file, "Qf", Qf_arr)
            write(save_file, "R", R_arr)
            write(save_file, "solver_type", Integer(solver_type))
            write(save_file, "cmax", cmax)
            write(save_file, "cmax_info", cmax_info)
        end
        if time_optimal
            # Sample the important metrics.
            (controls_sample, d2controls_dt2_sample, evolution_time_sample) = sample_controls(save_file_path)
            h5open(save_file_path, "r+") do save_file
                write(save_file, "controls_sample", controls_sample)
                write(save_file, "d2controls_dt2_sample", d2controls_dt2_sample)
                write(save_file, "evolution_time_sample", evolution_time_sample)
            end
        end
    end
end

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
# Integration using a function
function trapz(f, a, b, h)
    x = collect(Float64, a:h:b);
    suma = 0
    for j ∈ 2:(length(x) - 1)
        val = f(x[j])
        suma += val
    end
    I = (b - a) * (f(x[1]) + 2 * suma + f(x[end])) / (2 * (length(x) - 1));
    return I
end

# Integration of a table
function trapz(x::Array{Float64,1}, y::Array{Float64,1})
    suma = 0;
    for j ∈ 2:(length(x) - 1)
        suma += y[j]
    end
    I = (x[end] - x[1]) * (y[1] + 2 * suma + y[end]) / (2 * (length(x) - 1));
    return I
end


text:
using TimedAutomata
using Test, LightGraphs, GraphPlot, Compose, IntervalArithmetic

# Clock related test

@testset "All Tests" begin
    op1 = PartialOrder(3, (1,2)=>:<, (2,3)=>:(<=))
    op2 = PartialOrder(3, (2,1)=>:(>=), (1,3)=>:<)
    op3 = PartialOrder(3, (2,3)=>:<)
    ot1 = TotalOrder([1,2,3])
    ot2 = TotalOrder([1,1,2])
    ot3 = TotalOrder([1,2,2])

    t1 = Clk"0"
    t2 = Clk"1"
    t3 = Clk"3"
    t4 = Clk"(0,1)"
    t5 = Clk"[0,2)"
    t6 = Clk"(0,∞)"
    t7 = Clk"(∞,1]"
    t8 = Clk"[0,∞)"
    t9 = Clk""

    c1 = ClockRegion{3,Int}(1=>Clk"(0,1)", 2=>Clk"[0,1]", 3=>Clk"1"; order=op1)
    c2 = ClockRegion{3,Int}(1=>Clk"(0,1)", 2=>Clk"[0,1]", 3=>Clk"1"; order=op2)
    c3 = ClockRegion{3,Int}(1=>Clk"(0,1)", 2=>Clk"[0,1]"; order=op2)
    c4 = ClockRegion{3,Int}(1=>Clk"(0,1)", 2=>Clk"[0,1)"; order=op2)
    c5 = ClockRegion{3,Int}(1=>Clk"(0,1)", 2=>Clk"1"; order=op2)

    @testset "Equality tests" begin
        @test Clk"1" == Clk"1"
        @test Clk"1" != Clk"0"
        @test Clk"(0,1]" == Clk"(0,1]"
        @test Clk"(0,1]" != Clk"(0,1)"

        @test ClockRegion{1,Int64}(TimeSlice{Int64}[TimePoint{Int64}(0)], TotalOrder{1}([0])) == ClockRegion{1,Int64}(TimeSlice{Int64}[TimePoint{Int64}(0)], TotalOrder{1}([0]))
    end

    @testset "Order Relation tests" begin
        expected_results = [
            true  true  false false false false;
            false true  false false false false;
            false false true  false false false;
            true  true  true  true  false false;
            false true  true  false true  false;
            true  true  false false false true;
        ]
        results = zeros(Bool, size(expected_results))

        for (i,r1) ∈ enumerate((op1,op2,op3,ot1,ot2,ot3))
            for (j,r2) ∈ enumerate((op1,op2,op3,ot1,ot2,ot3))
                results[i,j] = (r1 ⊆ r2)
                @test results[i,j] == expected_results[i,j]
            end
        end

        failures = findall(results .!= expected_results)

        if ~isempty(failures)
            println("Failed for: $(failures)")
        end
    end

    @testset "Time slice tests" begin
        expected_results = [
            true  false false false true  false true  true  true;
            false true  false false true  true  true  true  true;
            false false true  false false true  false true  true;
            false false false true  true  true  true  true  true;
            false false false false true  false false true  true;
            false false false false false true  false true  true;
            false false false false false false true  false true;
            false false false false false false false true  true;
            false false false false false false false false true;
        ]
        results = zeros(Bool, size(expected_results))

        for (i,r1) ∈ enumerate((t1,t2,t3,t4,t5,t6,t7,t8,t9))
            for (j,r2) ∈ enumerate((t1,t2,t3,t4,t5,t6,t7,t8,t9))
                results[i,j] = (r1 ⊆ r2)
                @test results[i,j] == expected_results[i,j]
            end
        end

        failures = findall(results .!= expected_results)

        if ~isempty(failures)
            println("Failed for: $(failures)")
        end
    end

    @testset "Clock Region tests" begin
        expected_results = [
            true  true  true false false;
            false true  true false false;
            false false true false false;
            false false true true  false;
            false false true false true
        ]
        results = zeros(Bool, size(expected_results))

        for (i,r1) ∈ enumerate((c1,c2,c3,c4,c5))
            for (j,r2) ∈ enumerate((c1,c2,c3,c4,c5))
                results[i,j] = (r1 ⊆ r2)
                @test results[i,j] == expected_results[i,j]
            end
        end

        failures = findall(results .!= expected_results)

        if ~isempty(failures)
            println("Failed for: $(failures)")
        end
    end

end

#=
tma = TMA(
    [
        [0,0],
        [0,1],
        [1,0],
        [1,1]
    ],
    Vector{Int}[],
    2,
    Dict(
        ([0,0],:a)=>[TMATransition([0,1],[1],@condition true)],
        ([0,1],:b)=>[TMATransition([1,0],[2],@condition true)],
        ([1,0],:c)=>[TMATransition([1,1],Int[],@condition c[1] < 1.0), TMATransition([0,0],Int[],@condition c[1] == 1.0)],
        ([1,1],:d)=>[TMATransition([0,0],Int[],@condition c[2] > 2.0)],
    ),
    Dict([1,0]=>@condition c[1] < 1.0),
    [0,0]
)
@test collect_constants(tma) == [Set((1//1,)),Set((2//1,))]
=#

ma = MA(
    [
        0,
        1,
        2
    ],
    [0],
    Dict(
        0=>[MATransition(1,:a),MATransition(2,:b)],
        1=>[MATransition(2,:a),MATransition(0,:b)],
        2=>[MATransition(0,:a),MATransition(1,:b)],
    ),
    0
)

L = language(ma)
println("Language of simple, timeless automaton:")
println(repr(L))

tma = TMA(
    [
        0,
        1,
        2,
        3
    ],
    [3],
    1,
    Dict(
        0=>[TMATransition(1,:msg,[1],ClockRegion{1,Int}())],
        1=>[TMATransition(1,:msg,[1],ClockRegion{1,Int}(1=>Clk"1")),TMATransition(1,:msg,[1],ClockRegion{1,Int}(1=>Clk"[1,∞)")),TMATransition(2,:msg,[1],ClockRegion{1,Int}(1=>Clk"(0,1)"))],
        2=>[TMATransition(3,:alarm,Int[],ClockRegion{1,Int}())]
    ),
    #Dict{Int,ClockRegion{1,Int}}(),
    Dict(2=>ClockRegion{1,Int}(1=>Clk"(0,1)")),
    0
)

@test collect_constants(tma) == [Set([0,1])]

ma = untime(tma)
draw(SVG("test_untimed.svg", 16cm, 16cm), gplot(ma))

L = language(ma)
println("Language of untimed automaton")
println(repr(L))

draw(SVG("test_timed.svg", 16cm, 16cm), gplot(tma))

#=
for s ∈ ma.Σ
    println("State $(s):")
    for ((s′,a),ee) ∈ pairs(ma.E)
        if s′==s
            for e ∈ ee
                println("\t$(a)=>$(e.s)")
            end
        end
    end
end
 =#


text:
# Julia wrapper for header: /usr/include/scip/scip_message.h
# Automatically generated using Clang.jl wrap_c

function SCIPsetMessagehdlr(scip, messagehdlr)
    ccall((:SCIPsetMessagehdlr, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_MESSAGEHDLR}), scip, messagehdlr)
end

function SCIPgetMessagehdlr(scip)
    ccall((:SCIPgetMessagehdlr, libscip), Ptr{SCIP_MESSAGEHDLR}, (Ptr{SCIP_},), scip)
end

function SCIPsetMessagehdlrLogfile(scip, filename)
    ccall((:SCIPsetMessagehdlrLogfile, libscip), Cvoid, (Ptr{SCIP_}, Cstring), scip, filename)
end

function SCIPsetMessagehdlrQuiet(scip, quiet)
    ccall((:SCIPsetMessagehdlrQuiet, libscip), Cvoid, (Ptr{SCIP_}, UInt32), scip, quiet)
end

function SCIPgetVerbLevel(scip)
    ccall((:SCIPgetVerbLevel, libscip), SCIP_VERBLEVEL, (Ptr{SCIP_},), scip)
end


text:
# small scrip to compute the solution of Lippmann Schwinger equation
# We test the types introduced in FastConvolutionSlowGO.jl
# we test that the application is fast and that the construction
# is performed fast.
# we test a different formulation of the Lippmann-Schwinger equation

using PyPlot
using IterativeSolvers
using SpecialFunctions

include("../src/FastConvolutionDual.jl")
include("../src/FastConvolution.jl")
include("../src/Preconditioner.jl")

FFTW.set_num_threads(4);
BLAS.set_num_threads(4)

#Defining Omega
h = 2.0^(-8)
k = 200.0/4

# size of box
a = 1;
b = 1;
x = collect(-a/2:h:a/2)
y = collect(-b/2:h:b/2)
(n,m) = length(x), length(y)
N = n*m
X = repmat(x, 1, m)[:];
Y = repmat(y', n,1)[:];
# we solve \triangle u + k^2(1 + nu(x))u = 0

# We use the modified quadrature in Duan and Rohklin
(ppw,D) = referenceValsTrapRule();
D0 = D[1];

# window(y,alpha) = 1*(abs.(y).<=alpha/2) + (abs.(y).>alpha/2).*(abs.(y).<alpha).*exp.(2*exp.(-0.5*alpha./(abs.(y)-alpha/2))./ ((abs.(y)-alpha/2)./(0.5*alpha)-1) )

# window(y,alpha, beta) = 1*(abs.(y).<=beta) + (abs.(y).>beta).*(abs.(y).<alpha).*exp.(2*exp.(-(alpha- beta)./(abs.(y)-beta))./ ((abs.(y)-beta)./(alpha- beta)-1) )

# # Defining the smooth perturbation of the slowness
# nu(x,y) = 0.05*(1+sin(2*pi*x/(0.96))).*window(y,0.96*b/2, 0.48*b/2).*window(x,0.96*0.5, 0.3);

# # Defining the smooth perturbation of the slowness
 nu(x,y) = -0.3*exp.(-40*(x.^2 + y.^2)).*(abs.(x).<0.48).*(abs.(y).<0.48);

figure(1); clf();
imshow(real(reshape( (1 + nu(X,Y)),n,m)), extent=[y[1], y[end], x[1], x[end]]);cb =  colorbar();

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

# comparison with the true solution
fastconvDual = buildFastConvolutionDual(x,y,h,k,nu,quadRule ="Greengard_Vico" );
#fastconvDual = buildFastConvolutionDual(x,y,h,k,nu);

u_inc = exp.(k*im*X);

rhsDual = -k^2*nu(X,Y).*u_inc ;

# allocating the solution
sigma = zeros(Complex128,N);

# solving the system using GMRES
@time info =  gmres!(sigma, fastconvDual, rhsDual, log = true, maxiter = 10 )
# println(info[2].residuals[:])

figure(4)
clf()
imshow(real(reshape(sigma,n,m))); colorbar();
title("density \sigma")

uDual = FFTconvolution(fastconvDual,sigma);

figure(5)
clf()
imshow(real(reshape(uDual,n,m))); colorbar();
title("scattered wavefield")

figure(6)
clf()
imshow(real(reshape(uDual,n,m))); colorbar();

figure(7)
clf()
imshow(real(reshape(uDual+u_inc,n,m))); colorbar();

# testing the preconditioner (Not working yet )

@time (As, ADGs) = buildGSparseAandG(k,X,Y,D0, n ,m, nu);

Mapproxsp = -As + k^2*ADGs;

precondDual = SparsifyingPreconditioner(Mapproxsp, As);

# allocating the solution
sigmaprecond = zeros(Complex128,N);

# solving the system using GMRES
@time infoPrecond =  gmres!(sigmaprecond, fastconvDual, rhsDual,
                            Pl= precondDual, log = true)
# println(infoPrecond[2].residuals[:])

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
module PromiseExtractor

using Gumbo
using URIParser

const EXTRACT_RULE = [(:script, x -> true, "src"),
                      (:img, x -> true, "src"),
                      (:link, x -> haskey(attrs(x), "rel") && getattr(x, "rel") == "stylesheet", "href")]

function prefixed(url, prefix)
    @assert startswith(prefix, "http")
    prefix_uri = URI(prefix)

    if startswith(url, "http")
        if startswith(url, prefix)
            return Nullable(url)
        else
            return Nullable{AbstractString}()
        end
    elseif startswith(url, "//")
        p = "//" + prefix.host + (prefix.port == 0x0 ? "" : dec(prefix.port)) + prefix.path

        if startswith(url, p)
            return Nullable(url)
        else
            return Nullable{AbstractString}()
        end
    else
        return Nullable(url)
    end
end

function extract_url(elem, prefix)
    for rule in EXTRACT_RULE
        t, f, a = rule[1], rule[2], rule[3]

        if t == tag(elem) && f(elem) && haskey(attrs(elem), a)
            url = getattr(elem, a)
            p = prefixed(url, prefix)

            if !isnull(p)
                return p
            end
        end
    end
    return Nullable{AbstractString}()
end

function promise_for(raw_html, prefix)
    doc = parsehtml(raw_html)

    results = []
    for elem in preorder(doc.root)
        url = extract_url(elem, prefix)

        if !isnull(url)
            push!(results, get(url))
        end
    end

    return results
end

export promise_for

end # module


text:
function start(integrator::ODEIntegrator)
  0
end

@inline function next(integrator::ODEIntegrator,state)
  state += 1
  step!(integrator) # Iter updated in the step! header
  # Next is callbacks -> iterator  -> top
  integrator,state
end

done(integrator::ODEIntegrator) = done(integrator,integrator.iter)

@inline function done(integrator::ODEIntegrator,state)
  if integrator.iter > integrator.opts.maxiters
    if integrator.opts.verbose
      warn("Interrupted. Larger maxiters is needed.")
    end
    postamble!(integrator)
    return true
  end
  if integrator.opts.unstable_check(integrator.dt,integrator.t,integrator.u)
    if integrator.opts.verbose
      warn("Instability detected. Aborting")
    end
    postamble!(integrator)
    return true
  end
  if !integrator.opts.force_dtmin && integrator.opts.adaptive && abs(integrator.dt) <= abs(integrator.opts.dtmin)
    if integrator.opts.verbose
      warn("dt <= dtmin. Aborting. If you would like to force continuation with dt=dtmin, set force_dtmin=true")
    end
    postamble!(integrator)
    return true
  end
  if isempty(integrator.opts.tstops)
    postamble!(integrator)
    return true
  elseif integrator.just_hit_tstop
    integrator.just_hit_tstop = false
    if integrator.opts.stop_at_next_tstop
      postamble!(integrator)
      return true
    end
  end
  false
end

@inline function step!(integrator::ODEIntegrator)
  if integrator.opts.advance_to_tstop
    @inbounds while integrator.tdir*integrator.t < integrator.tdir*top(integrator.opts.tstops)
      loopheader!(integrator)
      perform_step!(integrator,integrator.cache)
      loopfooter!(integrator)
    end
  else
    @inbounds loopheader!(integrator)
    @inbounds perform_step!(integrator,integrator.cache)
    @inbounds loopfooter!(integrator)
    @inbounds while !integrator.accept_step
      loopheader!(integrator)
      perform_step!(integrator,integrator.cache)
      loopfooter!(integrator)
    end
  end
  @inbounds handle_tstop!(integrator)
end

eltype(integrator::ODEIntegrator) = typeof(integrator)


text:
using JLD2, FileIO, Test

fn = joinpath(mktempdir(), "test.jld2")

# test iotype fallback
#  no fallback specified → throw method error
@test_throws MethodError JLD2.openfile(ArgumentError, fn, true, true, false, nothing)
#  fallback specified → switch to fallback
fh = JLD2.openfile(ArgumentError, fn, true, true, false, JLD2.MmapIO)
@test fh isa JLD2.MmapIO
# To avoid an mmap error on mac, have to write something to the stream before closing
JLD2.ensureroom(fh, 8)
write(fh, 42)
JLD2.truncate_and_close(fh, 8)

# test file path checking
Sys.isunix() && @test_throws ArgumentError jldopen("/dev/null", "r")
@test_throws ArgumentError jldopen(dirname(fn), "r")

# Test load macros
jldopen(fn, "w") do f
    write(f, "loadmacrotestvar1", ['a', 'b', 'c'])
    write(f, "loadmacrotestvar2", 1)
end

@eval begin # wrapped in eval since @load with no args needs file at compile time
function func1()
    @load $fn
    @test loadmacrotestvar1 == ['a', 'b', 'c']
    @test loadmacrotestvar2 == 1
end
end

func1()

function func2()
    @load fn loadmacrotestvar1 loadmacrotestvar2
    @test loadmacrotestvar1 == ['a', 'b', 'c']
    @test loadmacrotestvar2 == 1
end

func2()

@test !isdefined(@__MODULE__, :loadmacrotestvar1) # should not be in global scope
@test (@eval @load $fn) == [:loadmacrotestvar1, :loadmacrotestvar2]
@test loadmacrotestvar1 == ['a', 'b', 'c']
@test loadmacrotestvar2 == 1

# Test save macros
hello = "world"
@save fn hello
jldopen(fn, "r") do f
    @test read(f, "hello") == "world"
end

cmd = """
using JLD2
hello = "there"
@save $(repr(fn))
"""
run(`$(Base.julia_cmd()) -e $cmd`)
jldopen(fn, "r") do f
    @test read(f, "hello") == "there"
end

# Test Dict save/load
save(format"JLD2", fn, Dict("the"=>"quick", "brown"=>"fox", "stuff"=>reshape(1:4, (2, 2))))
@test load(fn) == Dict("the"=>"quick", "brown"=>"fox", "stuff"=>reshape(1:4, (2, 2)))

# Test load/save with pairs
save(format"JLD2", fn, "jumps", "over", "the", "lazy", "dog", reshape(1:4, (2, 2)))
@test load(fn, "jumps", "the", "dog") == ("over", "lazy", reshape(1:4, (2, 2)))
@test load(fn, "jumps") == "over"

jldopen(fn, "r+") do f
    @test !isempty(f)
    @test haskey(f, "dog")
    @test !haskey(f, "notdog")
    @test keys(f) == ["jumps", "the", "dog"]
    x1 = rand(UInt8, 1024^2+1024)
    x2 = rand(UInt8, 1024)
    f["x1"] = x1
    f["x2"] = x2
    @test f["x1"] == x1
    @test f["x2"] == x2
end

# Issue #19
save(fn, Dict("a"=>[1,2,3]))
io = open(fn)
@info("The next error message (involving \"loading nothing\") is a sign of normal operation")
@test_throws MethodError load(io)
close(io)

# Issue #33
d = Dict("params/p1" => 1,
         "params/p2" => 2.,
         "params/p3/p1" => 94,
         "data" => [[1,2,3], [4.,5.,6]])
save(fn, d)
@test load(fn) == d

# Issue #106
mutable struct MyMutableTest
    a::Int
    b::Vector{Int}
end
Base.getproperty(df::MyMutableTest, s::Symbol) =
    throw(ArgumentError("should not be called"))
Base.setproperty!(df::MyMutableTest, s::Symbol, x::Int) =
    throw(ArgumentError("should not be called"))
Base.isequal(x::MyMutableTest, y::MyMutableTest) =
    isequal(getfield(x, :a), getfield(y, :a)) && isequal(getfield(x, :b), getfield(y, :b))
mmtd = Dict("A" => MyMutableTest(1, [10]))
save(fn, mmtd)
@test isequal(load(fn), mmtd)

# Issue #125
len = 2^16
longstring = prod(fill("*",len));
lsd = Dict("longstring" => longstring)
save(fn, lsd)
@test isequal(load(fn), lsd)

# Testing Save macro
hello = "world"
num = 1.5

@save fn hello num
@test load(fn) == Dict("hello"=>"world", "num"=>1.5)

@save fn {compress=true} hello num
@test load(fn) == Dict("hello"=>"world", "num"=>1.5)

@save fn {compress=true, mmaparrays=false} hello num
@test load(fn) == Dict("hello"=>"world", "num"=>1.5)

@save fn {iotype=IOStream} hello num
@test load(fn) == Dict("hello"=>"world", "num"=>1.5)

@save fn bye = hello num
@test load(fn) == Dict("bye"=>"world", "num"=>1.5)

@save fn bye = hello num = 10
@test load(fn) == Dict("bye"=>"world", "num"=>10)

@test_throws ArgumentError @save fn {compress} hello

@test_throws ArgumentError @save fn hello=>"error"

# Issue # 189
struct RecursiveStruct
    x::RecursiveStruct
    RecursiveStruct() = new()
    RecursiveStruct(x) = new(x)
end

@testset "Recursive Immutable Types" begin
    x = RecursiveStruct()
    y = RecursiveStruct(x)

    @save "out.jld2" x y
    JLD2.jldopen("out.jld2", "r") do f
        @test f["x"] == x
        @test f["y"] == y
    end
end

# Issue #131
# write/read a Union{T,Missing}
len = 10_000
vect = Vector{Union{Bool,Missing}}(undef,len)
vect .= true
jldopen(fn,"w") do f
  f["vect"] = vect
end
vect_read = jldopen(fn,"r") do f
  f["vect"]
end
@test !any(ismissing.(vect_read))

# Also related to issue #131, but more types
len = 10_000
vect = Vector{Union{Missing,Float32,Float64,Int32}}(missing,len)
vect[vcat(1:10,33,44,55)] .= Int32(21)
vect[vcat(11:20,66,77,88)] .= 33.
vect[vcat(21:30,99)] .= Float32(123.)
jldopen(fn,"w") do f
  f["vect"] = vect
end
vect_read = jldopen(fn,"r") do f
  f["vect"]
end
@test all(findall(ismissing,vect) == findall(ismissing,vect_read))
@test all( skipmissing(vect) .=== skipmissing(vect_read))

# Issue #183
jfn, _ = mktemp()
@test_throws SystemError jldopen(jfn, "r", fallback = nothing)

# PR #206 Allow serialization of UnionAll in Union
struct UA1{T}; x::T; end
struct UA2{T}; y::T; end
@testset "UnionAll in Union" begin
    fn = joinpath(mktempdir(), "test.jld2")

    U1 = Union{Float64, Int}
    U2 = Union{Int, Vector}
    U3 = Union{UA1, UA2, Int}

    # Test types
    jldopen(fn, "w") do f
        f["u1"] = U1
        f["u2"] = U2
        f["u3"] = U3
    end

    u1, u2, u3 = jldopen(fn, "r") do f
        f["u1"], f["u2"], f["u3"]
    end
    @test u1 === U1
    @test u2 === U2
    @test u3 === U3
    # Test Vector with that eltype
    jldopen(fn, "w") do f
        f["u1"] = U1[1.0, 2, 3.0]
        f["u2"] = U2[1, [2.0], 3, ["4"]]
        f["u3"] = U3[UA1(1), UA2(2.0), 3, UA1("4")]
    end

    u1, u2, u3 = jldopen(fn, "r") do f
        f["u1"], f["u2"], f["u3"]
    end

    @test u1 == U1[1.0, 2, 3.0]
    @test u2 == U2[1, [2.0], 3, ["4"]]
    @test u3 == U3[UA1(1), UA2(2.0), 3, UA1("4")]
end

# Test for Issue #247
@testset "Tuple of Empty Structs" begin
    fn = joinpath(mktempdir(), "test.jld2")
    @save fn tup=(EmptyImmutable(), EmptyImmutable())
    @load fn tup

    @test tup == (EmptyImmutable(), EmptyImmutable())
end


text:
using Distributed
addprocs(10)
@everywhere using BioXMLSciMLTest, CellMLModelRepository, CellMLToolkit

# curl_cellml_models()
p = joinpath(@__DIR__, "C:/Users/Anand/.julia/dev/CellMLModelRepository/data/cellml_models")
fns = readdir(p;join=true)
# # @show fns
# @test !isempty(fns)

Pkg.add(url="https://github.com/anandijain/BioXMLSciMLTest.jl")
using BioXMLSciMLTest

f(fn) = getsys(CellModel(fn))

sciml(fns[1], f)

files_to_sciml(fns, f; pmap=true)
df = results_df("out/")
@test df isa DataFrame