text:
using DiffEqBiological
using DifferentialEquations
using Plots
using Colors

logistic_model = @reaction_network logistic begin
  r, cell + nutrient --> 2cell
end r

r = 1.5
cell0 = 1.0
nutrient0 = 2000.0

tspan = (0.0,0.008)
p =  (r)
u0 = [cell0, nutrient0]

oprob = ODEProblem(logistic_model, u0, tspan, p)
osol = solve(oprob, reltol=1e-8,abstol=1e-8)
Kbio = nutrient0+cell0
n = [((Kbio)*(cell0*exp(r*t)))/((Kbio)+(cell0*(exp(r*t)-1))) for t in osol.t]
#plot(osol)
plot(osol.t, [osol.u[i][1] for i in 1:length(osol.t)], labels="Cells", color= "blue", lw=2.0)
plot!(osol.t, [osol.u[i][2] for i in 1:length(osol.t)], labels="Nutrients", color="red", lw=2.0)

# https://nextjournal.com/sosiris-de/diffeqbiological
# What is the final value of cell(t)?  How does it relate to K?
# Can you overlay the analytical solution?

prob = DiscreteProblem(u0, tspan ,p)
jump_prob = JumpProblem(prob,Direct(),logistic_model)
sol = solve(jump_prob,FunctionMap())
#plot(sol)
plot(sol.t, [sol.u[i][1] for i in 1:length(sol.t)], labels="", color="blue", lw=2)
plot!(sol.t, [sol.u[i][2] for i in 1:length(sol.t)], labels="", color="red", lw=2)
# Can you plot the stochastic solution (sol) on top of the deterministic solution (osol)?

plot(osol)
plot!(sol)

#or

p1=plot(osol)
p2=plot(sol)

plot(p1,p2, xaxis="Time", yaxis="Population Size")

# Can you see how step-like the results are?
# What happens if you increase nutrient0 to 2000.0?
# Do the stochastic results still look very different to the deterministic ones?
# If you run line 27 above several times, saving the results in a list comprehension,
# then drawing them all on the same plot, will they be the same?

nsins=100
prob = DiscreteProblem(u0, tspan ,p)
jump_prob = JumpProblem(prob,Direct(),logistic_model)
solutions =[solve(jump_prob,FunctionMap()) for i in 1:nsins]
plot(solutions[1].t, [solutions[1].u[i][1] for i in 1:length(solutions[1].t)], labels="Cells", color= "blue", xaxis="Time", yaxis="Population Size", title="Discrete Stochastic Populations Dynamics", lw=0.3)
plot!(solutions[1].t,[solutions[1].u[i][2] for i in 1:length(solutions[1].t)], labels="Nutrients", color="red", lw=0.3)
for j in 2:nsins
    p=plot!(solutions[j].t, [solutions[j].u[i][1] for i in 1:length(solutions[j].t)], labels="", color= "blue", lw=0.3)
    q=plot!(solutions[j].t,[solutions[j].u[i][2] for i in 1:length(solutions[j].t)], labels="", color="red", lw=0.3)
    display(p)
    display(q)
end
oprob = ODEProblem(logistic_model, u0, tspan, p)
osol = solve(oprob, reltol=1e-8,abstol=1e-8)
Kbio = nutrient0+cell0
n = [((Kbio)*(cell0*exp(r*t)))/((Kbio)+(cell0*(exp(r*t)-1))) for t in osol.t]
plot!(osol.t, [osol.u[i][1] for i in 1:length(osol.t)], labels="Deterministic Continuous Solution", color="black", lw=2)
plot!(osol.t, [osol.u[i][2] for i in 1:length(osol.t)], labels="", color="black", legend=:right, lw=2)

# Can you simulate from the discrete stochastic exponential model?

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

using DiffEqBiological
using DifferentialEquations
using Plots
using Colors

exponential_model = @reaction_network exponential begin
  r, cell  --> 2cell
end r

r = 5.0
cell0 = 10.0

tspan = (0.0,1.0)
p = (r)
u0 = [cell0]

oprobe = ODEProblem(exponential_model, u0, tspan, p)
osole = solve(oprobe)
plot(osole, labels="Deterministic Solution", lw=2)

probe = DiscreteProblem(u0, tspan ,p)
jump_prob = JumpProblem(probe,Direct(),exponential_model)
sole = solve(jump_prob,FunctionMap())
#plot(sole)
plot!(sole, title="Deterministic and Stochastic Exponential Model", xaxis="Time", yaxis="Popuplation Size", legend=:left, lw=2, labels="Stochastic Solution")

nsims=100
probe = DiscreteProblem(u0, tspan ,p)
jump_prob = JumpProblem(probe,Direct(),exponential_model)
solutions =[solve(jump_prob,FunctionMap()) for i in 1:nsims]
plot(solutions[1].t, [solutions[1].u[i][1] for i in 1:length(solutions[1].t)], labels="Cells", color="blue", lw=0.2)
for j in 2:nsims
    p=plot!(solutions[j].t, [solutions[j].u[i][1] for i in 1:length(solutions[j].t)], labels="", color="blue", legend=:left, title="Exponential Model w/ Stochastic and Deterministic Solutions", lw=0.2)
    display(p)
end
oprobe = ODEProblem(exponential_model, u0, tspan, p)
osole = solve(oprobe)
plot!(osole, labels="Deterministic Solution", xaxis="Time", yaxis="Population Size",color="black", lw=2)

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

# Can you simulate from the discrete stochastic Lotka-Volterra model?

lv_model = @reaction_network lv begin
  k1, prey --> 2prey
  k2, prey + predator --> 2predator
  k3, predator --> null
 # k4, prey --> null
end k1 k2 k3 #k4

k1 = 1.0
k2 = 0.1
k3 = 0.1
#k4 = 0.001

prey0 = 40.0
predator0 = 100.0
null0 = 0.0
tspan = (0.0,200.0)
p =  (k1,k2,k3)
u0 = [prey0, predator0]

oproblv = ODEProblem(lv_model, u0, tspan, p)
osollv = solve(oproblv, reltol=1e-5,abstol=1e-5)
#plot(osollv)
plot(osollv.t, [osollv.u[i][1] for i in 1:length(osollv.t)], labels="D-Prey", color= "blue", lw=2)
plot!(osollv.t, [osollv.u[i][2] for i in 1:length(osollv.t)], labels="D-Predators", color="orange", lw=2)

prob = DiscreteProblem(u0, tspan ,p)
jump_prob = JumpProblem(prob,Direct(),lv_model)
sollv = solve(jump_prob,FunctionMap())
#plot(sollv)
plot!(sollv.t, [sollv.u[i][1] for i in 1:length(sollv.t)], labels="S-Prey", color= "green", lw=2.0)
plot!(sollv.t, [sollv.u[i][2] for i in 1:length(sollv.t)], labels="S-Predators", color="red", lw=2.0, title="Lotka-Volterra Deterministic and Stochastic Models", xaxis="Time", yaxis="Population Size")

# at larger pops the stochastic solution looks very similar to the deterministic solution but it doesnt have the advantage of being continuous

# Can you simulate from the discrete stochastic birth-death model?

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

bd_model = @reaction_network bd begin
  lambda, cell  --> 2cell
  mu, cell  --> null
end lambda mu

lambda = 2.0
#lambda = 0.8
mu= 1.5
#mu= 0.5
cell0 = 100.0
null0 = 0.0

tspan = (0.0,5.0)
p =  (lambda, mu)
u0 = [cell0, null0]

oprobbd = ODEProblem(bd_model, u0, tspan, p)
osolbd = solve(oprobbd, reltol=1e-8,abstol=1e-8)
#plot(osolbd)
plot(osolbd.t, [osolbd.u[i][1] for i in 1:length(osolbd.t)], labels="", color= "blue", ls=:dash)
plot!(osolbd.t, [osolbd.u[i][2] for i in 1:length(osolbd.t)], labels="", color="red", ls=:dash)

probbd = DiscreteProblem(u0, tspan ,p)
jump_prob = JumpProblem(probbd,Direct(),bd_model)
solbd = solve(jump_prob,FunctionMap())
#plot(solbd)
plot!(solbd.t, [solbd.u[i][1] for i in 1:length(solbd.t)], labels="Cells", color= "blue")
plot!(solbd.t, [solbd.u[i][2] for i in 1:length(solbd.t)], labels="Dead Cells", color="red", legend=:left, title="Birth-Death Model w/ Deterministic and Stochastic Solutions", xaxis="Time", yaxis="Popuplation Size")

nsins=100
probbd = DiscreteProblem(u0, tspan ,p)
jump_prob = JumpProblem(probbd,Direct(),bd_model)
solutions =[solve(jump_prob,FunctionMap()) for i in 1:nsins]
plot(solutions[1].t, [solutions[1].u[i][1] for i in 1:length(solutions[1].t)], labels="Cells", color= "blue", xaxis="Time", yaxis="Population Size", title="Discrete Stochastic Populations Dynamics", lw=0.3)
plot!(solutions[1].t,[solutions[1].u[i][2] for i in 1:length(solutions[1].t)], labels="Dead Cells", color="red", lw=0.2)
for j in 2:nsins
    p=plot!(solutions[j].t, [solutions[j].u[i][1] for i in 1:length(solutions[j].t)], labels="", color= "blue", lw=0.2)
    q=plot!(solutions[j].t,[solutions[j].u[i][2] for i in 1:length(solutions[j].t)], labels="", color="red", lw=0.2)
    display(p)
    display(q)
end
oprobbd = ODEProblem(bd_model, u0, tspan, p)
osolbd = solve(oprobbd, reltol=1e-10,abstol=1e-10)
plot!(osolbd.t, [osolbd.u[i][1] for i in 1:length(osolbd.t)], labels="Deterministic Solution", color="black", lw=3)
plot!(osolbd.t, [osolbd.u[i][2] for i in 1:length(osolbd.t)], labels="", color="black", legend=:left, lw=3, title="Birth-Death Model Showing Stochastic Behavior")

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

# What is the difference between the birth-death model and the exponential model?

using DiffEqBiological
using DifferentialEquations
using Plots
using Colors
#using ParameterizedFunctions

mtdna_model = @reaction_network mtdna begin
  b, wildtype --> 2wildtype
  d, wildtype  --> null
  bmut, mut --> 2mut
  dmut, mut --> null
  m, wildtype --> wildtype + mut
  # can we have wildtype --> 2mut
end b d bmut dmut m

# parameter values

b = 0.001875*365
d = 0.001875*365
bmut = 0.001875*365
dmut = 0.001875*365
m = (1e-5)*365

# initial conditions

wildtype0 = 100.0
mut0 = 0.0
null0 = 0.0
dummy0 = 1.0
tspan= (0.0, 100.0)
p = (b,d,bmut,dmut,m)
u0 = [wildtype0, mut0, null0]

sim_to_tot = function(res, totmax, b_small, b_big, replicative_advantage)
  u0 = res.u[end]
  tot_mtDNA0 = u0[1] + u0[3]
  pnew = collect(p)
  if tot_mtDNA0 >= totmax
    pnew[1] = b_small
    pnew[3] = b_small * replicative_advantage
  else
    pnew[1] = b_big
    pnew[3] = b_big * replicative_advantage
  end
  pnewer = (pnew...,)
  oprobmtdna = ODEProblem(mtdna_model, u0, tspan, pnewer)#, callback=cset)
  osolmtdna = solve(oprobmtdna)
  totmtDNA = [osolmtdna.u[i][1] + osolmtdna.u[i][3] for i in 1:length(osolmtdna.t)];
  triggered = [t>=totmax for t in totmtDNA];
  lastres = findfirst(triggered)-1;
  times = [t for t in osolmtdna.t[1:lastres]]
  species = [u for u in osolmtdna.u[1:lastres]]
  res = (t = times, u = species)
  return(res)
end

res = (t = [0.0], u = [u0])
current_time = res.t[1]
allres = []
while current_time <= tspan[2]
  global current_time, allres, b
  resnew = sim_to_tot(res, 105.0, 0.0, b, 1.0)
  ttotal = [t + current_time for t in resnew.t]
  current_time = ttotal[end]
  updatedres = (t=ttotal, u=resnew.u)
  push!(allres,updatedres)
end

times = collect(Iterators.flatten([res.t for res in allres]))
species = collect(Iterators.flatten([res.u for res in allres]))
osolmtdna = (t=times,u=species)

plot(osolmtdna.t, [osolmtdna.u[i][1] for i in 1:length(osolmtdna.t)], labels="Wildtype", color= "blue")
plot!(osolmtdna.t, [osolmtdna.u[i][3] for i in 1:length(osolmtdna.t)], labels="Mutant", color="red")
plot!(osolmtdna.t, [osolmtdna.u[i][1]+osolmtdna.u[i][3] for i in 1:length(osolmtdna.t)], labels="Total", color="black")

using DiffEqBiological
using DifferentialEquations
using Plots
using Colors

mtdna_model = @reaction_network mtdna begin
  b, wildtype --> 2wildtype
  d, wildtype  --> null
  bmut, mut --> 2mut
  dmut, mut --> null
  m, wildtype --> wildtype + mut
  # can we have wildtype --> 2mut
end b d bmut dmut m

# parameter values

b = 0.001875*365
d = 0.001875*365
bmut = 0.001875*365
dmut = 0.001875*365
m = (1e-5)*365

# initial conditions

wildtype0 = 100.0
mut0 = 10.0
null0 = 0.0
dummy0 = 1.0
tspan= (0.0, 100.0)
p = (b,d,bmut,dmut,m)
u0 = [wildtype0, mut0, null0]

sim_to_tot2 = function(res, totmax, b_small, b_big, replicative_advantage)
  u0 = res.u[end]
  tot_mtDNA0 = u0[1] + u0[3]
  pnew = collect(p)
  if tot_mtDNA0 >= totmax
    pnew[1] = b_small
    pnew[3] = b_small * replicative_advantage
  else
    pnew[1] = b_big
    pnew[3] = b_big * replicative_advantage
  end
  pnewer = (pnew...,)
  probmtdna = DiscreteProblem(u0, tspan ,pnewer)
  jump_prob = JumpProblem(probmtdna,Direct(),mtdna_model)
  solmtdna = solve(jump_prob,FunctionMap())
  totmtDNA = [solmtdna.u[i][1] + solmtdna.u[i][3] for i in 1:length(solmtdna.t)];
  triggered = [t>=totmax for t in totmtDNA];
  lastres = findfirst(triggered)-1;
  times = [t for t in solmtdna.t[1:lastres]]
  species = [u for u in solmtdna.u[1:lastres]]
  res = (t = times, u = species)
  return(res)
end

res = (t = [0.0], u = [u0])
current_time = res.t[1]
allres = []
while current_time <= tspan[2]
  global current_time, allres, b
  resnew = sim_to_tot2(res, 105.0, 0.0, b, 1.0)
  ttotal = [t + current_time for t in resnew.t]
  current_time = ttotal[end]
  updatedres = (t=ttotal, u=resnew.u)
  push!(allres,updatedres)
end

times = collect(Iterators.flatten([res.t for res in allres]))
species = collect(Iterators.flatten([res.u for res in allres]))
solmtdna = (t=times,u=species)

plot(solmtdna.t, [solmtdna.u[i][1] for i in 1:length(solmtdna.t)], labels="Wildtype", color= "blue")
plot!(solmtdna.t, [solmtdna.u[i][3] for i in 1:length(solmtdna.t)], labels="Mutant", color="red")
plot!(solmtdna.t, [solmtdna.u[i][1] + solmtdna.u[i][3] for i in 1:length(solmtdna.t)], labels="Total", color="black", legend=:outerright, title="mtDNA Population Dynamics Within a Cell", xaxis="Time (Years)", yaxis="Population Size")

nsins=10
probmtdna = DiscreteProblem(u0, tspan ,p)
jump_prob = JumpProblem(probmtdna,Direct(),mtdna_model)
solutions =[solve(jump_prob,FunctionMap()) for i in 1:nsins]
plot(solutions[1].t, [solutions[1].u[i][1] for i in 1:length(solutions[1].t)], labels="Wildtype", color= "blue", xaxis="Time", yaxis="Population Size", title="mtDNA Populations Dynamics", lw=0.3)
plot!(solutions[1].t,[solutions[1].u[i][2] for i in 1:length(solutions[2].t)], labels="Mutant", color="red", lw=0.3)
for j in 2:nsins
    p=plot!(solutions[j].t, [solutions[j].u[i][1] for i in 1:length(solutions[j].t)], labels="", color= "blue", lw=0.3)
    q=plot!(solutions[j].t,[solutions[j].u[i][3] for i in 1:length(solutions[j].t)], labels="", color="red", lw=0.3)
    display(p)
    display(q)
end
oprobmtdna = ODEProblem(mtdna_model, u0, tspan, p)
osolmtdna = solve(oprobmtdna)
plot!(osolmtdna.t, [osolmtdna.u[i][1] for i in 1:length(osolmtdna.t)], labels="Deterministic Continuous Solution", color="black", lw=2)
plot!(osolmtdna.t, [osolmtdna.u[i][3] for i in 1:length(osolmtdna.t)], labels="", color="black", legend=:right, lw=2)

# What happens if a cells is born with no mutations?

# What happens if a cell inherits a lot of mutations?  70%?

# What happens in both of the scenarios above if bmut > b?


text:
@testset "Energy expressions" begin
    ome(a, op, b) = EnergyExpressions.OrbitalMatrixElement(a, op, b)
    ome(a::AbstractOrbital, op, b::AbstractOrbital) = ome([a], op, [b])
    c(a, op, b) = ContractedOperator(Tuple(e for e in a), op, Tuple(e for e in b))
    c(a::AbstractOrbital, op, b::AbstractOrbital) = c([a], op, [b])
    eq(orb, op, args...) = NBodyEquation(orb, op, args...)
    ceq(orb, op, args...) = eq(Conjugate(orb), op, args...)
    ov(a,b) = OrbitalOverlap(a,b)
    nbt(args...) = EnergyExpressions.NBodyTerm([args...], 1)

    @testset "Average-of-configuration" begin
        for gst = (c"1s²", rc"1s²")
            cfgs = [gst]

            os = sort(unique(reduce(vcat, orbitals.(cfgs))))
            o = os[1]
            v = csfs(cfgs)

            h = FieldFreeOneBodyHamiltonian()
            g = CoulombInteraction()
            H = h + g

            E = Matrix(H, cfgs)
            eqs = diff(E, conj.(os))

            @test length(eqs.equations) == 1
            eq1 = eqs.equations[1]
            @test eq1.orbital == o
            @test size(eq1.equation) == (1,1)
            eq1 = eq1.equation[1]

            @test eq1 == 2eq(o, h) + 2eq(o, c(o, CoulombInteractionMultipole(0, g), o))
        end
    end
end

ext: jl


text:
"""
    AbstractFlowAlgorithm

Abstract type that allows users to pass in their preferred algorithm
"""
abstract type AbstractFlowAlgorithm end

"""
    EdmondsKarpAlgorithm <: AbstractFlowAlgorithm

Forces the maximum_flow function to use the Edmonds–Karp algorithm.
"""
struct EdmondsKarpAlgorithm <: AbstractFlowAlgorithm end

"""
    DinicAlgorithm <: AbstractFlowAlgorithm

Forces the maximum_flow function to use Dinic's algorithm.
"""
struct DinicAlgorithm <: AbstractFlowAlgorithm end

"""
    BoykovKolmogorovAlgorithm <: AbstractFlowAlgorithm

Forces the maximum_flow function to use the Boykov-Kolmogorov algorithm.
"""
struct BoykovKolmogorovAlgorithm <: AbstractFlowAlgorithm end

"""
Forces the maximum_flow function to use the Push-Relabel algorithm.
"""
struct PushRelabelAlgorithm <: AbstractFlowAlgorithm end

"""
    DefaultCapacity{T}

Structure that returns `1` if a forward edge exists in `flow_graph`, and `0` otherwise.
"""
struct DefaultCapacity{T<:Integer} <: AbstractMatrix{T}
    flow_graph::DiGraph
    nv::T
end

@traitfn DefaultCapacity(flow_graph::::IsDirected) =
    DefaultCapacity(DiGraph(flow_graph), nv(flow_graph))

getindex(d::DefaultCapacity{T}, s::Integer, t::Integer) where T = if has_edge(d.flow_graph, s, t) one(T) else zero(T) end
# isassigned{T<:Integer}(d::DefaultCapacity{T}, u::T, v::T) = (u in 1:d.nv) && (v in 1:d.nv)
size(d::DefaultCapacity) = (Int(d.nv), Int(d.nv))
transpose(d::DefaultCapacity) = DefaultCapacity(reverse(d.flow_graph))
ctranspose(d::DefaultCapacity) = DefaultCapacity(reverse(d.flow_graph))

"""
    residual(flow_graph)

Return a directed residual graph for a directed `flow_graph`.

The residual graph comprises the same node list as the orginal flow graph, but
ensures that for each edge (u,v), (v,u) also exists in the graph. This allows
flow in the reverse direction.

If only the forward edge exists, a reverse edge is created with capacity 0.
If both forward and reverse edges exist, their capacities are left unchanged.
Since the capacities in [`DefaultDistance`](@ref) cannot be changed, an array of ones
is created.
"""
function residual end
@traitfn residual(flow_graph::::IsDirected) = DiGraph(Graph(flow_graph))

# Method for Edmonds–Karp algorithm

@traitfn function maximum_flow(
    flow_graph::::IsDirected,                   # the input graph
    source::Integer,                       # the source vertex
    target::Integer,                       # the target vertex
    capacity_matrix::AbstractMatrix,   # edge flow capacities
    algorithm::EdmondsKarpAlgorithm        # keyword argument for algorithm
    )
    residual_graph = residual(flow_graph)
    return edmonds_karp_impl(residual_graph, source, target, capacity_matrix)
end

# Method for Dinic's algorithm

@traitfn function maximum_flow(
    flow_graph::::IsDirected,                   # the input graph
    source::Integer,                       # the source vertex
    target::Integer,                       # the target vertex
    capacity_matrix::AbstractMatrix,   # edge flow capacities
    algorithm::DinicAlgorithm              # keyword argument for algorithm
    )
    residual_graph = residual(flow_graph)
    return dinic_impl(residual_graph, source, target, capacity_matrix)
end

# Method for Boykov-Kolmogorov algorithm

@traitfn function maximum_flow(
    flow_graph::::IsDirected,                   # the input graph
    source::Integer,                       # the source vertex
    target::Integer,                       # the target vertex
    capacity_matrix::AbstractMatrix,   # edge flow capacities
    algorithm::BoykovKolmogorovAlgorithm   # keyword argument for algorithm
    )
    residual_graph = residual(flow_graph)
    return boykov_kolmogorov_impl(residual_graph, source, target, capacity_matrix)
end

# Method for Push-relabel algorithm

@traitfn function maximum_flow(
    flow_graph::::IsDirected,                   # the input graph
    source::Integer,                       # the source vertex
    target::Integer,                       # the target vertex
    capacity_matrix::AbstractMatrix,   # edge flow capacities
    algorithm::PushRelabelAlgorithm        # keyword argument for algorithm
    )
    residual_graph = residual(flow_graph)
    return push_relabel(residual_graph, source, target, capacity_matrix)
end

"""
    maximum_flow(flow_graph, source, target[, capacity_matrix][, algorithm][, restriction])

Generic maximum_flow function for `flow_graph` from `source` to `target` with
capacities in `capacity_matrix`.
Uses flow algorithm `algorithm` and cutoff restriction `restriction`.

- If `capacity_matrix` is not specified, `DefaultCapacity(flow_graph)` will be used.
- If `algorithm` is not specified, it will default to [`PushRelabelAlgorithm`](@ref).
- If `restriction` is not specified, it will default to `0`.

Return a tuple of (maximum flow, flow matrix). For the Boykov-Kolmogorov
algorithm, the associated mincut is returned as a third output.

### Usage Example:

```jldoctest
julia> flow_graph = DiGraph(8) # Create a flow-graph
julia> flow_edges = [
(1,2,10),(1,3,5),(1,4,15),(2,3,4),(2,5,9),
(2,6,15),(3,4,4),(3,6,8),(4,7,16),(5,6,15),
(5,8,10),(6,7,15),(6,8,10),(7,3,6),(7,8,10)
]

julia> capacity_matrix = zeros(Int, 8, 8)  # Create a capacity matrix

julia> for e in flow_edges
    u, v, f = e
    add_edge!(flow_graph, u, v)
    capacity_matrix[u,v] = f
end

julia> f, F = maximum_flow(flow_graph, 1, 8) # Run default maximum_flow without the capacity_matrix

julia> f, F = maximum_flow(flow_graph, 1, 8) # Run default maximum_flow with the capacity_matrix

julia> f, F = maximum_flow(flow_graph,1,8,capacity_matrix,algorithm=EdmondsKarpAlgorithm()) # Run Edmonds-Karp algorithm

julia> f, F = maximum_flow(flow_graph,1,8,capacity_matrix,algorithm=DinicAlgorithm()) # Run Dinic's algorithm

julia> f, F, labels = maximum_flow(flow_graph,1,8,capacity_matrix,algorithm=BoykovKolmogorovAlgorithm()) # Run Boykov-Kolmogorov algorithm

```
"""
function maximum_flow(
    flow_graph::AbstractGraph,                   # the input graph
    source::Integer,                       # the source vertex
    target::Integer,                       # the target vertex
    capacity_matrix::AbstractMatrix =  # edge flow capacities
    DefaultCapacity(flow_graph);
    algorithm::AbstractFlowAlgorithm  =    # keyword argument for algorithm
    PushRelabelAlgorithm(),
    restriction::Real = 0               # keyword argument for restriction max-flow
    )
    if restriction > 0
        return maximum_flow(flow_graph, source, target, min.(restriction, capacity_matrix), algorithm)
    end
    return maximum_flow(flow_graph, source, target, capacity_matrix, algorithm)
end


text:
import Random: randn, randexp, bitrand, randstring, randsubseq, shuffle, randperm, randcycle

@info "Testing Generic functions"
r = PCG.PCGStateOneseq(123)
@test randn(r) == -0.4957408739873887
@test randn(r, Float64, 8) == [
    -0.38659391873108484, -0.12391575881839223, 0.305709506865529, 1.0234147128314572,
    1.1925044460767074, 0.8448640519165571, -0.2982275988025108, 0.9615482011555472
]
@test randexp(r) == 1.0279939998988223
@test bitrand(r, 3, 3) == [false false false; false true false; false true true]
@test randstring(r) == "KuSFMMEc"
a = 1:100
@test randsubseq(r, a, 0.1) == [1, 11, 12, 17, 19, 21, 27, 38, 44, 54, 58, 78, 80]
a = 1:10
@test shuffle(r, a) == [3, 7, 5, 10, 2, 6, 1, 4, 9, 8]
@test randperm(r, 10) == [2, 10, 1, 6, 4, 3, 7, 8, 9, 5]
@test randcycle(r, 10) == [8, 4, 5, 1, 10, 2, 3, 9, 7, 6]


text:
const PyJuliaTypeValue_Type__ref = Ref(PyNULL)
PyJuliaTypeValue_Type() = begin
    ptr = PyJuliaTypeValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaAnyValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.TypeValue",
                base = base,
                as_mapping = (
                    subscript = pyjltype_getitem,
                    ass_subscript = pyjltype_setitem,
                ),
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaTypeValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaTypeValue_New(x::Type) = PyJuliaValue_New(PyJuliaTypeValue_Type(), x)
PyJuliaValue_From(x::Type) = PyJuliaTypeValue_New(x)

pyjltype_getitem(xo::PyPtr, ko::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::Type
    if PyTuple_Check(ko)
        r = PyObject_Convert(ko, Tuple)
        ism1(r) && return PyNULL
        k = takeresult(Tuple)
        try
            PyObject_From(x{k...})
        catch err
            PyErr_SetJuliaError(err)
            PyNULL
        end
    else
        r = PyObject_Convert(ko, Any)
        ism1(r) && return PyNULL
        k = takeresult(Any)
        try
            PyObject_From(x{k})
        catch err
            PyErr_SetJuliaError(err)
            PyNULL
        end
    end
end

pyjltype_setitem(xo::PyPtr, ko::PyPtr, vo::PyPtr) = begin
    PyErr_SetString(PyExc_TypeError(), "Not supported.")
    PyErr()
end


text:
# A simple NN graph implementation

struct DescentGraph{S, V <: AbstractVector{S}, M <: SemiMetric}
    data::Vector{V}
    metric::M
    indices::Matrix{Int}
    distances::Matrix{S}
end

"""
    DescentGraph(data::Vector{V}, n_neighbors::Integer[, metric::SemiMetric = Euclidean()]; <keyword arguments>)

Build an approximate kNN graph of `data` using nearest neighbor descent.
# Arguments
- `max_iters::Integer = 10`: Limits the number of iterations to refine candidate
nearest neighbors. Higher values trade off speed for accuracy. Note that graph
construction may terminate early if little progress is being made.
- `sample_rate::AbstractFloat = 1.`: The sample rate for calculating *local joins*
around each point. Lower values trade off accuracy for speed.
- `precision::AbstractFloat = .001`: The threshold for early termination,
where precision is "roughly the fraction of true kNN allowed to be missed due to
early termination". Lower values take longer but return more accurate results.
"""
function DescentGraph(data::Vector{V},
                      n_neighbors::Integer,
                      metric::SemiMetric = Euclidean();
                      max_iters::Integer = 10,
                      sample_rate::R = 1.,
                      precision::R = 0.001
                     ) where {V <: AbstractVector, R <: AbstractFloat}
    length(data) >= 2 || error("data must contain at least 2 elements")
    n_neighbors <= length(data) - 1 || error("n_neighbors must be 1 less than length(data)=", length(data))
    max_iters >= 1 || error("max_iters must be greater than 0")
    0. < sample_rate ≤ 1. || error("sample_rate must be in (0., 1.]")
    0. ≤ precision ≤ 1. || error("precision must be in [0., 1.]")
    ids, dists = build_graph(data, n_neighbors, metric, max_iters, sample_rate, precision)
    DescentGraph(data, metric, ids, dists)
end

"""
    DescentGraph(data::AbstractMatrix, n_neighbors::Integer[, metric::SemiMetric = Euclidean()]; <keyword arguments>)
"""
function DescentGraph(data::AbstractMatrix,
                      n_neighbors::Integer,
                      metric::SemiMetric = Euclidean();
                      max_iters::Integer = 10,
                      sample_rate::AbstractFloat = 1.,
                      precision::AbstractFloat = 0.001
                     )
    data_vectors = [@view data[:, i] for i in 1:size(data, 2)]
    DescentGraph(data_vectors, n_neighbors, metric; max_iters=max_iters, sample_rate=sample_rate, precision=precision)
end

"""
Return a kNN graph for the input data according to the given metric.
"""
function build_graph(data::Vector{V},
                     k::Integer,
                     metric::M,
                     max_iters::Integer,
                     sample_rate::R,
                     precision::R
                    ) where {V <: AbstractArray,
                             M <: SemiMetric,
                             R <: AbstractFloat}

    np = length(data)
    # initialize with random neighbors
    Dtype = result_type(metric, data[1], data[1])
    knn_heaps = make_knn_heaps(data, k, metric)

    # until no further updates
    for i = 1:max_iters
        # get the fw and bw neighbors of each point
        old_fw, fw, old_bw, bw = _neighbors(knn_heaps, sample_rate)
        old_neighbors = [union(old_fw[i], old_bw[i]) for i in 1:np]
        new_neighbors = [union(fw[i], bw[i]) for i in 1:np]
        c = 0
        # calculate local join around each point
        for i in 1:np
            for u₁ ∈ new_neighbors[i]
                for u₂ ∈ new_neighbors[i]
                    # both points are new
                    if i ≠ u₁ && i ≠ u₂ && u₁ < u₂
                        d = evaluate(metric, data[u₁], data[u₂])
                        c += _heappush!(knn_heaps[u₁], NNTuple(u₂, d))
                        c += _heappush!(knn_heaps[u₂], NNTuple(u₁, d))
                    end
                end
                for u₂ ∈ old_neighbors[i]
                    # one point is new
                    if i ≠ u₁ && i ≠ u₂
                        d = evaluate(metric, data[u₁], data[u₂])
                        c += _heappush!(knn_heaps[u₁], NNTuple(u₂, d))
                        c += _heappush!(knn_heaps[u₂], NNTuple(u₁, d))
                    end
                end
            end

        end
        if c <= precision*k*np
            break
        end
    end

    return deheap_knns(knn_heaps, k)
end

"""
    search(graph::DescentGraph, queries::Vector{V}, n_neighbors::Integer, queue_size=1.) -> indices, distances

Search the kNN `graph` for the nearest neighbors of the points in `queries`.
`queue_size` controls how large the candidate queue should be as a multiple of
`n_neighbors`. Larger values increase accuracy at the cost of speed.
"""
function search(graph::DescentGraph,
                queries::Vector{V},
                n_neighbors::Integer,
                queue_size::Real = 1.,
                ) where {V <: AbstractArray}
    length(queries) ≥ 1 || error("queries must have at least 1 element")
    n_neighbors ≥ 1 || error("n_neighbors must be at least 1")
    queue_size ≥ 1. || error("queue_size must be at least 1.")

    max_candidates = trunc(Int, n_neighbors*queue_size)
    Dtype = result_type(graph.metric, queries[1], queries[1])
    candidates = [BinaryMaxHeap{NNTuple{Int, Dtype}}() for _ in 1:length(queries)]
    for i in eachindex(queries)
        # init
        seen = fill(false, length(graph.data))
        j =  rand(1:length(graph.data))
        d = evaluate(graph.metric, queries[i], graph.data[j])
        _heappush!(candidates[i], NNTuple(j, d, false), max_candidates)
        seen[j] = true

        while true
            unexp = min_flagged(candidates[i])
            if unexp.idx == -1
                break
            end
            # expand closest unexpanded neighbor
            unexp.flag = true
            for k in view(graph.indices, :, unexp.idx)
                if !seen[k]
                    seen[k] = true
                    d = evaluate(graph.metric,
                                 queries[i],
                                 graph.data[k])
                    _unchecked_heappush!(candidates[i],
                                         NNTuple(k, d, false),
                                         max_candidates)
                end
            end
        end
    end
    return deheap_knns(candidates, n_neighbors)
end

"""
    search(graph::DescentGraph, queries::AbstractMatrix, n_neighbors::Integer, queue_size::Real = 1.)
"""
function search(graph::DescentGraph,
                queries::AbstractMatrix,
                n_neighbors::Integer,
                queue_size::Real = 1.,
                )
    queries_array = [@view queries[:, i] for i in 1:size(queries, 2)]
    search(graph, queries_array, n_neighbors, queue_size)
end


text:
module Garlic

using Cropbox

include("atmosphere/atmosphere.jl")
include("rhizosphere/rhizosphere.jl")
include("phenology/phenology.jl")
include("morphology/morphology.jl")
include("physiology/physiology.jl")

include("examples/examples.jl")

end


text:
using VortexDistributions, JLD2, Interpolations;
using NearestNeighbors, LinearAlgebra, Distances;
using GLMakie, Colors

include("../util.jl");
# @load "Vortices3D/data/box_vorts.jld2" psi_chaos1 psi_chaos2 psi_knots1 psi_knots2 psi_knots3 psi_knots4 psi_tubes1 psi_tubes2 psi_ringtube X
@load "data/box_rings.jld2" psi_ring1 psi_ring2 psi_ring3 X # make sure julia repl is in base directory

psi = psi_ring1;
plot_iso(psi,X)

## Find vortex points
N = 100
@time vorts_3d = findvortices3D_itp(psi, X, N); # Find vortex points with interpolation depth N
v_matrix = vcat(vorts_3d'...)[:,1:3]' # Convert to matrix for kdtree

scatterVortsOnIso(vorts_3d)
## Find filaments
x = X[1]; y = X[2]; z = X[3];
Δx = x[2]-x[1]; Δy = y[2]-y[1]; Δz = z[2]-z[1];

α = .2
ϵ = (1+α)*sqrt(Δx^2+Δy^2+Δz^2)/N
@time sets3 = setMethod3(v_matrix, ϵ)
@time setsPeriodic = setMethodPeriodic(v_matrix, X, ϵ, true)

# @time nearest = uniqueVortices(ϵ, psi, X, N)
# @time sets1 = setMethod(psi, X, N, ϵ)
# @time sets2 = setMethod2(v_matrix, ϵ);

plot_iso(psi, X)
using Colors
colors = distinguishable_colors(length(sets3),[RGB(0,0.35,0.25)],dropseed=true)
for i in 1:length(setsPeriodic)
    vi = v_matrix[:, collect(setsPeriodic[i])]
    scatter!(vi[1,:],vi[2,:],vi[3,:],markersize=200,color=colors[i])
end

sortslices(v_matrix, dims=2, by=x->x[1])

A = [v_matrix[:, i] for i in 1:length(v_matrix[1,:]) if v_matrix[1, i] == -7.75 || v_matrix[1,i] == 8.0]

for i in 1:length(A)
    println(A[i])
end

## x[end] points
scatter!([8.0], [-0.875], [1.2947], markersize=400, color="red") #A[end-9]

scatter!([8.0], [-3.34913], [0.741825], markersize=400, color="red") #A[end-11]

## x[1] points

scatter!([-7.75], [-1.03938], [0.9677467], markersize=400, color="blue") #A[9]

scatter!([-7.75], [-3.707936], [0.317925], markersize=400, color="blue") #A[6]

euclidean(A[end-9][2:3], A[9][2:3])
euclidean(A[end-11][2:3], A[6][2:3])

sqrt(3)Δx

@time psi_itp1 = interpolate(psi, BSpline(Cubic(Periodic(OnCell()))));
@time psi_itp2 = interpolate(psi, BSpline(Quadratic(Periodic(OnCell()))));

psi_exp = extrapolate(psi_itp, Periodic(OnCell()))

psi_exp[0.1,0:65,0:65]

psi_itp[1.1,1.1,1.1]
psi_exp[65.1,65.1,65.1]

x_itp = interpolate(X[1], BSpline(Linear()))
x_etp = extrapolate(x_itp, Line())
x_etp[64]

x_etp[1:65]

r1 = LinRange(1,10,10)
r2 = LinRange(10,11,10)
vcat(r1, r2)

X[1]

x = collect(LinRange(1,2,64))
X[1]


text:
export Quad

mutable struct Quad <: AbstractElement
    id::String
    hid::Int

    node1::Node
    node2::Node
    node3::Node
    node4::Node
    material::Material
    t::Float64
    t2::Float64
    t3::Float64
    t4::Float64

    elm_type::String
    mass_type::String

    center::Vector{Float64}
    A::Float64
    T::Matrix{Float64}
end

function Quad(id,hid,node1,node2,node3,node4,material,t,t2=0,t3=0,t4=0,elm_type="TMGQ",mass_type="concentrate")
    #Initialize local CSys,could be optimized by using a MSE plane
    o=(node1.loc+node2.loc+node3.loc+node4.loc)/4
    pt1 = (node2.loc+node3.loc)/2
    pt2 = (node3.loc+node4.loc)/2
    csys = CSys(o, pt1, pt2)

    vertices=[n.loc for n in (node1,node2,node3,node4)]
    #check aspect ratio
    a=vertices[1]-vertices[3]
    b=vertices[2]-vertices[4]
    if !(0.1<norm(a)/norm(b)<10)
        @warn "Quad "*string(id)*" has an aspect ratio larger than 10, may lead to element illness."
    end
    #check wrong sequence
    for i in 1:4
        j=i%4+1
        if (vertices[i]-o)×(vertices[j]-o)⋅csys.T[:,3]<0
            @warn  "Quad "*string(id)*" has a wrong sequence of nodes, may lead to element illness."
            break
        end
    end
    #check devolution vertex
    for i in 1:4
        j=i%4+1
        if norm(vertices[i]-vertices[j])<1e-9
            @warn  "Quad "*string(id)*" has only 3 independent nodes, may lead to element illness."
            break
        end
    end
    #check distorsion
    for i in 1:4
        v=vertices[i]-o
        tol=cos(5*π/180)
        if abs(v⋅csys.T[:,3])/norm(csys.T[:,3])/norm(v)>tol
            @warn  "Quad "*string(id)*" is too distorted, may lead to element illness."
            break
        end
    end
    #check convexness
    for i in 1:4
        j=(i+1)%4+1
        k=i%4+1
        if norm((vertices[i]+vertices[j])/2-o)>=norm(vertices[k]-o)
            @warn  "Quad "*string(id)*" is not convex, may lead to element illness."
            break
        end
    end

    #area is considered as the average of trangles generated by splitting the quand with diagonals
    A=0.5*det([[1;1;1] node1.loc-o node2.loc-o]')
    A+=0.5*det([[1;1;1]  node2.loc-o node3.loc-o]')
    A+=0.5*det([[1;1;1] node3.loc-o node4.loc-o]')
    A+=0.5*det([[1;1;1] node4.loc-o node1.loc-o]')

    #3D to local 2D
    T=zeros(24,24)
    for i in 1:8
        T[3i-2:3i,3i-2:3i]=csys.T
    end

    Quad(id,hid,node1,node2,node3,node4,material,t,t2,t3,t4,elm_type,mass_type,o,A,T)
end

for (root,dirs,files) in walkdir(joinpath(@__DIR__,"quads"))
    for file in files
        if file[end-2:end]==".jl"
            include(joinpath(@__DIR__,"quads",file))
        end
    end
end

function integrateK(elm::Quad)::SparseMatrixCSC{Float64}
    if elm.elm_type=="GQ12"
        K=K_GQ12(elm)
    elseif elm.elm_type=="TMQ"
        K=K_TMQ(elm)
    elseif elm.elm_type=="DKQ"
        K=K_DKQ(elm)
    elseif elm.elm_type=="DKGQ"
        K=K_DKGQ(elm)
    elseif elm.elm_type=="TMGQ"
        K=K_TMGQ(elm)
    else
        throw("Quad element type error!")
    end
    return K
end

function integrateKσ(elm::Quad,σ)::Vector{Float64}
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    x₄,y₄,z₄=T*(elm.node4.loc-center)
    K=Matrix{Float64}(undef,8,8)
    J=Matrix{Float64}(undef,2,2)

    o=(node1.loc+node2.loc+node3.loc+node4.loc)/4
    n1=node1.loc+[u[1];u[2];0]
    n2=node2.loc+[u[3];u[4];0]
    n3=node3.loc+[u[5];u[6];0]
    n4=node4.loc+[u[7];u[8];0]
    pt1 = n1+n2
    pt2 = n2+n3
    csys = CSys(o, pt1, pt2)
    V=csys.T
    T̄ᵉ=zeros(8,8)
    T̄ᵉ[1:2,1:2]=V
    T̄ᵉ[3:4,3:4]=V
    T̄ᵉ[5:6,5:6]=V
    T̄ᵉ[7:8,7:8]=V

    function GtΣG(ξ,η)
        σ[1,1]= G*(-1 + (1/16)*(1 - ξ)^2*(1 - η)^2 + (1/16)*(1 - ξ)^2*(1 + η)^2 + (1/16)*(1 + ξ)^2*(1 - η)^2 + (1/16)*(1 + ξ)^2*(1 + η)^2) + (1/2)*(-1 + (1/16)*(1 - ξ)^2*(1 - η)^2 + (1/16)*(1 - ξ)^2*(1 + η)^2 + (1/16)*(1 + ξ)^2*(1 - η)^2 + (1/16)*(1 + ξ)^2*(1 + η)^2)*λ
        σ[2,2]= G*(-1 + (1/16)*(1 - ξ)^2*(1 - η)^2 + (1/16)*(1 - ξ)^2*(1 + η)^2 + (1/16)*(1 + ξ)^2*(1 - η)^2 + (1/16)*(1 + ξ)^2*(1 + η)^2) + (1/2)*(-1 + (1/16)*(1 - ξ)^2*(1 - η)^2 + (1/16)*(1 - ξ)^2*(1 + η)^2 + (1/16)*(1 + ξ)^2*(1 - η)^2 + (1/16)*(1 + ξ)^2*(1 + η)^2)*λ

        Σᵉ=[σ[1,1] 0;
            0 σ[2,2]]

        G=       [0  (-1/4)*(1 - η)              0  (-1/4)*(1 + η)               0  (1/4)*(1 - η)              0  (1/4)*(1 + η);
     (-1/4)*(1 - ξ)               0  (1/4)*(1 - ξ)               0  (-1/4)*(1 + ξ)              0  (1/4)*(1 + ξ)              0]

        J[1,1]=(-1/4)*(1 - η)*x₁ + (1/4)*(1 - η)*x₃ + (-1/4)*(1 + η)*x₂ + (1/4)*(1 + η)*x₄
        J[1,2]=(-1/4)*(1 - η)*y₂ + (1/4)*(1 - η)*y₃ + (-1/4)*(1 + η)*y₂ + (1/4)*(1 + η)*y₄
        J[2,1]=(-1/4)*(1 - ξ)*x₁ + (1/4)*(1 - ξ)*x₂ + (-1/4)*(1 + ξ)*x₃ + (1/4)*(1 + ξ)*x₄
        J[2,2]=(-1/4)*(1 + ξ)*y₃ + (1/4)*(1 + ξ)*y₄
        Gᵉ'*Σᵉ*Gᵉ*det(J)
    end
    Kσ=hquadrature(GtΣG,[-1,1],[1,1])
end

function integrateKu(elm::Quad,u::Vector{Float64})
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    x₄,y₄,z₄=T*(elm.node4.loc-center)
    K=Matrix{Float64}(undef,8,8)
    J=Matrix{Float64}(undef,2,2)

    o=(node1.loc+node2.loc+node3.loc+node4.loc)/4
    n1=node1.loc+[u[1];u[2];0]
    n2=node2.loc+[u[3];u[4];0]
    n3=node3.loc+[u[5];u[6];0]
    n4=node4.loc+[u[7];u[8];0]
    pt1 = n1+n2
    pt2 = n2+n3
    csys = CSys(o, pt1, pt2)
    V=csys.T
    T̄ᵉ=zeros(8,8)
    T̄ᵉ[1:2,1:2]=V
    T̄ᵉ[3:4,3:4]=V
    T̄ᵉ[5:6,5:6]=V
    T̄ᵉ[7:8,7:8]=V
end

# function integrateK(elm,p_delta=false,σ=nothing,large_deform=false,u=nothing)
#     K=integrateK(elm)
#     if p_delta
#         K+=integrateKσ(elm,σ)
#     elseif large_deform
#         K+=integrateKu(elm,u)
#     end
#     return K
# end

function integrateM(elm::Quad)::SparseMatrixCSC{Float64}
    ρ=elm.material.ρ
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    x₄,y₄,z₄=T*(elm.node4.loc-center)
    M=Matrix{Float64}(undef,24,24)
    J=Matrix{Float64}(undef,2,2)
    function ρNtN(x)
        ξ,η=x[1],x[2]
        # include(joinpath(dirname(@__FILE__),"jacobi.jl"))
        J[1,1]=(-1/4)*(1 - η)*x₁ + (1/4)*(1 - η)*x₃ + (-1/4)*(1 + η)*x₂ + (1/4)*(1 + η)*x₄
        J[1,2]=(-1/4)*(1 - η)*y₂ + (1/4)*(1 - η)*y₃ + (-1/4)*(1 + η)*y₂ + (1/4)*(1 + η)*y₄
        J[2,1]=(-1/4)*(1 - ξ)*x₁ + (1/4)*(1 - ξ)*x₂ + (-1/4)*(1 + ξ)*x₃ + (1/4)*(1 + ξ)*x₄
        J[2,2]=(-1/4)*(1 + ξ)*y₃ + (1/4)*(1 + ξ)*y₄

        M*det(J)
    end
    return hcubature(ρNtN,[-1,1],[1,1])[1]*t
end

function integrateP!(elm::Quad,elm_force)
    if elm.elm_type=="GQ12"
        P=P_GQ12(elm)
    elseif elm.elm_type=="TMQ"
        P=P_TMQ(elm)
    elseif elm.elm_type=="DKQ"
        P=P_DKQ(elm)
    elseif elm.elm_type=="DKGQ"
        P=P_DKGQ(elm)
    elseif elm.elm_type=="TMGQ"
        P=P_TMGQ(elm)
    else
        throw("Quad element type error!")
    end
    elm.Pᵉ=P
end


text:
using RE2
using Test

function collect_eachmatch(re, str; overlap=false)
    [m.match for m in collect(eachmatch(re, str, overlap = overlap))]
end

for f in [matchall #=, collect_eachmatch =#]
    @test f(rr"a?b?", "asbd") == ["a","","b","",""] == f(rr"""a?b?""", "asbd")
    @test f(rr"a?b?", "asbd", overlap=true) == ["a","","b","",""]
    @test f(rr"\w+", "hello", overlap=true) == ["hello","ello","llo","lo","o"]
    @test f(rr".\s", "x \u2200 x \u2203 y") == ["x ", "∀ ", "x ", "∃ "]
    @test f(rr"(\w+)(\s*)", "The dark side of the moon") ==
          ["The ", "dark ", "side ", "of ", "the ", "moon"]
    @test f(rr"", "") == [""]
    @test f(rr"", "", overlap=true) == [""]
    @test f(rr"aa", "aaaa") == ["aa", "aa"]
    @test f(rr"aa", "aaaa", overlap=true) == ["aa", "aa", "aa"]
    @test f(rr"", "aaa") == ["", "", "", ""]
    @test f(rr"", "aaa", overlap=true) == ["", "", "", ""]
    @test f(rr"GCG","GCGCG") == ["GCG"]
    @test f(rr"GCG","GCGCG",overlap=true) == ["GCG","GCG"]
end

# Issue 8278
target = """71.163.72.113 - - [30/Jul/2014:16:40:55 -0700] "GET emptymind.org/thevacantwall/wp-content/uploads/2013/02/DSC_006421.jpg HTTP/1.1" 200 492513 "http://images.search.yahoo.com/images/view;_ylt=AwrB8py9gdlTGEwADcSjzbkF;_ylu=X3oDMTI2cGZrZTA5BHNlYwNmcC1leHAEc2xrA2V4cARvaWQDNTA3NTRiMzYzY2E5OTEwNjBiMjc2YWJhMjkxMTEzY2MEZ3BvcwM0BGl0A2Jpbmc-?back=http%3A%2F%2Fus.yhs4.search.yahoo.com%2Fyhs%2Fsearch%3Fei%3DUTF-8%26p%3Dapartheid%2Bwall%2Bin%2Bpalestine%26type%3Dgrvydef%26param1%3D1%26param2%3Dsid%253Db01676f9c26355f014f8a9db87545d61%2526b%253DChrome%2526ip%253D71.163.72.113%2526p%253Dgroovorio%2526x%253DAC811262A746D3CD%2526dt%253DS940%2526f%253D7%2526a%253Dgrv_tuto1_14_30%26hsimp%3Dyhs-fullyhosted_003%26hspart%3Dironsource&w=588&h=387&imgurl=occupiedpalestine.files.wordpress.com%2F2012%2F08%2F5-peeking-through-the-wall.jpg%3Fw%3D588%26h%3D387&rurl=http%3A%2F%2Fwww.stopdebezetting.com%2Fwereldpers%2Fcompare-the-berlin-wall-vs-israel-s-apartheid-wall-in-palestine.html&size=49.0KB&name=...+%3Cb%3EApartheid+wall+in+Palestine%3C%2Fb%3E...+%7C+Or+you+go+peeking+through+the+%3Cb%3Ewall%3C%2Fb%3E&p=apartheid+wall+in+palestine&oid=50754b363ca991060b276aba291113cc&fr2=&fr=&tt=...+%3Cb%3EApartheid+wall+in+Palestine%3C%2Fb%3E...+%7C+Or+you+go+peeking+through+the+%3Cb%3Ewall%3C%2Fb%3E&b=0&ni=21&no=4&ts=&tab=organic&sigr=13evdtqdq&sigb=19k7nsjvb&sigi=12o2la1db&sigt=12lia2m0j&sign=12lia2m0j&.crumb=.yUtKgFI6DE&hsimp=yhs-fullyhosted_003&hspart=ironsource" "Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/36.0.1985.125 Safari/537.36"""
pat = rr"""([\d\.]+) ([\w.-]+) ([\w.-]+) (\[.+\]) "([^"\r\n]*|[^"\r\n\[]*\[.+\][^"]+|[^"\r\n]+.[^"]+)" (\d{3}) (\d+|-) ("(?:[^"]|\")+)"? ("(?:[^"]|\")+)"?"""
rm = match(pat, target)
@test length(rm.captures) == 9
@test rm.offsets == [1, 15, 17, 19, 49, 133, 137, 144, 1454]

# Issue 9545 (32 bit)
buf = PipeBuffer()
show(buf, rr"")
@test read(buf, String) == "rr\"\""

# see #10994, #11447: RE2 allows NUL chars in the pattern
@test contains("a\0b", Regex2("^a\0b\$"))

# regex match / search string must be a String
@test_throws ArgumentError match(rr"test", GenericString("this is a test"))
@test_throws ArgumentError findfirst(rr"test", GenericString("this is a test"))

# Named subpatterns
let m = match(rr"(?P<a>.)(.)(?P<b>.)", "xyz")
    @test (m[:a], m[2], m["b"]) == ("x", "y", "z")
    @test sprint(show, m) == "Regex2Match(\"xyz\", a=\"x\", 2=\"y\", b=\"z\")"
end

# Backcapture reference in substitution string
@test replace("xaay", rr"x(a+)" => s"!\1!") == "!aa!y"
@test replace("abcde", rr"(..)(?P<byname>d)" => s"\g<byname>xy\\\1") == "adxy\\bce"
@test_throws ErrorException replace("a", rr"(?P<x>)" => s"\g<y>")

# Proper unicode handling
@test  match(rr"∀∀", "∀x∀∀∀").match == "∀∀"

# TODO missing tests for different flags
# TODO test for error exceptions
# TODO test different behavior of RE2 compared to PCRE2


text:
using PCHIPInterpolation
using Interpolations

include("../../src/noise_src_jl/ioaspl.jl")
include("../../src/noise_src_jl/ipnlt.jl")
include("../../src/noise_src_jl/epnl.jl")

# Inputs
n_t = 100
t_o = range(1,100, length=100)
pnlt = 100 .- (range(1, 10, length=100).-9).^2

@time f_ipnlt(t_o, pnlt)
@time f_ipnlt(t_o, pnlt)

@time f_ioaspl(t_o, pnlt)
@time f_ioaspl(t_o, pnlt)

@time f_epnl(t_o, pnlt)
@time f_epnl(t_o, pnlt)

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
import AQFED
using AQFED.Random, Statistics
import AQFED.MonteCarlo: BrownianBridgeConstruction, transform!, transformByDim!

function simulateGBM(rng::AbstractSeq, nSim, nSteps; withBB = false, size=ceil(Int, log2(nSteps))*4)
    tte = 1.0
    genTimes = collect(LinRange(0.0, tte, ceil(Int, nSteps * tte) + 1))
    logpayoffValues = Vector{Float64}(undef, nSim)
    logpayoffValues .= 0.0
    t0 = genTimes[1]
    local payoffValues
    u = Vector{Float64}(undef, nSim)
    z = u

    cache =  AQFED.MonteCarlo.BBCache{Int, Vector{Float64}}(size) #Dict{Int, Vector{Float64}}()
    local bb
    if withBB
        bb = BrownianBridgeConstruction(genTimes[2:end])
    end

    @inbounds for j = 2:length(genTimes)
        t1 = genTimes[j]
        h = t1 - t0
        dim = j - 1
        if withBB
            transformByDim!(bb, rng, 1, dim, z, cache)
        else
            skipTo(rng, dim, 1) #don't use first point for Sobol
            nextn!(rng, dim, u)
            @. z = u * sqrt(h)
        end

        @. logpayoffValues += -h / 2 + z
        if (t1 == tte)
            payoffValues = @. exp(logpayoffValues)
            #payoffValues = @. max(payoffValues, 1/payoffValues)
        end
        t0 = t1
    end
    payoffMean = mean(payoffValues)
    return payoffMean, stdm(payoffValues, payoffMean) / sqrt(length(payoffValues))
end

function simulateGBMIter(rng::AbstractSeq, nSim, nSteps; withBB = false)
    tte = 1.0
    genTimes = collect(LinRange(0.0, tte, ceil(Int, nSteps * tte) + 1))
    #println(length(genTimes) - 1)
    payoffValues = Vector{Float64}(undef, nSim)
    u = Vector{Float64}(undef, length(genTimes) - 1)
    z = u
    w = Vector{Float64}(undef, length(u))
    sqrth = Vector{Float64}(undef, length(u))
    @. sqrth = sqrt(genTimes[2:end] - genTimes[1:end-1])
    local bb
    if withBB
        bb = BrownianBridgeConstruction(genTimes[2:end])
    end

    #mean = 0.0
    @inbounds for i = 1:nSim
        nextn!(rng, z)
        if withBB
            transform!(bb, z, w)
        else
            @. w = z * sqrth
        end
        #println(u," ",z)
        logpayoffValues = 0.0
        t0 = genTimes[1]
        @inbounds for j = 2:length(genTimes)
            t1 = genTimes[j]
            h = t1 - t0
            logpayoffValues += -h / 2 + w[j-1]
            t0 = t1
        end
        payoffValues[i] = exp(logpayoffValues)
        #mean += payoffValue
    end
    #@. payoffValues = exp(payoffValues)
    #payoffMean = mean/nSim
    payoffMean = mean(payoffValues)
    return payoffMean, stdm(payoffValues, payoffMean) / sqrt(length(payoffValues))
end

@testset "GBM-Randomized" begin
    n = 64 * 1024
    ms = [8, 16, 32, 64]
    nSteps = 1000
    smax = 1
    for m in ms
        valArray = Vector{Float64}(undef, m)
        rng =     AQFED.Random.Blabla8(
                [
                    0xB105F00DBAAAAAAD,
                    0xdeadbeefcafebabe,
                    0xFEE1DEADFEE1DEAD,
                    0xdeadbeefcafebabe,
                ],
                UInt32,
            )
        #Chacha8SIMD(UInt32)
        sArray = Vector{Float64}(undef, smax)
    for s=1:smax
            for i=1:m
           #seq = DigitalSobolSeq(nSteps, trunc(Int,n/m), rng)
            seq = ScrambledSobolSeq(nSteps, 1 << 29, Owen(30, ScramblingRngAdapter(rng)))
            value, _ = simulateGBM(seq, trunc(Int, n/m), nSteps)
            valArray[i] = value
            end
            vMean = mean(valArray)
            vErr = stdm(valArray, vMean) / sqrt(length(valArray))
            println(m, " ", vMean, " ",vErr)
            sArray[s] = vErr
        end
    #println(m, " ", mean(sArray), " ", minimum(sArray)," ", median(sArray)," ",maximum(sArray))
    end
end

@testset "GBM" begin
    n = 32 * 2 * 1024 - 1
    nSteps = 1000
    rng = ScrambledSobolSeq(nSteps, n, NoScrambling())
    value, mcerr,mcerr2 = simulateGBMIter(rng, n, nSteps)
    println(typeof(rng.scrambling), " ", value, " ", mcerr, " ",mcerr2)
    @test isapprox(value, 0.968395, atol = 1e-6)

    rng = ScrambledSobolSeq(nSteps, n, NoScrambling())
    value, mcerr = simulateGBM(rng, n, nSteps)
    println(typeof(rng.scrambling), " ", value, " ", mcerr)
    @test isapprox(value, 0.968395, atol = 1e-6)

    rng = ScrambledSobolSeq(nSteps, n, NoScrambling())
    value, mcerr,mcerr2 = simulateGBMIter(rng, n, nSteps, withBB = true)
    println(typeof(rng.scrambling), " ", value, " ", mcerr," ",mcerr2)
    @test isapprox(0.999562, value, atol= 1e-6)

    maxd = 30 #min(32,ceil(Int,log2(n)+10))
    rng = ScrambledSobolSeq(nSteps, n, Owen(maxd, OriginalScramblingRng()))
    value, mcerr,mcerr2 = simulateGBMIter(rng, n, nSteps)
    println(typeof(rng.scrambling), " ", value, " ", mcerr," ",mcerr2)
    #    @test isapprox(value, 1.000178, atol=1e-6)

    rng = ScrambledSobolSeq(
        nSteps,
        n,
        Owen(
            maxd,
            ScramblingRngAdapter(AQFED.Random.Blabla8(
                [
                    0xB105F00DBAAAAAAD,
                    0xdeadbeefcafebabe,
                    0xFEE1DEADFEE1DEAD,
                    0xdeadbeefcafebabe,
                ],
                UInt64,
            )),
        ),
    )
    value, mcerr, mcerr2 = simulateGBMIter(rng, n, nSteps)
    println(typeof(rng.scrambling), " ", value, " ", mcerr," ", mcerr2)
    #    @test isapprox(value, 1.000178, atol=1e-6)
    rng = ScrambledSobolSeq(nSteps, 1 << (maxd - 1), FaureTezuka(OriginalScramblingRng()))
    value, mcerr, mcerr2 = simulateGBMIter(rng, n, nSteps)
    println(typeof(rng.scrambling), " ", value, " ", mcerr," ", mcerr2)
    #    @test isapprox(value, 1.000178, atol=1e-6)

    rng = ScrambledSobolSeq(
        nSteps,
        1 << (maxd - 1),
        FaureTezuka(ScramblingRngAdapter(AQFED.Random.Blabla8(
            [
                0xB105F00DBAAAAAAD,
                0xdeadbeefcafebabe,
                0xFEE1DEADFEE1DEAD,
                0xdeadbeefcafebabe,
            ],
            UInt64,
        ))),
    )
    value, mcerr, mcerr2 = simulateGBMIter(rng, n, nSteps)
    println(typeof(rng.scrambling), " ", value, " ", mcerr," ", mcerr2)

    rng = DigitalSobolSeq(
        nSteps,
        n,
        AQFED.Random.Blabla8(
            [
                0xB105F00DBAAAAAAD,
                0xdeadbeefcafebabe,
                0xFEE1DEADFEE1DEAD,
                0xdeadbeefcafebabe,
            ],
            UInt32,
        ),
    )
    #AQFED.Random.Chacha8SIMD(UInt32))
    value, mcerr, mcerr2 = simulateGBMIter(rng, n, nSteps)
    println(typeof(rng), " ", value, " ", mcerr," ", mcerr2)

    rng = ModuloSobolSeq(
        nSteps,
        n,
        AQFED.Random.Blabla8(
            [
                0xB105F00DBAAAAAAD,
                0xdeadbeefcafebabe,
                0xFEE1DEADFEE1DEAD,
                0xdeadbeefcafebabe,
            ],
            Float64,
        ),
    )
    # AQFED.Random.Chacha8SIMD(Float64))
    value, mcerr, mcerr2 = simulateGBMIter(rng, n, nSteps)
    println(typeof(rng), " ", value, " ", mcerr," ", mcerr2)

end

@testset "ffz" begin
    n = AQFED.Random.ffz(1)
    println(n)
    @test isequal(1, n)
end

@testset "joekuofirst" begin
    nDim = 3
    n = 10
    #found via external Fortran implementation
    ref = [
        [0.0, 0.0, 0.0],
        [0.5, 0.5, 0.5],
        [0.75, 0.25, 0.25],
        [0.25, 0.75, 0.75],
        [0.375, 0.375, 0.625],
        [0.875, 0.875, 0.125],
        [0.625, 0.125, 0.875],
        [0.125, 0.625, 0.375],
        [0.1875, 0.3125, 0.9375],
        [0.6875, 0.8125, 0.4375],
    ]
    sobol = ScrambledSobolSeq(nDim, n, NoScrambling())
    points = zeros(nDim)
    for i = 2:n
        next!(sobol, points)
        println(points)
        for j = 1:nDim
            @test isapprox(ref[i][j], points[j], atol = 1e-15)
        end
    end
    skipTo(sobol, 0)
    for i = 1:n
        next!(sobol, points)
        #println(points)
        for j = 1:nDim
            @test isapprox(ref[i][j], points[j], atol = 1e-15)
        end
    end
end

@testset "joekuofirstbydim" begin
    nDim = 3
    n = 10
    #found via external Fortran implementation
    ref = [
        [0.0, 0.0, 0.0],
        [0.5, 0.5, 0.5],
        [0.75, 0.25, 0.25],
        [0.25, 0.75, 0.75],
        [0.375, 0.375, 0.625],
        [0.875, 0.875, 0.125],
        [0.625, 0.125, 0.875],
        [0.125, 0.625, 0.375],
        [0.1875, 0.3125, 0.9375],
        [0.6875, 0.8125, 0.4375],
    ]
    sobol = ScrambledSobolSeq(nDim, n, NoScrambling())
    points = zeros(n)
    for j = 1:nDim
        skipTo(sobol, j, 0)
        next!(sobol, j, points)
        #println(points)
        for i = 1:n
            @test isapprox(ref[i][j], points[i], atol = 1e-15)
        end
    end
end

@testset "owenfirst" begin
    nDim = 2
    n = 1 << 29
    #found via external Fortran implementation
    ref = [
        [0.102750, 0.861468],
        [0.705979, 0.113136],
        [0.978561, 0.532342],
        [0.329897, 0.282107],
        [0.408047, 0.647400],
        [0.775438, 0.397518],
        [0.534102, 0.943155],
        [0.149600, 0.191606],
        [0.214946, 0.586519],
        [0.601325, 0.337305],
        [0.842920, 0.761363],
        [0.473622, 0.010117],
        [0.270487, 0.933886],
        [0.921027, 0.185250],
        [0.648217, 0.733849],
        [0.043080, 0.482944],
        [0.015017, 0.694562],
        [0.676498, 0.444115],
        [0.887070, 0.894538],
        [0.304228, 0.146482],
        [0.444827, 0.785095],
        [0.871443, 0.033268],
        [0.567125, 0.610189],
        [0.249417, 0.360518],
        [0.181935, 0.981963],
        [0.501550, 0.230932],
        [0.806098, 0.686269],
        [0.377605, 0.436783],
        [0.361989, 0.509107],
        [0.946739, 0.258475],
        [0.735908, 0.838295],
        [0.072550, 0.089444],
        [0.078716, 0.559680],
        [0.729743, 0.309457],
        [0.956690, 0.857568],
        [0.352040, 0.109226],
        [0.398905, 0.962679],
        [0.784797, 0.211146],
        [0.526634, 0.635692],
        [0.156850, 0.385789],
        [0.224332, 0.773072],
        [0.592210, 0.021847],
        [0.850143, 0.566996],
        [0.466128, 0.317766],
        [0.294277, 0.737751],
        [0.897020, 0.486855],
        [0.670332, 0.906549],
        [0.021182, 0.157901],
        [0.053517, 0.882825],
        [0.637781, 0.134757],
    ]
    rng = OriginalScramblingRng()
    sobol = ScrambledSobolSeq(nDim, n, Owen(30, rng))
    points = zeros(nDim)
    for i = 1:length(ref)
        next!(sobol, points)
        #println(points)
        for j = 1:nDim
            @test isapprox(ref[i][j], points[j], atol = 1e-5)
        end
    end
    skipTo(sobol, 0)
    for i = 1:length(ref)
        next!(sobol, points)
        #println(points)
        for j = 1:nDim
            @test isapprox(ref[i][j], points[j], atol = 1e-5)
        end
    end
end

@testset "owenmaxd20" begin
    nDim = 2
    rng = OriginalScramblingRng()
    points = zeros(nDim)
    #found via external Fortran implementation
    eiRefs =
        [0.99712205451487534, 0.99931521505777710, 0.99886019692965533, 0.99998703162784253]
    n = 1024 * 4 + 1
    sobol = ScrambledSobolSeq(nDim, n, Owen(20, rng))
    sum = 0.0
    twoK = 512
    indexTwoK = 1
    for i = 1:n
        f = 1.0
        next!(sobol, points)
        for j = 1:nDim
            f = f * abs(4 * points[j] - 2)
        end
        if (i % twoK == 0)
            println("i=", i)
            println("ei=", (sum / i))
            @test isapprox(eiRefs[indexTwoK], sum / i, atol = 1e-15)
            indexTwoK += 1
            twoK *= 2
        end
        sum += f
    end
end

@testset "owenblabla" begin
    nDim = 2
    rng = ScramblingRngAdapter(AQFED.Random.Blabla8())
    points = zeros(nDim)
    #non reg values
    eiRefs =
        [0.9986280015086777, 0.9989421327867944, 0.9997705799667642, 0.9995579290533669]
    n = 1024 * 4 + 1
    sobol = ScrambledSobolSeq(nDim, n, Owen(20, rng))
    sum = 0.0
    twoK = 512
    indexTwoK = 1
    for i = 1:n
        f = 1.0
        next!(sobol, points)
        for j = 1:nDim
            f = f * abs(4 * points[j] - 2)
        end
        if (i % twoK == 0)
            println("i=", i)
            println("ei=", (sum / i))
            @test isapprox(eiRefs[indexTwoK], sum / i, atol = 1e-15)
            indexTwoK += 1
            twoK *= 2
        end
        sum += f
    end
end

@testset "fauretezukamaxd20" begin
    nDim = 2
    rng = OriginalScramblingRng()
    points = zeros(nDim)
    n = 1024 * 4
    sobol = ScrambledSobolSeq(nDim, n, FaureTezuka(rng))
    #found via external Fortran implementation
    firstNumbers = [
        [0.1610107421875, 0.0618896484375],
        [0.6610107421875, 0.5618896484375],
        [0.9110107421875, 0.3118896484375],
        [0.4110107421875, 0.8118896484375],
        [0.7860107421875, 0.9368896484375],
        [0.2860107421875, 0.4368896484375],
        [0.0360107421875, 0.6868896484375],
        [0.5360107421875, 0.1868896484375],
        [0.2235107421875, 0.9993896484375],
    ]
    eiRefs =
        [0.99863761523738503, 0.99997501377947628, 0.99842982890550047, 0.99996583949541673]
    sum = 0.0
    twoK = 512
    indexTwoK = 1
    for i = 1:n
        f = 1.0
        next!(sobol, points)
        if (i < 10)
            for j = 1:nDim
                @test isapprox(firstNumbers[i][j], points[j], atol = 1e-15)
            end
        end
        for j = 1:nDim
            f = f * abs(4 * points[j] - 2)
        end
        if (i % twoK == 0)
            println("i=", i)
            println("ei=", (sum / i))
            @test isapprox(eiRefs[indexTwoK], sum / i, atol = 1e-15)
            indexTwoK += 1
            twoK *= 2
        end
        sum += f
    end
end

@testset "owenfauretezukafirst" begin
    nDim = 2
    rng = OriginalScramblingRng()
    points = zeros(nDim)
    n = 1 << 29
    sobol = ScrambledSobolSeq(nDim, n, OwenFaureTezuka(30, rng))
    #found via external Fortran implementation
    ref = [
        [0.171337, 0.890245],
        [0.520179, 0.139048],
        [0.793268, 0.715480],
        [0.398030, 0.466157],
        [0.723320, 0.805123],
        [0.093221, 0.056202],
        [0.351145, 0.604769],
        [0.965125, 0.354211],
        [0.852922, 0.062542],
        [0.455807, 0.814408],
        [0.228854, 0.268739],
        [0.579604, 0.518367],
        [0.283907, 0.242083],
        [0.899795, 0.990335],
        [0.657731, 0.411004],
        [0.025754, 0.661258],
        [0.251729, 0.967130],
        [0.931755, 0.218420],
        [0.627906, 0.638114],
        [0.055796, 0.387280],
        [0.822363, 0.853245],
        [0.486095, 0.101836],
        [0.196431, 0.557142],
        [0.612299, 0.308095],
        [0.695144, 0.017341],
        [0.121127, 0.765866],
        [0.317320, 0.315412],
        [0.999222, 0.565451],
        [0.136776, 0.162291],
        [0.554521, 0.913885],
        [0.764845, 0.489338],
        [0.426669, 0.739181],
        [0.387960, 0.933425],
        [0.803337, 0.181814],
        [0.514136, 0.727452],
        [0.177378, 0.477631],
        [0.943946, 0.761955],
        [0.372326, 0.013440],
        [0.068013, 0.592802],
        [0.748529, 0.342749],
        [0.573562, 0.113545],
        [0.234897, 0.864975],
        [0.445738, 0.288573],
        [0.862992, 0.537603],
        [0.000546, 0.191084],
        [0.682939, 0.939780],
        [0.878615, 0.391182],
        [0.305087, 0.642026],
        [0.035106, 0.978609],
        [0.648597, 0.230371],
    ]
    for i = 1:length(ref)
        next!(sobol, points)
        for j = 1:nDim
            @test isapprox(ref[i][j], points[j], atol = 1e-5)
        end
    end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libexif"
version = v"0.6.21"

# Collection of sources required to complete build
sources = [
    "https://jaist.dl.sourceforge.net/project/libexif/libexif/0.6.21/libexif-0.6.21.tar.gz" =>
    "edb7eb13664cf950a6edd132b75e99afe61c5effe2f16494e6d27bc404b287bf",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libexif-*/
if [[ "${target}" == powerpc64le-* ]] || [[ "${target}" == *-freebsd* ]]; then
    # Install `autopoint` and other tools needed by `autoreconf`
    apk add gettext-dev
    # Rebuild the configure script to convince it to build the shared library
    autoreconf -vi
fi
update_configure_scripts
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libexif", :libexif)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import Oceananigans.Grids: required_halo_size
using Oceananigans.Utils: prettysummary

"""
    struct ScalarBiharmonicDiffusivity{F, N, K} <: AbstractScalarBiharmonicDiffusivity{F}

Holds viscosity and diffusivities for models with prescribed isotropic diffusivities.
"""
struct ScalarBiharmonicDiffusivity{F, N, K} <: AbstractScalarBiharmonicDiffusivity{F}
    ν :: N
    κ :: K

    function ScalarBiharmonicDiffusivity{F}(ν::N, κ::K) where {F, N, K}
        return new{F, N, K}(ν, κ)
    end
end

# Aliases that allow specify the floating type, assuming that the discretization is Explicit in time
          ScalarBiharmonicDiffusivity(FT::DataType;         kwargs...) = ScalarBiharmonicDiffusivity(ThreeDimensionalFormulation(), FT; kwargs...)
  VerticalScalarBiharmonicDiffusivity(FT::DataType=Float64; kwargs...) = ScalarBiharmonicDiffusivity(VerticalFormulation(), FT; kwargs...)
HorizontalScalarBiharmonicDiffusivity(FT::DataType=Float64; kwargs...) = ScalarBiharmonicDiffusivity(HorizontalFormulation(), FT; kwargs...)

required_halo_size(::ScalarBiharmonicDiffusivity) = 2

"""
    ScalarBiharmonicDiffusivity([formulation=ThreeDimensionalFormulation(), FT=Float64;]
                                ν=0, κ=0,
                                discrete_form = false)

Return a scalar biharmonic diffusivity turbulence closure with viscosity coefficient `ν` and tracer
diffusivities `κ` for each tracer field in `tracers`. If a single `κ` is provided, it is applied to
all tracers. Otherwise `κ` must be a `NamedTuple` with values for every tracer individually.

Arguments
=========

* `formulation`:
  - `HorizontalFormulation()` for diffusivity applied in the horizontal direction(s)
  - `VerticalFormulation()` for diffusivity applied in the vertical direction,
  - `ThreeDimensionalFormulation()` (default) for diffusivity applied isotropically to all directions

* `FT`: the float datatype (default: `Float64`)

Keyword arguments
=================

  - `ν`: Viscosity. `Number`, `AbstractArray`, or `Function(x, y, z, t)`.

  - `κ`: Diffusivity. `Number`, `AbstractArray`, or `Function(x, y, z, t)`, or
         `NamedTuple` of diffusivities with entries for each tracer.

  - `discrete_form`: `Boolean`.
"""
function ScalarBiharmonicDiffusivity(formulation=ThreeDimensionalFormulation(), FT=Float64;
                                     ν=0, κ=0,
                                     discrete_form = false)

    ν = convert_diffusivity(FT, ν; discrete_form)
    κ = convert_diffusivity(FT, κ; discrete_form)
    return ScalarBiharmonicDiffusivity{typeof(formulation)}(ν, κ)
end

function with_tracers(tracers, closure::ScalarBiharmonicDiffusivity{F}) where {F}
    κ = tracer_diffusivities(tracers, closure.κ)
    return ScalarBiharmonicDiffusivity{F}(closure.ν, κ)
end

@inline viscosity(closure::ScalarBiharmonicDiffusivity, K) = closure.ν
@inline diffusivity(closure::ScalarBiharmonicDiffusivity, K, ::Val{id}) where id = closure.κ[id]

calculate_diffusivities!(diffusivities, closure::ScalarBiharmonicDiffusivity, args...) = nothing

function Base.summary(closure::ScalarBiharmonicDiffusivity)
    F = summary(formulation(closure))
    return string("ScalarBiharmonicDiffusivity{$F}(ν=", prettysummary(closure.ν), ", κ=", prettysummary(closure.κ), ")")
end

Base.show(io::IO, closure::ScalarBiharmonicDiffusivity) = print(io, summary(closure))


text:
using Plots
using OceanParameterizations
using WindMixing

reconstruct_fluxes = false
println("Reconstruct fluxes? $(reconstruct_fluxes)")

enforce_surface_fluxes = true
println("Enforce surface fluxes? $(enforce_surface_fluxes)")

output_directory=pwd()*"/ProfileVisuals/reconstruct_$(reconstruct_fluxes)/enforce_surface_fluxes_$(enforce_surface_fluxes)"
mkpath(output_directory)

files =  ["free_convection", "strong_wind", "strong_wind_no_coriolis", "weak_wind_strong_cooling",
          "strong_wind_weak_cooling", "strong_wind_weak_heating"]

file_labels = Dict(
    "free_convection" => "Free convection",
    "strong_wind" => "Strong wind",
    "strong_wind_no_coriolis" => "Strong wind, no rotation",
    "weak_wind_strong_cooling" => "Weak wind, strong cooling",
    "strong_wind_weak_cooling" => "Strong wind, weak cooling",
    "strong_wind_weak_heating" => "Strong wind, weak heating"
)

Ts = Dict()

for file in files
    Ts[file] = data(file, reconstruct_fluxes=reconstruct_fluxes,
                    enforce_surface_fluxes=enforce_surface_fluxes) # <: OceananigansData
end

x_lims = Dict(
    "uw" => (-10,0),
    "vw" => (-4,4.5),
    "wT" => (-1.5,0.7),
    "T" => (19.6,20)
)

f = Dict(
    "u" => 𝒟 -> 𝒟.u.coarse,
    "v" => 𝒟 -> 𝒟.v.coarse,
    "uw" => 𝒟 -> 𝒟.uw.coarse,
    "vw" => 𝒟 -> 𝒟.vw.coarse,
    "wT" => 𝒟 -> 𝒟.wT.coarse,
    "T"  => 𝒟 -> 𝒟.T.coarse
)

zs = Dict(
    "u" => 𝒟 -> 𝒟.u.z,
    "v" => 𝒟 -> 𝒟.v.z,
    "uw" => 𝒟 -> 𝒟.uw.z,
    "vw" => 𝒟 -> 𝒟.vw.z,
    "wT" => 𝒟 -> 𝒟.wT.z,
    "T"  => 𝒟 -> 𝒟.T.z
)

# legend_placement = Dict(
#     "uw" => :bottomleft,
#     "vw" => :bottomright,
#     "wT" => :right
# )

legend_placement = Dict(
    "u" => false,
    "v" => false,
    "T" => false,
    "uw" => false,
    "vw" => false,
    "wT" => false,
    "T"  => false,
)

scaling_factor = Dict(
    "uw" => 1e4,
    "vw" => 1e4,
    "wT" => 1e4,
    "u" => 1,
    "v" => 1,
    "T" => 1
)

x_labels = Dict(
    "u" => "U (m/s)",
    "v" => "V (m/s)",
    "uw" => "U'W' x 10⁴ (m²/s²)",
    "vw" => "V'W' x 10⁴ (m²/s²)",
    "wT" => "W'T' x 10⁴ (C⋅m/s)",
    "T" => "T (C)"
)

titles = Dict(
    "uw" => "Zonal momentum flux, U'W'",
    "vw" => "Meridional momentum flux, V'W'",
    "wT" => "Temperature flux, W'T'",
    "u" => "Zonal momentum (m/s)",
    "v" => "Meridional momentum (m/s)",
    "T" => "Temperature, T",
)

function plot_frame_i(name, i)
    p = plot(xlabel=x_labels[name], ylabel="Depth (m)", palette=:Paired_6, legend=legend_placement[name], foreground_color_grid=:white, plot_titlefontsize=20)
    # p = plot(xlabel=x_labels[name], xlims = x_lims[name], ylabel="Depth (m)", palette=:Paired_6, legend=legend_placement[name], foreground_color_grid=:white, plot_titlefontsize=20)
    for (file, T) in Ts
        plot!(f[name](T)[:,i].*scaling_factor[name], zs[name](T), title = titles[name], label="$(file)", linewidth=3)
    end
    plot!(size=(400,500))
    p
end

p1 = plot_frame_i("uw", 288)
savefig(p1, output_directory*"/uw_last_frame.pdf")

p2 = plot_frame_i("vw", 288)
savefig(p2, output_directory*"/vw_last_frame.pdf")

p3 = plot_frame_i("wT", 288)
savefig(p3, output_directory*"/wT_last_frame.pdf")

pu = plot_frame_i("u", 288)
savefig(pu, output_directory*"/u_last_frame.pdf")

pv = plot_frame_i("v", 288)
savefig(pv, output_directory*"/v_last_frame.pdf")

pT = plot_frame_i("T", 288)
savefig(pT, output_directory*"/T_last_frame.pdf")

p4 = plot(grid=false, showaxis=false, palette=:Paired_6, ticks=nothing)
for (file, T) in Ts
    plot!(1, label=file_labels[file], legend=:left, size=(200,600))
end
p4
savefig(p4, output_directory*"/legend_last_frame.pdf")

layout = @layout [a b c d]
p = plot(p1,p2,p3,p4,layout=layout, size=(1600,400), tickfontsize=12)
savefig(p, output_directory*"/all_last_frame_new_suite.pdf")

layout = @layout [a b c d e]
p = plot(p1,p2,p3,pT,p4,layout=layout, size=(1600,400), tickfontsize=12)
savefig(p, output_directory*"/all_last_frame_new_suite_with_T.pdf")

## ANIMATION

# add x_lims for animations
function plot_frame_i(name, i)
    p = plot(xlabel=x_labels[name], xlims = x_lims[name], ylabel="Depth (m)", palette=:Paired_6, legend=legend_placement[name], foreground_color_grid=:white, plot_titlefontsize=20)
    for (file, T) in Ts
        plot!(f[name](T)[:,i].*scaling_factor[name], zs[name](T), title = titles[name], label="$(file)", linewidth=3)
    end
    plot!(size=(400,500))
    p
end

function animate_all(name, Ts)
    anim = @animate for i in 1:288
        plot_frame_i(name, i)
    end
    return anim
end

save_animation(anim, filename) = gif(anim, filename, fps=20)
save_animation(animate_all("uw", Ts), output_directory*"/uw.gif")
save_animation(animate_all("vw", Ts), output_directory*"/vw.gif")
save_animation(animate_all("wT", Ts), output_directory*"/wT.gif")
save_animation(animate_all("T", Ts), output_directory*"/T.gif")


text:
# -------------------------------- #
# Utilities to handle state vector #
# -------------------------------- #
abstract AbstractState{T}

"""
State vector at time t. This is the state vector that can be used to evaluate
the value functions and only contains elements at a single time period.

This is NOT the complete markets state vector used by the planner when updating
the policy function. That vector includes exogenous states at time t+1
"""
type TimeTState{T} <: AbstractState{T}
    l♠::T
    lzbar::T
end

"""
Full state vector used by the planner to apply the policy rules.

Includes endogenous and exogenous states at time t as well as exogenous states
at time t+1. We include the t+1 exog states because we have assumed complete
markets, so the planner should be able to choose the allocation based on all
possible shocks.

Notice that we force all arguments to be of the same type. We have constructors
to do conversions as needed
"""
type FullState{T} <: AbstractState{T}
    l♠::T
    lzbar::T
    lzbarp::T
    lgp::T
end

@inline function FullState(st::TimeTState, lzbarp::Vector, lgp::Vector)
    FullState(st.l♠, st.lzbar, lzbarp, lgp)
end

@inline function FullState(l♠::Real, lzbar::Real, lzbarp::Vector, lgp::Vector)
    N = length(lzbarp)
    FullState(fill(l♠, N), fill(lzbar, N), lzbarp, lgp)
end

# --------------------------- #
# Interface for AbstractState #
# --------------------------- #

# # Define the iterator protocol for the types
for T in (:TimeTState, :FullState)
    @eval begin
        Base.getindex{T<:Real}(st::$T{Vector{T}}, i::Int) =
            $T([getfield(st, nm)[i] for nm in fieldnames(st)]...)
        Base.length(x::$T) = length(x.l♠)
        Base.size(x::$T) = (length(x.l♠), length(fieldnames(x)))
        Base.start{T<:Real}(x::$T{Vector{T}}) = 1
        Base.next{T<:Real}(x::$T{Vector{T}}, i::Int) = (x[i], i+1)
        Base.done{T<:Real}(x::$T{Vector{T}}, i::Int) = i == length(x)+1
        Base.eltype{T<:Real}(::$T{Vector{T}}) = $T{T}
    end
end

# for easy conversion to Vector and Matrix
Base.convert{T<:Real}(::Type{Vector}, st::AbstractState{T}) =
    vcat([getfield(st, nm) for nm in fieldnames(st)]...)
Base.convert(::Type{Matrix}, st::AbstractState) =
    hcat([getfield(st, nm) for nm in fieldnames(st)]...)

# Convert a FullState to a TimeTState -- just extract first to fields
TimeTState(fst::FullState) = TimeTState(fst.l♠, fst.lzbar)

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
#using Test

#=
Test mcrnn operations:
=#

@testset "MaskRCNN" begin
end


text:
"""
HTS_Engine is the main interface for hts_engine_API
"""
mutable struct HTS_Engine
    condition::HTS_Condition
    audio::HTS_Audio
    ms::HTS_ModelSet
    label::HTS_Label
    sss::HTS_SStreamSet
    pss::HTS_PStreamSet
    gss::HTS_GStreamSet
    function HTS_Engine()
        p = new(HTS_Condition(), HTS_Audio(), HTS_ModelSet(), HTS_Label(),
                HTS_SStreamSet(), HTS_PStreamSet(), HTS_GStreamSet())
        initialize(p)
        finalizer(p, clear)
        p
    end
end

function HTS_Engine(voices)
    engine = HTS_Engine()
    return load(engine, voices)
end

for name in [:initialize, :refresh, :clear]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine)
            @htscall($fsymbol, Void, (Ptr{HTS_Engine},), &engine)
        end
    end
end

### Setter ###

# TODO: remove dupulicate code

for (name, argtype) in [
                         (:set_sampling_frequency, Csize_t),
                         (:set_fperiod, Csize_t),
                         (:set_audio_buff_size, Csize_t),
                         (:set_stop_flag, HTS_Boolean),
                         (:set_volume, Cdouble),
                         (:set_speed, Cdouble),
                         (:set_phoneme_alignment_flag, HTS_Boolean),
                         (:set_alpha, Cdouble),
                         (:set_beta, Cdouble),
                         (:add_half_tone, Cdouble)
                         ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, val)
            @htscall($fsymbol, Void, (Ptr{HTS_Engine}, $argtype), &engine, val)
        end
    end
end

### Setter for stream

for (name, argt1, argt2) in [
                             (:set_msd_threshold, Csize_t, Cdouble),
                             (:set_gv_weight, Csize_t, Cdouble)
                             ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, stream_index, val)
            @assert stream_index > 0 && stream_index <= get_nstream(engine)
            @htscall($fsymbol, Void, (Ptr{HTS_Engine}, $argt1, $argt2),
                     &engine, stream_index-1, val)
        end
    end
end

### Setter for voice

for (name, argt1, argt2) in [
                             (:set_duration_interpolation_weight, Csize_t, Cdouble)
                             ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, voice_index, val)
            @assert voice_index > 0 && voice_index <= get_nvoices(engine)
            @htscall($fsymbol, Void, (Ptr{HTS_Engine}, $argt1, $argt2),
                     &engine, voice_index-1, val)
        end
    end
end

### Setter for stream and voice

for (name, argt1, argt2, argt3) in [
                                    (:set_parameter_interpolation_weight,
                                     Csize_t, Csize_t, Cdouble),
                                    (:set_gv_interpolation_weight,
                                     Csize_t, Csize_t, Cdouble)
                                    ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, voice_index, stream_index, val)
            @assert voice_index > 0 && voice_index <= get_nvoices(engine)
            @assert stream_index > 0 && stream_index <= get_nstream(engine)
            @htscall($fsymbol, Void, (Ptr{HTS_Engine}, $argt1, $argt2, $argt3),
                     &engine, voice_index-1, stream_index-1, val)
        end
    end
end

### Getter ###

for (name, rettype) in [
                        (:get_sampling_frequency, Csize_t),
                        (:get_fperiod, Csize_t),
                        (:get_audio_buff_size, Csize_t),
                        (:get_volume, Cdouble),
                        (:get_alpha, Cdouble),
                        (:get_beta, Cdouble),
                        (:get_total_state, Csize_t),
                        (:get_nvoices, Csize_t),
                        (:get_nstream, Csize_t),
                        (:get_nstate, Csize_t),
                        (:get_total_frame, Csize_t),
                        (:get_nsamples, Csize_t)
                        ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine)
            r = @htscall($fsymbol, $rettype, (Ptr{HTS_Engine},), &engine)
            @signed r # for convenience
        end
    end
end

###  Getter for stream

for (name, rettype, argtype) in [
                                 (:get_msd_threshold, Cdouble, Csize_t),
                                 (:get_gv_weight, Cdouble, Csize_t)
                                 ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, stream_index)
            @assert stream_index > 0 && stream_index <= get_nstream(engine)
            @htscall($fsymbol, $rettype, (Ptr{HTS_Engine}, $argtype),
                     &engine, stream_index-1)
        end
    end
end

###  Getter for voice

for (name, rettype, argtype) in [
                                 (:get_duration_interpolation_weight,
                                  Cdouble, Csize_t)
                                 ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, voice_index)
            @assert voice_index > 0 && voice_index <= get_nvoices(engine)
            @htscall($fsymbol, $rettype, (Ptr{HTS_Engine}, $argtype),
                     &engine, voice_index-1)
        end
    end
end

###  Getter for state

# NOTE: get_state_duration must be called after synthesis
for (name, rettype, argtype) in [
                                 (:get_state_duration, Csize_t, Csize_t)
                                 ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, state_index)
            @assert state_index > 0 && state_index <= get_nstate(engine)
            @htscall($fsymbol, $rettype, (Ptr{HTS_Engine}, $argtype),
                     &engine, state_index-1)
        end
    end
end

### Getter for voice and stream

for (name, rettype, argt1, argt2) in [
                                      (:get_parameter_interpolation_weight,
                                       Cdouble, Csize_t, Csize_t),
                                      (:get_gv_interpolation_weight,
                                       Cdouble, Csize_t, Csize_t)
                                      ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, voice_index, stream_index)
            @assert voice_index > 0 && voice_index <= get_nvoices(engine)
            @assert stream_index > 0 && stream_index <= get_nstream(engine)
            @htscall($fsymbol, $rettype, (Ptr{HTS_Engine}, $argt1, $argt2),
                     &engine, voice_index-1, stream_index-1)
        end
    end
end

function load{T<:AbstractString}(engine::HTS_Engine, voices::Vector{T})
    for voice in voices
        if !isfile(voice)
            error("$voice doesn't exists")
        end
    end
    c_voices = Vector{Vector{Cchar}}(length(voices))
    for i in 1:length(voices)
        c_voices[i] = collect(voices[i])
    end
    ret = @htscall(:HTS_Engine_load, HTS_Boolean,
                 (Ptr{HTS_Engine}, Ptr{Ptr{Cchar}}, Csize_t),
                 &engine, c_voices, length(c_voices))
    success = convert(Bool, ret)
    if !success
        error("failed to load hts voice(s)")
    end

    engine
end

function load{T<:AbstractString}(engine::HTS_Engine, voice::T)
    load(engine, [voice])
end

function get_stop_flag(engine::HTS_Engine)
    stop_flag = @htscall(:HTS_Engine_get_stop_flag, HTS_Boolean,
                         (Ptr{HTS_Engine},), &engine)
    convert(Bool, stop_flag)
end

function get_fullcontext_label_format(engine::HTS_Engine)
    ret = @htscall(:HTS_Engine_get_fullcontext_label_format,
                   Ptr{Cchar}, (Ptr{HTS_Engine},), &engine)
    unsafe_string(ret)
end

function get_fullcontext_label_version(engine::HTS_Engine)
    ret = @htscall(:HTS_Engine_get_fullcontext_label_version,
                   Ptr{Cchar}, (Ptr{HTS_Engine},), &engine)
    unsafe_string(ret)
end

function get_generated_parameter(engine::HTS_Engine,
                                 stream_index, frame_index, vector_index)
    @htscall(:HTS_Engine_get_generated_parameter, Cdouble,
             (Ptr{HTS_Engine}, Csize_t, Csize_t, Csize_t),
             &engine, stream_index, frame_index, vector_index)
end

function get_generated_speech(engine::HTS_Engine, index)
    @assert index > 0 # 1-origin for Julia interface
    @htscall(:HTS_Engine_get_generated_speech, Cdouble,
             (Ptr{HTS_Engine}, Csize_t), &engine, index-1)
end

function get_generated_speech(engine::HTS_Engine)
    speech = Vector{Float64}(get_nsamples(engine))

    @inbounds for i in eachindex(speech)
        speech[i] = get_generated_speech(engine, i)
    end

    speech
end

### Synthesize ###

function synthesize_from_fn(engine::HTS_Engine, labelpath)
    if !isfile(labelpath)
        error("Input lable file doesn't exists")
    end
    ret = @htscall(:HTS_Engine_synthesize_from_fn, HTS_Boolean,
                   (Ptr{HTS_Engine}, Ptr{Cchar}), &engine, labelpath)
    success = convert(Bool, ret)
    if !success
        error("failed to synthesize waveform")
    end
end

function synthesize_from_strings(engine::HTS_Engine, lines)
    clines = Vector{Vector{Cchar}}(length(lines))
    for i in 1:length(lines)
        clines[i] = collect(lines[i])
    end

    ret = @htscall(:HTS_Engine_synthesize_from_strings, HTS_Boolean,
                   (Ptr{HTS_Engine}, Ptr{Ptr{Cchar}}, Csize_t),
                   &engine, clines, length(clines))
    success = convert(Bool, ret)
    if !success
        error("failed to synthesize waveform")
    end
end

### Save ###

for (name, mode) = [
                    (:save_information, "wt"),
                    (:save_label, "wt"),
                    (:save_generated_speech, "wb"),
                    (:save_riff, "wb")
                    ]
    fsymbol = QuoteNode(Symbol(:HTS_Engine_, name))
    @eval begin
        function $name(engine::HTS_Engine, path)
            fp = ccall(:fopen, Ptr{Void}, (Ptr{Cchar}, Ptr{Cchar}),
                       path, $mode)
            @assert fp != C_NULL
            @htscall($fsymbol, Void,
                     (Ptr{HTS_Engine}, Ptr{Void}), &engine, fp)
            ccall(:fclose, Void, (Ptr{Void},), fp)
        end
    end
end

function save_generated_parameter(engine::HTS_Engine, stream_index, path)
    @assert stream_index > 0 && stream_index <= get_nstream(engine)
    fp = ccall(:fopen, Ptr{Void}, (Ptr{Cchar}, Ptr{Cchar}), path, "wb")
    @assert fp != C_NULL
    @htscall(:HTS_Engine_save_generated_parameter, Void,
             (Ptr{HTS_Engine}, Csize_t, Ptr{Void}),
             &engine, stream_index-1, fp)
    ccall(:fclose, Void, (Ptr{Void},), fp)
end


text:
# BigWig Reader
# =============

struct Reader <: BioCore.IO.AbstractReader
    stream::IO
    header::BBI.Header
    zooms::Vector{BBI.Zoom}
    summary::BBI.Summary
    # chrom name => (ID, length)
    chroms::Dict{String,Tuple{UInt32,Int}}
    # ID => chrom name
    chrom_names::Dict{UInt32,String}
    # data index
    index::BBI.RTree
end

function Base.eltype(::Type{Reader})
    return Record
end

function BioCore.IO.stream(reader::Reader)
    return reader.stream
end

"""
    BigWig.Reader(input::IO)

Create a reader for bigWig file format.

Note that `input` must be seekable.
"""
function Reader(input::IO)
    # read header
    header = read(input, BBI.Header)
    if header.magic != BBI.WIG_MAGIC
        error("invalid bigWig magic")
    elseif header.version < 3
        error("not a supported version of bigWig")
    end
    # read zoom objects
    zoom_headers = Vector{BBI.ZoomHeader}(undef, header.zoom_levels)
    read!(input, zoom_headers)
    zooms = [BBI.Zoom(input, h, header.uncompress_buf_size) for h in zoom_headers]
    sort!(zooms, by=z->z.header.reduction_level)
    # read summary, B tree, and R tree
    seek(input, header.total_summary_offset)
    summary = read(input, BBI.Summary)
    chromindex = BBI.BTree(input, header.chromosome_tree_offset)
    chroms = Dict(name => (id, Int(len)) for (name, id, len) in BBI.chromlist(chromindex))
    chrom_names = Dict(id => name for (name, (id, _)) in chroms)
    index = BBI.RTree(input, header.full_index_offset)
    return Reader(input, header, zooms, summary, chroms, chrom_names, index)
end

"""
    chromlist(reader::BigWig.Reader)::Vector{Tuple{String,Int}}

Get the `(name, length)` pairs of chromosomes/contigs.
"""
function chromlist(reader::Reader)::Vector{Tuple{String,Int}}
    return sort!([(name, Int(len)) for (name, (id, len)) in reader.chroms], by=x->x[1])
end

"""
    values(reader::BigWig.Reader, interval::Interval)::Vector{Float32}

Get a vector of values within `interval` from `reader`.

This function fills missing values with `NaN32`.
"""
function values(reader::Reader, interval::Interval)
    return values(reader, seqname(interval), leftposition(interval):rightposition(interval))
end

"""
    values(reader::BigWig.Reader, chrom::AbstractString, range::UnitRange)::Vector{Float32}

Get a vector of values within `range` of `chrom` from `reader`.

This function fills missing values with `NaN32`.
"""
function values(reader::Reader, chrom::AbstractString, range::UnitRange)::Vector{Float32}
    values = Vector{Float32}(undef, length(range))
    if isempty(range)
        return values
    end
    fill!(values, NaN32)
    offset = first(range) - 1
    for record in GenomicFeatures.eachoverlap(reader, Interval(chrom, first(range), last(range)))
        rstart = clamp(record.chromstart+1-offset, 1, lastindex(values))
        rend = clamp(record.chromend-offset, 1, lastindex(values))
        values[rstart:rend] .= record.value
    end
    return values
end

# Record
# ------

mutable struct Record
    chromstart::UInt32
    chromend::UInt32
    value::Float32
    header::SectionHeader
    reader::Reader

    function Record(chromstart, chromend, value)
        return new(chromstart, chromend, value)
    end
end

# Iterator
# --------

mutable struct IteratorState
    stream::BufferedStreams.BufferedInputStream
    done::Bool
    header::SectionHeader
    record::Record
    n_sections::UInt64
    current_section::UInt64
    n_records::UInt16
    current_record::UInt16
end

function Base.iterate(reader::Reader)
    seek(reader.stream, reader.header.full_data_offset)
    # this is defined as UInt32 in the spces but actually UInt64
    section_count = read(reader.stream, UInt64)
    # dummy header
    header = SectionHeader(0, 0, 0, 0, 0, 0, 0, 0)
    state = IteratorState(Libz.ZlibInflateInputStream(reader.stream), false, header, Record(), section_count, 0, 0, 0)
    return iterate(reader, state)
end

function Base.iterate(reader::Reader, state)
    advance!(reader, state)
    if state.done
        return nothing
    end
    return copy(state.record), state
end

function advance!(reader::Reader, state::IteratorState)
    # find a section that has at least one record
    while state.current_section < state.n_sections && state.current_record == state.n_records
        state.header = read(state.stream, SectionHeader)
        state.current_section += 1
        state.n_records = state.header.itemcount
        state.current_record = 0
    end
    if state.current_record == state.n_records
        state.done = true
        return state
    end

    # read a new record
    _read!(reader, state, state.record)
    return state
end

function _read!(reader::Reader, state, record::Record)
    @assert !state.done
    header = state.header
    if isbedgraph(header.datatype)
        chromstart = read(state.stream, UInt32)
        chromend   = read(state.stream, UInt32)
    elseif isvarstep(header.datatype)
        chromstart = read(state.stream, UInt32)
        chromend   = chromstart + header.itemspan
    elseif isfixedstep(header.datatype)
        chromstart = state.current_record == 0 ? header.chromstart : state.record.chromstart + header.itemstep
        chromend   = chromstart + header.itemspan
    else
        throw(ArgumentError("invalid data type"))
    end
    value = read(state.stream, Float32)
    record.chromstart = chromstart
    record.chromend = chromend
    record.value = value
    record.header = header
    record.reader = reader
    state.current_record += 1
    return record
end