text:
using CyclicArrays
using Plots

nx=ny=100
nfaces=1;

faces=zeros(nfaces,2,2,4);
faces[1,1,1,:]=[1,1,2,0];
faces[1,1,2,:]=[1,1,1,0];
faces[1,2,1,:]=[1,2,2,0];
faces[1,2,2,:]=[1,2,1,0];

g=CyclicArray(faces);

ind=range(1, 360, length = 360)
x = y = ind/180*pi
f(y,x) = sin(x) + sin(y)

X = repeat(reshape(x, 1, :), length(y), 1)
Y = repeat(y, 1, length(x))

F=f.(Y,X)

#f = sin.(2*x .+ 10 .* sin.(x)) .+ sin.(3*y)

phi=CyclicArray(F,g)
contourf(phi.data, aspect_ratio=1,xlim=(0,360),ylim=(0,360))
xg=CyclicArray(repeat(ind, 1, length(x))',g)
yg=CyclicArray(repeat(ind, 1, length(x)),g)
u=diff(phi,dims=2)
v=-diff(phi,dims=1)

N=1500
np=25
xi=zeros(N+1,np)
yi=zeros(N+1,np)

x0=collect(0:20:490).%100 .+ 170
y0=reshape(repeat(collect(0:20:80),1,5)',25) .+ 160
xi[1,:]=x0;
yi[1,:]=y0;

dt=100
#for t in 1:N
plt=contourf(yg.data[:,1],xg.data[1,:],phi.data, aspect_ratio=1,xlim=(0,360),ylim=(0,360))
scatter!(xi[1,:],yi[1,:],leg=false)
anim = @animate for t ∈ 1:N
  println(t)
  i=Int.(floor.(xg[1,Int.(floor.(xi[t,:]))]))
  j=Int.(floor.(yg[Int.(floor.(yi[t,:])),1]))
  U=[u[i[ii],j[ii]] for ii in 1:np]
  V=[v[i[ii],j[ii]] for ii in 1:np]
  #println(t,' ',U,' ',V)
  xi[t+1,:]=xi[t,:]+U*dt
  yi[t+1,:]=yi[t,:]+V*dt
  xi[t+1,:]=floor.(xg[1,Int.(sign.(xi[t+1,:]).*floor.(abs.(xi[t+1,:])))])+rem.(xi[t+1,:],1)
  yi[t+1,:]=floor.(yg[Int.(sign.(yi[t+1,:]).*floor.(abs.(yi[t+1,:]))),1])+rem.(yi[t+1,:],1)
  plt=contourf(yg.data[:,1],xg.data[1,:],phi.data, aspect_ratio=1,xlim=(0,360),ylim=(0,360))
  scatter!([-1,xi[t+1,:]],[-1,yi[t+1,:]],leg=false, c=:green)
  #display(plt)
end
gif(anim, "RandomFlow_2D.gif", fps = 30)


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

function SCIPincludeConshdlrLinking(scip)
    ccall((:SCIPincludeConshdlrLinking, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsLinking(scip, cons, name, intvar, binvars, vals, nbinvars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsLinking, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cint}, Cint, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, intvar, binvars, vals, nbinvars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicLinking(scip, cons, name, intvar, binvars, vals, nbinvars)
    ccall((:SCIPcreateConsBasicLinking, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cint}, Cint), scip, cons, name, intvar, binvars, vals, nbinvars)
end

function SCIPexistsConsLinking(scip, intvar)
    ccall((:SCIPexistsConsLinking, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, intvar)
end

function SCIPgetConsLinking(scip, intvar)
    ccall((:SCIPgetConsLinking, libscip), Ptr{SCIP_CONS}, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, intvar)
end

function SCIPgetIntvarLinking(scip, cons)
    ccall((:SCIPgetIntvarLinking, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetBinvarsLinking(scip, cons, binvars, nbinvars)
    ccall((:SCIPgetBinvarsLinking, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Cint}), scip, cons, binvars, nbinvars)
end

function SCIPgetNBinvarsLinking(scip, cons)
    ccall((:SCIPgetNBinvarsLinking, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetValsLinking(scip, cons)
    ccall((:SCIPgetValsLinking, libscip), Ptr{Cint}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

ext: jl


text:
# convenience functions for extracting data from DataFrames

"""
    med_and_errs(p) -> (; med, err1_l, err1_u, err2_l, err2_u)
Convenience function for extracting plottable data from a distribution or
an uncertain object created by the packages `MonteCarloMeasurements` or `Measurements`.
Returns the median `med` and the lower `err1_l` and upper `err1_u` standard error
(for 1σ or 68% confidence inteval). `err2_l` and `err2_u` provide the lower and upper error
limits for the 2σ or 95% confidence interval.

**Example:**
```julia
julia> results = [Particles(d) for d in datasets] # Particles[]
julia> res_w_errs = med_and_errs.(results) # Vector of NamedTuple's with standard errors
julia> res_df = DataFrame(res_w_errs) # results as DataFrame with lower an upper error
1×5 DataFrame
 Row │ med      err1_l     err1_u     err2_l    err2_u
     │ Float64  Float64    Float64    Float64   Float64
─────┼────────────────────────────────────────────────────
   1 │ 1.01325  0.0173805  0.0183057  0.034042  0.0366713
```
**Note:** This function is deprecated and will be removed soon. Use [`val_and_errs()`](@ref) instead.
"""
function med_and_errs(p)
    @warn "med_and_errs() is deprecated and will be removed soon. Use val_and_errs() instead!" maxlog=1
    q = pquantile(p, [0.025, 0.16, 0.5, 0.84, 0.975])
    med = q[3]
    err1_l = med - q[2]
    err1_u = q[4] - med
    err2_l = med - q[1]
    err2_u = q[5] - med
    return (; med, err1_l, err1_u, err2_l, err2_u)
end
function med_and_errs(p::Measurements.Measurement)
    @warn "med_and_errs() is deprecated and will be removed soon. Use val_and_errs() instead!" maxlog=1
    med = Measurements.value(p)
    err1_l = err1_u = Measurements.uncertainty(p)
    err2_l = err2_u = 2err1_l
    return (; med, err1_l, err1_u, err2_l, err2_u)
end

"""
    ratio_with_errs(r::RatioBlockingResult)
    -> (;ratio=med, err1_l, err1_u, err2_l, err2_u, f, σ_f, δ_y, k, success)
Convenience function for extracting plottable data from [`RatioBlockingResult`](@ref).
Returns `NamedTuple` with median and standard error of `r` extracted by
[`p_to_errs()`](@ref). See also [`ratio_of_means()`](@ref).

**Example:**
```julia
julia> results = [ratio_of_means(n[i], d[i]; args...) for i in datasets]
julia> res_w_errs = ratio_with_errs.(results) # Vector of NamedTuple's with standard errors
julia> res_df = DataFrame(res_w_errs) # results as DataFrame with lower an upper error
1×10 DataFrame
 Row │ ratio    err1_l     err1_u     err2_l    err2_u     f        σ_f        δ_y        k      success
     │ Float64  Float64    Float64    Float64   Float64    Float64  Float64    Float64    Int64  Bool
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ 1.01325  0.0173805  0.0183057  0.034042  0.0366713  1.01361  0.0181869  0.0128806      2     true
```
**Note:** This function is deprecated and will be removed soon. Use [`val_and_errs()`](@ref)
instead.
"""
function ratio_with_errs(r::RatioBlockingResult)
    med, err1_l, err1_u, err2_l, err2_u = med_and_errs(r.ratio)
    @warn "ratio_with_errs() is deprecated and will be removed soon. Use NamedTuple() instead!" maxlog=1
    return (;
        ratio=med,
        err1_l,
        err1_u,
        err2_l,
        err2_u,
        f=r.f,
        σ_f=r.σ_f,
        δ_y=r.δ_y,
        k=r.k,
        success=r.success,
    )
end

"""
    to_measurement(p::MonteCarloMeasurements.Particles) -> ::Measurements.measurement
Convert an uncertain number from `MonteCarloMeasurements` to `Measurements` format using the median as the central
point. The new `±` boundaries will include the 68% quantile around the median.
"""
function to_measurement(p::MonteCarloMeasurements.Particles)
        q = pquantile(p, [0.16, 0.5, 0.84])
        σ = max(q[3]-q[2],q[2]-q[1]) # take the larger one
        return Measurements.measurement(q[2],σ)
end

"""
    val(x)
Return the best estimate value for an uncertain `x`. Defaults to the `median`
for uncertain `x` represented by a (sampled) distribution. Supports `MonteCarloMeasurements`
and `Measurements`.

See [`errs`](@ref), [`BlockingResult`](@ref), [`RatioBlockingResult`](@ref).
"""
val(x) = x
val(m::Measurements.Measurement) = Measurements.value(m)
val(p::AbstractParticles) = pmedian(p)
val(r::BlockingResult) = r.mean
val(r::RatioBlockingResult) = pmedian(r)

import Measurements: Measurement, uncertainty, value
"""
    val_and_errs(x; n=1, p=nothing, name=:val) -> (;val, val_l, val_u)
Return the median and the lower and upper error bar for the uncertain value `x` as a
`NamedTuple`. This is useful for plotting scripts. The
interval `[val - val_l, val + val_u]` represents the confidence interval at level `n*σ`,
or at probability `p`. Setting `p` overrides `n`. Supports `MonteCarloMeasurements`
and `Measurements`. The names in the `NamedTuple` can be changed with `name`.

**Example:**
```jldoctest
julia> results = [blocking_analysis(i:0.1:2i+20) for i in 1:3]; # mock results

julia> v = val_and_errs.(results, name="res"); # Vector of NamedTuple's with standard errors

julia> DataFrame(v)
3×3 DataFrame
 Row │ res      res_l    res_u
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │    11.5  1.7282   1.7282
   2 │    13.0  1.7282   1.7282
   3 │    14.5  1.78885  1.78885
```

See [`NamedTuple`](@ref), [`val`](@ref), [`errs`](@ref), [`BlockingResult`](@ref),
[`RatioBlockingResult`](@ref).
"""
function val_and_errs(x; name=:val, kwargs...)
    return (; Symbol(name) => x, Symbol(name, :_l) => zero(x), Symbol(name, :_u) => zero(x))
end
function val_and_errs(
    m::T; p=nothing, n=1, name=:val
) where T <:Union{
    Measurement,
    Complex{<:Measurement},
    AbstractParticles,
    Complex{<:AbstractParticles}
}
    return _errs(p, m, n, name)
end
function _errs(::Nothing, m::Measurement, n, name)
    σ = uncertainty(m)
    return (; Symbol(name) => value(m), Symbol(name, :_l) => n*σ, Symbol(name, :_u) => n*σ)
end
function _errs(p, m::Complex, n, name)
    tr = Tuple(_errs(p, real(m), n, name))
    ti = Tuple(_errs(p, imag(m), n, name))
    nt = (;
        Symbol(name) => tr[1] +  ti[1]*im,
        Symbol(name, :_l) => tr[2] +  ti[2]*im,
        Symbol(name, :_u) => tr[3] +  ti[3]*im,
    )
    return nt
end
function _errs(p::AbstractFloat, m::Measurement, _, name)
    d = Normal(value(m), uncertainty(m))
    cp1 = (1-p)/2
    cp2 = 1-cp1
    q1, q2 = quantile(d, cp1), quantile(d, cp2)
    return (;
        Symbol(name) => value(m),
        Symbol(name, :_l) => value(m) - q1,
        Symbol(name, :_u) => q2 - value(m),
    )
end
function _errs(::Nothing, m::AbstractParticles, n, name)
    p = erf(n/√2)
    return _errs(p, m, n, name)
end
function _errs(p::T, m::AbstractParticles, _, name) where T <: AbstractFloat
    cp1 = (1-p)/2
    cp2 = 1-cp1
    q1, q2, q3 = pquantile(m, (cp1, one(T)/2, cp2))
    return (;Symbol(name) => q2, Symbol(name, :_l) => q2 - q1, Symbol(name, :_u) => q3 - q2)
end
function val_and_errs(r::BlockingResult; kwargs...)
    return val_and_errs(measurement(r); kwargs...)
end
val_and_errs(r::RatioBlockingResult; kwargs...) = val_and_errs(r.ratio; kwargs...)

"""
    errs(x; n=1, p=nothing, name=:err) -> (; err_l, err_u)
Return the lower and upper error bar for the uncertain value `x`.

See [`val_and_errs`](@ref).
"""
function errs(args...; name=:err, kwargs...)
    _, err_l, err_u = val_and_errs(args...; kwargs...)
    return (; Symbol(name, :_l) => err_l, Symbol(name, :_u) => err_u)
end

"""
    NamedTuple(x::BlockingResult; n=1, p=nothing, name=:val)
    NamedTuple(x::RatioBlockingResult; n=1, p=nothing, name=:val)
Return a named tuple with value and error bars (see [`val_and_errs`](@ref)) as well
as additional numerical fields relevant for `x`.

**Example:**
```jldoctest
julia> results = [blocking_analysis(i:0.1:2i+20) for i in 1:3]; # mock results

julia> df = NamedTuple.(results, name=:res)|>DataFrame
3×7 DataFrame
 Row │ res      res_l    res_u    res_err_err  res_p_cov  res_k  res_blocks
     │ Float64  Float64  Float64  Float64      Float64    Int64  Int64
─────┼──────────────────────────────────────────────────────────────────────
   1 │    11.5  1.7282   1.7282      0.352767    2.98667      5          13
   2 │    13.0  1.7282   1.7282      0.352767    2.98667      5          13
   3 │    14.5  1.78885  1.78885     0.350823    3.2          5          14
```
```julia-repl
julia> rbs = ratio_of_means(1 .+sin.(1:0.1:11),2 .+sin.(2:0.1:12)); # more mock results

julia> [NamedTuple(rbs),]|>DataFrame
1×9 DataFrame
 Row │ val       val_l      val_u      val_f     val_σ_f    val_δ_y    val_k  val_blocks  val_success
     │ Float64   Float64    Float64    Float64   Float64    Float64    Int64  Int64       Bool
─────┼────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ 0.581549  0.0925669  0.0812292  0.560532  0.0875548  0.0875548      4          12         true

```

See [`val_and_errs`](@ref), [`val`](@ref), [`errs`](@ref), [`BlockingResult`](@ref),
[`RatioBlockingResult`](@ref).
"""
function Base.NamedTuple(r::BlockingResult; name=:val, kwargs...)
    vae =  val_and_errs(r; name, kwargs...)
    br_info = (;
        Symbol(name, :_err_err) => r.err_err,
        Symbol(name, :_p_cov) => r.p_cov,
        Symbol(name, :_k) => r.k,
        Symbol(name, :_blocks) => r.blocks,
    )
    return (; vae..., br_info...)
end
function Base.NamedTuple(r::RatioBlockingResult; name=:val, kwargs...)
    vae =  val_and_errs(r; name, kwargs...)
    br_info = (;
        Symbol(name, :_f) => r.f,
        Symbol(name, :_σ_f) => r.σ_f,
        Symbol(name, :_δ_y) => r.σ_f,
        Symbol(name, :_k) => r.k,
        Symbol(name, :_blocks) => r.blocks,
        Symbol(name, :_success) => r.success,
    )
    return (; vae..., br_info...)
end


text:
function update_and_get_max_abs_diff!(  fssrb::PotentialSimulationSetupRB{T, N1, N2},
                                        depletion_handling::Val{depletion_handling_enabled}, only2d::Val{only_2d} = Val{false}(), is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
                                        use_nthreads::Int = Base.Threads.nthreads())::T where {T, N1, N2, depletion_handling_enabled, only_2d, _is_weighting_potential}
    tmp_potential::Array{T, N2} = copy(fssrb.potential)
    update!(fssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
    max_diff::T = maximum(abs.(tmp_potential - fssrb.potential))
    return max_diff
end

function update_till_convergence!(  fssrb::PotentialSimulationSetupRB{T, N1, N2},
                                    convergence_limit::Real, bias_voltage::AbstractFloat;
                                    n_iterations_between_checks = 500,
                                    depletion_handling::Val{depletion_handling_enabled} = Val{false}(), only2d::Val{only_2d} = Val{false}(), is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
                                    use_nthreads::Int = Base.Threads.nthreads(), max_n_iterations::Int = -1)::Real where {T, N1, N2, depletion_handling_enabled, only_2d, _is_weighting_potential}
    n_iterations::Int = 0
    cf::T = Inf
    cl::T = abs(convergence_limit * bias_voltage) # to get relative change in respect to bias voltage
    prog = ProgressThresh(cl, 0.1, "Convergence: ")
    while cf > cl
        for i in 1:n_iterations_between_checks
            update!(fssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
        end
        cf = update_and_get_max_abs_diff!(fssrb, depletion_handling, only2d, is_weighting_potential, use_nthreads)
        ProgressMeter.update!(prog, cf)
        n_iterations += n_iterations_between_checks
        if max_n_iterations > 0 && n_iterations > max_n_iterations
            @info "Maximum number of iterations reached. (`n_iterations = $(n_iterations)`)"
            break
        end
    end
    ProgressMeter.finish!(prog)
    return cf
end

function check_for_refinement( potential::Array{T, 3}, grid::Grid{T, 3, :cylindrical}, max_diff::AbstractArray, minimum_distance::AbstractArray)::NTuple{3, Vector{Int}} where {T}
    inds_r::Vector{Int} = Int[]
    inds_φ::Vector{Int} = Int[]
    inds_z::Vector{Int} = Int[]

    ax1::DiscreteAxis{T} = grid[1]
    ax2::DiscreteAxis{T} = grid[2]
    ax3::DiscreteAxis{T} = grid[3]
    int1::Interval = ax1.interval
    int2::Interval = ax2.interval
    int3::Interval = ax3.interval

    refinement_value_r::T = max_diff[1]
    refinement_value_φ::T = max_diff[2]
    refinement_value_z::T = max_diff[3]

    minimum_distance_r::T = minimum_distance[1]
    minimum_distance_φ::T = minimum_distance[2]
    minimum_distance_z::T = minimum_distance[3]

    minimum_distance_φ_deg::T = rad2deg(minimum_distance_φ)

    axr::Vector{T} = collect(grid.axes[1])
    axφ::Vector{T} = collect(grid.axes[2])
    axz::Vector{T} = collect(grid.axes[3])
    cyclic::T = grid.axes[2].interval.right

    for iz in 1:size(potential, 3)
        z = axz[iz]
        for iφ in 1:size(potential, 2)
            φ = axφ[iφ]
            isum = iz + iφ
            for ir in 1:size(potential, 1)
                r = axr[ir]

                # r
                if ir != size(potential, 1)
                    if abs(potential[ir + 1, iφ, iz] - potential[ir, iφ, iz]) >= refinement_value_r
                        if axr[ir + 1] - axr[ir] > minimum_distance_r
                            if !in(ir, inds_r) push!(inds_r, ir) end
                        end
                    end
                end

                # φ
                if iφ == size(potential, 2)
                    if typeof(int2).parameters[2] == :open # ugly solution for now. should be done over multiple dispatch..
                        if abs(potential[ir, 1, iz] - potential[ir, iφ, iz]) >= refinement_value_φ
                            if cyclic - axφ[iφ] > minimum_distance_φ
                                if !in(iφ, inds_φ) push!(inds_φ, iφ) end
                            end
                        end
                    end
                else
                    if abs(potential[ir, iφ + 1, iz] - potential[ir, iφ, iz]) >= refinement_value_φ
                        if axφ[iφ + 1] - axφ[iφ] > minimum_distance_φ
                            if !in(iφ, inds_φ) push!(inds_φ, iφ) end
                        end
                    end
                end

                # z
                if iz != size(potential, 3)
                    if abs(potential[ir, iφ, iz + 1] - potential[ir, iφ, iz]) >= refinement_value_z
                        if axz[iz + 1] - axz[iz] > minimum_distance_z
                            if !in(iz, inds_z) push!(inds_z, iz) end
                        end
                    end
                end

            end
        end
    end

    if typeof(int2).parameters[2] == :open
        if isodd(length(inds_φ))
            for iφ in 1:size(potential, 2)
                if !in(iφ, inds_φ)
                    if iφ != size(potential, 2)
                        if axφ[iφ + 1] - axφ[iφ] > minimum_distance_φ
                            append!(inds_φ, iφ)
                            break
                        end
                    else
                        if cyclic - axφ[iφ] > minimum_distance_φ
                            append!(inds_φ, iφ)
                            break
                        end
                    end
                end
            end
        end
        if isodd(length(inds_φ))
            b::Bool = true
            diffs = diff(axφ)
            while b
                if length(diffs) > 0
                    idx = findmax(diffs)[2]
                    if !in(idx, inds_φ)
                        append!(inds_φ, idx)
                        b = false
                    else
                        deleteat!(diffs, idx)
                    end
                else
                    for i in eachindex(axφ)
                        if !in(i, inds_φ) push!(inds_φ, i) end
                    end
                    b = false
                end
            end
        end
        @assert iseven(length(inds_φ)) "Refinement would result in uneven grid in φ."
    end

    if isodd(length(inds_z))
        for iz in 1:size(potential, 3)
            if !in(iz, inds_z)
                if iz != size(potential, 3)
                    if axz[iz + 1] - axz[iz] > minimum_distance_z
                        append!(inds_z, iz)
                        break
                    end
                end
            end
        end
    end
    if isodd(length(inds_z))
        b = true
        diffs = diff(axz)
        while b
            if length(diffs) > 0
                idx = findmax(diffs)[2]
                if !in(idx, inds_z)
                    append!(inds_z, idx)
                    b = false
                else
                    deleteat!(diffs, idx)
                end
            else
                for i in eachindex(axz[1:end-1])
                    if !in(i, inds_z) push!(inds_z, i) end
                end
                b = false
            end
        end
    end
    if isodd(length(inds_z)) inds_z = inds_z[1:end-1] end
    @assert iseven(length(inds_z)) "Refinement would result in uneven grid in z."

    return sort(inds_r), sort(inds_φ), sort(inds_z)
end

function check_for_refinement( potential::Array{T, 3}, grid::Grid{T, 3, :cartesian}, max_diff::AbstractArray, minimum_distance::AbstractArray)::NTuple{3, Vector{Int}} where {T}
    inds_x::Vector{Int} = Int[]
    inds_y::Vector{Int} = Int[]
    inds_z::Vector{Int} = Int[]

    ax1::DiscreteAxis{T} = grid[1]
    ax2::DiscreteAxis{T} = grid[2]
    ax3::DiscreteAxis{T} = grid[3]
    int1::Interval = ax1.interval
    int2::Interval = ax2.interval
    int3::Interval = ax3.interval

    refinement_value_x::T = max_diff[1]
    refinement_value_y::T = max_diff[2]
    refinement_value_z::T = max_diff[3]

    minimum_distance_x::T = minimum_distance[1]
    minimum_distance_y::T = minimum_distance[2]
    minimum_distance_z::T = minimum_distance[3]

    ax_x::Vector{T} = collect(grid.axes[1])
    ax_y::Vector{T} = collect(grid.axes[2])
    ax_z::Vector{T} = collect(grid.axes[3])

    for iz in 1:size(potential, 3)
        z = ax_z[iz]
        for iy in 1:size(potential, 2)
            y = ax_y[iy]
            isum = iz + iy
            for ix in 1:size(potential, 1)
                x = ax_x[ix]

                # x
                if ix != size(potential, 1)
                    if abs(potential[ix + 1, iy, iz] - potential[ix, iy, iz]) >= refinement_value_x
                        if ax_x[ix + 1] - ax_x[ix] > minimum_distance_x
                            if !in(ix, inds_x) push!(inds_x, ix) end
                        end
                    end
                end

                # y
                if iy != size(potential, 2)
                    if abs(potential[ix, iy + 1, iz] - potential[ix, iy, iz]) >= refinement_value_y
                        if ax_y[iy + 1] - ax_y[iy] > minimum_distance_y
                            if !in(iy, inds_y) push!(inds_y, iy) end
                        end
                    end
                end

                # z
                if iz != size(potential, 3)
                    if abs(potential[ix, iy, iz + 1] - potential[ix, iy, iz]) >= refinement_value_z
                        if ax_z[iz + 1] - ax_z[iz] > minimum_distance_z
                            if !in(iz, inds_z) push!(inds_z, iz) end
                        end
                    end
                end

            end
        end
    end

    if isodd(length(inds_x))
        for ix in 1:size(potential, 1)
            if !in(ix, inds_x)
                if ix != size(potential, 1)
                    if ax_x[ix + 1] - ax_x[ix] > minimum_distance_x
                        append!(inds_x, ix)
                        break
                    end
                end
            end
        end
    end
    if isodd(length(inds_x))
        b = true
        diffs = diff(ax_x)
        while b
            if length(diffs) > 0
                idx = findmax(diffs)[2]
                if !in(idx, inds_x)
                    append!(inds_x, idx)
                    b = false
                else
                    deleteat!(diffs, idx)
                end
            else
                for i in eachindex(ax_x[1:end-1])
                    if !in(i, inds_x) push!(inds_x, i) end
                end
                b = false
            end
        end
    end
    if isodd(length(inds_x)) inds_x = inds_x[1:end-1] end
    @assert iseven(length(inds_x)) "Refinement would result in uneven grid in x. This is not allowed since this is the red black dimension."

    return sort(inds_x), sort(inds_y), sort(inds_z)
end


text:
using FlightSims
using Test
using Plots
using Transducers

function main()
    env = TwoDimensionalNonlinearOscillator()
    x10, x20 = 0.5, 0.5
    x0 = State(env)(x10, x20)
    tf = 10.0
    simulator = Simulator(x0,
                          apply_inputs(Dynamics!(env);
                          u=(x, p, t) -> FSimZoo.OptimalControl(env)(x));
                          tf=tf,
                         )
    df = solve(simulator)
    ts = df.time
    states = df.sol |> Map(datum -> datum.state) |> collect
    inputs = df.sol |> Map(datum -> datum.input) |> collect
    xs = states
    us = inputs
    fig_x = plot(ts, hcat(xs...)';
                 label=["x1" "x2"],
                 color=[:blue :red],
                )
    fig_u = plot(ts, hcat(us...)';
                 label="u",
                )
    fig = plot(fig_x, fig_u; layout=(2, 1))
    display(fig)
end

@testset "two_dimensional_nonlinear_oscillator" begin
    main()
end


text:
function maxModelLinear(::Type{T} ; resources::Dict, alternates = String[], allowMultiRecipes = true) where T <: Product

    allowedRecipes = union(baseRecipes, filter(r -> any(alt -> occursin(alt, r.name), alternates), allRecipes))

    m = Model(Cbc.Optimizer)
    @variable(m, 0 <= x[p in subtypes(Product)] <= 100_000_000)
    @variable(m, 0 <= y[r in allRecipes] <= 100_000_000)

    @objective(m, Max, x[T] - sum(qty * y[r] for (r, qty) in dependantRecipes(T))
        - 1e-6 * sum(y[r] for r in allowedRecipes)
        - 1e-6 * sum(y[r] for r in allRecipes if occursin("Alternate", r.name)))

    # Limit on recipes for all base resources (-> recipes with empty input)
    @constraint(m, [r in filter(r -> isempty(r.in), allRecipes)], y[r] * only(r.out)[2] <= get(resources, only(r.out)[1], 0))

    # Can't have more of a product than we have recipes producing it
    for p in subtypes(Product)
        if p in baseResources
            @constraint(m, x[p] <= sum(qty * y[r] for (r, qty) in recipes(p)))
        else
            @constraint(m, x[p] <= sum(qty * y[r] for (r, qty) in recipes(p)) + get(resources, p, 0))
        end
        if !allowMultiRecipes
            @constraint(m, [y[r] for (r, qty) in recipes(p)] in SOS1()) # only use 1 recipe to produce a specific product
        end
    end

    # For each Product, make sure we have enough of that product for all recipes depending on it
    @constraint(m, [p in subtypes(Product)], x[p] >= sum(qty * y[r] for (r, qty) in dependantRecipes(p)))

    @constraint(m, [r in allRecipes ; !(r in allowedRecipes)], y[r] == 0)

    return m, x, y
end

function maxModelMIP(::Type{T}, frac = 1/4 ; resources::Dict, alternates = String[], allowMultiRecipes = true, minProductionThreshold = 1) where T <: Product

    allowedRecipes = union(baseRecipes, filter(r -> any(alt -> occursin(alt, r.name), alternates), allRecipes))

    m = Model(Cbc.Optimizer)
    @variable(m, 0 <= x[p in subtypes(Product)] <= 100_000_000 / frac)
    @variable(m, 0 <= y[r in allRecipes] <= 100_000_000 / frac, Int)
    @variable(m, is_used[r in allRecipes], Bin)

    M = 100_000_000 / frac ; minProductionThreshold /= frac
    @constraint(m, [r in allRecipes], y[r] >= minProductionThreshold * is_used[r])
    @constraint(m, [r in allRecipes], y[r] <= M * is_used[r])

    for r in filter(r -> isempty(r.in), allowedRecipes)
        unset_integer(y[r])
    end

    @objective(m, Max, x[T] - sum(qty * y[r] for (r, qty) in dependantRecipes(T)) - 1e-6 * sum(r.duration * y[r] for r in allowedRecipes))

    # Limit on recipes for all base resources (-> recipes with empty input)
    @constraint(m, [r in filter(r -> isempty(r.in), allRecipes)], frac * y[r] * only(r.out)[2] <= get(resources, only(r.out)[1], 0))

    # Can't have more of a product than we have recipes producing it
    for p in subtypes(Product)
        if p in (Limestone, IronOre, CopperOre, CateriumOre, Coal, RawQuartz, Sulfur, Bauxite, SAMOre, Uranium, Water, CrudeOil)
            @constraint(m, x[p] <= sum(frac * qty * y[r] for (r, qty) in recipes(p)))
        else
            @constraint(m, x[p] <= sum(frac * qty * y[r] for (r, qty) in recipes(p)) + get(resources, p, 0))
        end
        if !allowMultiRecipes
            @constraint(m, [y[r] for (r, qty) in recipes(p)] in SOS1()) # only use 1 recipe to produce a specific product
        end
    end

    # For each Product, make sure we have enough of that product for all recipes depending on it
    @constraint(m, [p in subtypes(Product)], x[p] >= sum(frac * qty * y[r] for (r, qty) in dependantRecipes(p)))

    @constraint(m, [r in allRecipes ; !(r in allowedRecipes)], y[r] == 0)

    return m, x, y
end


text:
#!/usr/bin/env julia

using Luxor

function get_png_files(folder)
    cd(folder)
    imagelist = filter(f -> !startswith(f, ".") && endswith(f, "png"), readdir(folder))
    return map(realpath, imagelist)
end

function addimagetile(imgfile, xcenter, ycenter, tilewidth, tileheight; cropping=true)
    gsave()
    box(xcenter, ycenter, tilewidth, tileheight, :clip)
    img = readpng(imgfile)
    w = img.width
    h = img.height
    if cropping == true
        # expand image to fill tile
        if w < h
            scalefactor = max(tilewidth, tileheight)/w
        else
            scalefactor = max(tilewidth, tileheight)/h
        end
    else
        # shrink image to fit inside tile
        if w > h
            scalefactor = min(tilewidth, tileheight)/w
        else
            scalefactor = min(tilewidth, tileheight)/h
        end
    end
    translate(xcenter, ycenter)
    scale(scalefactor, scalefactor)
    placeimage(img, -w/2, -h/2)
    clipreset()
    grestore()
end

imagelist = get_png_files(dirname(@__FILE__) * "/../docs/figures")

shuffle!(imagelist)

width, height = 1600, 1000
fname = "/tmp/tiled-images.png"
Drawing(width, height, fname)
origin()
background("grey50")
setopacity(0.5)
pagetiles = Tiler(width, height, 8, 8, margin=50)
for (pos, n) in pagetiles
  if n > length(imagelist) # run out of images
    break
  end
  addimagetile(imagelist[n], pos.x, pos.y, pagetiles.tilewidth, pagetiles.tileheight, cropping=true)
end
finish()
println("finished test: output in $(fname)")


text:
using ReTest, Random, AdvancedHMC

@testset "Metric" begin
    @testset "Sample momentum variables from metric via vector of RNGs" begin
        D = 10
        n_chains = 5
        rng = [MersenneTwister(1) for _ in 1:n_chains]
        for metric in [
            UnitEuclideanMetric((D, n_chains)),
            DiagEuclideanMetric((D, n_chains)),
            # DenseEuclideanMetric((D, n_chains)) # not supported ATM
        ]
            r = rand(rng, metric)
            all_same = true
            for i in 2:n_chains
                all_same = all_same && r[:,i] == r[:,1]
            end
            @test all_same
        end
    end

    @testset "Resize metric" begin
        D = 10
        rng = MersenneTwister(1)
        θ = randn(rng, D)
        ℓπ(θ) = 1
        for metric in [
            UnitEuclideanMetric(1),
            DiagEuclideanMetric(1),
            DenseEuclideanMetric(1)
        ]
            h = Hamiltonian(metric, ℓπ, ℓπ)
            h = AdvancedHMC.resize(h, θ)
            @test size(h.metric) == size(θ)
        end
    end
end


text:
immutable MLPStorage{F, LayerDescriptor, N} <: Storage{MultilayerPerceptron{F, LayerDescriptor}, N}
  layer_buffer::Vector{F}
end

immutable MLPBackpropStorage{F, LayerDescriptor, N} <: BackpropStorage{MultilayerPerceptron{F, LayerDescriptor}, N}
  layer_buffer::Vector{F}
  backprop_buffer::Vector{F}
end

#implement allocations that allow creation of the default nonparametric allocations.
function (::Type{Storage}){F, LD}(::Type{MultilayerPerceptron{F, LD}}, N = :v)
  layercount = sum(LD[2:end-1]) * ((N == :v) ? 1 : N)
  MLPStorage{F, LD, N}(Vector{F}(layercount))
end
(::Type{Storage}){T <: MultilayerPerceptron}(::T, N = :v) = Storage(T, N)

function (::Type{BackpropStorage}){F, LD}(::Type{MultilayerPerceptron{F, LD}}, N = :v)
  layercount = sum(LD[2:end-1]) * ((N == :v) ? 1 : N)
  MLPBackpropStorage{F, LD, N}(Vector{F}(layercount), Vector{F}(layercount))
end
(::Type{BackpropStorage}){T <: MultilayerPerceptron}(::T, N = :v) = BackpropStorage(T)

doc"""
  `GenML.MLP.layer(::Storage{MultilayerPerceptron}, Val{N})`

  finds a reference to a layer in the execution storage for an MLP
"""
@generated function layer{F, LD, N, L}(s::Storage{MultilayerPerceptron{F, LD}, N}, ::Type{Val{L}})
  #retrieves the view on the array which represents layer n.  Note:
  #layer 1 is NOT allowed, since that should generically be the "input"
  #layer.
  (L < 2) && throw(ArgumentError("no storage for layers lower than 2"))
  #calculate the batch width
  n = (N == :v) ? 1 : N
  if (L == 2)
    start = 1
    finish = LD[2] * n
  else
    start = sum(LD[2:L-1]) * n + 1
    finish = sum(LD[2:L]) * n
  end

  if N == :v
    :(view(s.layer_buffer, $start:$finish))
  else
    layersize = LD[L]
    :(reshape(view(s.layer_buffer, $start:$finish), $layersize, N))
  end
end

doc"""
  `GenML.MLP.layer(::Storage{MultilayerPerceptron}, Val{N})`

  finds a reference to backpropagation data in the execution storage for an MLP
"""
@generated function backprop_layer{F, LD, N, L}(s::MLPBackpropStorage{F, LD, N}, ::Type{Val{L}})
  #retrieves the view on the array which represents layer n.  Note:
  #layer 1 is NOT allowed, since that should generically be the "input"
  #layer.
  (L < 2) && throw(ArgumentError("no storage for layers lower than 2"))
  n = (N == :v) ? 1 : N
  if (L == 2)
    start = 1
    finish = LD[2] * n
  else
    start = sum(LD[2:L-1]) * n + 1
    finish = sum(LD[2:L]) * n
  end

  if N == :v
    :(view(s.backprop_buffer, $start:$finish))
  else
    layersize = LD[L]
    :(reshape(view(s.backprop_buffer, $start:$finish), layersize, N))
  end
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
function compute_likelihood_save(time_phot,param,data,aex,al_small,indx,logdeta,p)

# First call prior:

#log_prior = compute_prior(param)
log_prior = 0.0

# Then call model:

model = model_tyc3559_v2(time_phot,param[1:16])

# Then compute log-likelihood:

# Subtract off the model from the data:
y = data - model

log_like = lorentz_likelihood_hermitian_band_save(p,y,aex,al_small,indx,logdeta)
# Apply prior on flux ratio:
log_like -= 0.5*(param[15]-0.692)^2/0.041^2

return log_like
end


text:
@testset "Crimisini" begin
  if visualtests
    # Diagonal
    img  = [i + j > 10 ? 0.1 : 0.2 for i in 1:10, j in 1:10]
    mask = [1.2cos(i) + i > j for i in 1:10, j in 1:10]
    out  = inpaint(img, mask, Criminisi(5,5))
    @plottest plot_before_after(img,mask,out) joinpath(datadir,"Diagonal.png") !istravis

    # Cross
    img = zeros(50,50) .+ 0.5
    for i in 1:50
      img[i,i] = 1.0
      img[i,50-i+1] = 1.0
    end
    mask = falses(size(img))
    mask[20:30,15:40] .= true
    out = inpaint(img, mask, Criminisi(9,9))
    @plottest plot_before_after(img,mask,out) joinpath(datadir,"Cross.png") !istravis

    # Blobs
    img  = Float64.(Gray.(testimage("blobs")))
    mask = falses(size(img))
    mask[50:150,50:150] .= true
    out = inpaint(img, mask, Criminisi(11,11))
    @plottest plot_before_after(img,mask,out) joinpath(datadir,"Blobs.png") !istravis

    # LightHouse
    img  = Float64.(Gray.(testimage("lighthouse")))
    mask = falses(size(img))
    mask[50:350,300:400] .= true
    out = inpaint(img, mask, Criminisi(30,30))
    @plottest plot_before_after(img,mask,out) joinpath(datadir,"LightHouse.png") !istravis
  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

# Collection of sources required to build LCIOBuilder
sources = [
    "https://github.com/iLCSoft/LCIO/archive/v02-12-01.tar.gz" =>
    "8a3d2e66c2f2d4489fc2e1c96335472728d913d4090327208a1d93b3b0728737",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mkdir build
cd build
ln -s /opt/x86_64-apple-darwin14/x86_64-apple-darwin14/sys-root/usr/include/sys /usr/include
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain ../LCIO-02-12-01/
VERBOSE=ON cmake --build . --config Release --target install
rm /usr/include/sys
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7,:cxx11)),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8,:cxx11)),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)),
]
# platforms = expand_gcc_versions(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
	"https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.3/build_Zlib.v1.2.11.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "LCIOBuilder", VersionNumber("v02.12.01-04"), sources, script, platforms, products, dependencies)


text:
import JSON

AnyIntType = Union{Int8, Int32, Int64, Int128, Int}

@enum BrainFlowMetrics begin

    RELAXATION = 0
    CONCENTRATION = 1

end

@enum BrainFlowClassifiers begin

    REGRESSION = 0
    KNN = 1

end

mutable struct BrainFlowModelParams

    metric::Int32
    classifier::Int32
    file::String
    other_info::String

    function BrainFlowModelParams(metric::AnyIntType, classifier::AnyIntType)
        new(Int32(metric), Int32(classifier), "", "")
    end

end

struct MLModel

    input_json::String

    function MLModel(params::BrainFlowModelParams)
        new(JSON.json(params))
    end

end

function prepare(ml_model::MLModel)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:prepare, "MLModule.dll"), Cint, (Ptr{UInt8},), ml_model.input_json)
    elseif Sys.isapple()
        ec = ccall((:prepare, "libMLModule.dylib"), Cint, (Ptr{UInt8},), ml_model.input_json)
    else
        ec = ccall((:prepare, "libMLModule.so"), Cint, (Ptr{UInt8},), ml_model.input_json)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in prepare ", ec), ec))
    end
end

function release(ml_model::MLModel)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:release, "MLModule.dll"), Cint, (Ptr{UInt8},), ml_model.input_json)
    elseif Sys.isapple()
        ec = ccall((:release, "libMLModule.dylib"), Cint, (Ptr{UInt8},), ml_model.input_json)
    else
        ec = ccall((:release, "libMLModule.so"), Cint, (Ptr{UInt8},), ml_model.input_json)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in release ", ec), ec))
    end
end

function predict(data, ml_model::MLModel)
    val = Vector{Float64}(undef, 1)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:predict, "MLModule.dll"), Cint, (Ptr{Float64}, Cint, Ptr{Float64}, Ptr{UInt8}),
            data, length(data), val, ml_model.input_json)
    elseif Sys.isapple()
        ec = ccall((:predict, "libMLModule.dylib"), Cint, (Ptr{Float64}, Cint, Ptr{Float64}, Ptr{UInt8}),
            data, length(data), val, ml_model.input_json)
    else
        ec = ccall((:predict, "libMLModule.so"), Cint, (Ptr{Float64}, Cint, Ptr{Float64}, Ptr{UInt8}),
            data, length(data), val, ml_model.input_json)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in predict ", ec), ec))
    end
    value = val[1]
    value
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

function gen_rtic_grid!{S<:Union{Seg2,Seg3}}(thetas, element::Element{S}; npts=2)
    xi1 = midpoints(linspace(-1, 1, npts+1))
    for xi in xi1
        push!(thetas, [0.0, xi])
    end
end

function gen_rtic_grid!(thetas, element::Element{NSeg}; npts=2)
    knots_u = element.properties.knots
    u = midpoints(linspace(minimum(knots_u), maximum(knots_u), npts+1))
    for ui in u
        push!(thetas, [0.0, ui])
    end
end

function gen_rtic_grid!(thetas, element::Element{NSurf}; npts=2)
    knots_u = element.properties.knots_u
    knots_v = element.properties.knots_v
    u = midpoints(linspace(minimum(knots_u), maximum(knots_u), npts+1))
    v = midpoints(linspace(minimum(knots_v), maximum(knots_v), npts+1))
    for ui in u
        for vi in v
            push!(thetas, [0.0, ui, vi])
        end
    end
end

"""
Find intersection between element surface X(ξ) and ray x(t) = s + t*n by
solving equation F(t, ξ) = x(t) - X(ξ) = 0

Returns
-------
t, ξ

References
----------
[1] https://en.wikipedia.org/wiki/Ray_tracing_%28graphics%29

"""
function find_intersection{S<:Union{Seg2, Seg3, NSeg, NSurf}}(element::Element{S}, s, n, time;
    max_iterations=10, tolerance=1.0e-12, info_output=false, deformed=true, secondary=false, npts=2)

    x = element["geometry"](time)
    if deformed && haskey(element, "displacement")
        x += element["displacement"](time)
    end

    function calc_intersection!(theta)
        i = 0
        dtheta = zeros(theta)
        for i=1:max_iterations
            t = theta[1]
            xi = theta[2:end]
            N = get_basis(element, xi, time)
            dN = get_dbasis(element, xi, time)
            A = [n -dN*x]
            r = s + t*n - N*x
            try
                dtheta = A \ -r
            catch err
                info("failed to solve equation.")
                dump(theta)
                dump(A)
                dump(r)
                throw(err)
            end
            theta[:] += dtheta
            info_output && info("iter $i, theta = $theta, norm(dtheta) = $(norm(dtheta))")
            norm(dtheta) < tolerance && return theta
            isnan(theta[1]) && break
        end
        theta[1] = NaN
        return theta
        # error("didn't converge in $i iterations")
    end

    thetas = Vector{Float64}[]
    gen_rtic_grid!(thetas, element; npts=npts)
    map(calc_intersection!, thetas)
    filter!(t -> !isnan(t[1]), thetas)
    info_output && info("thetas: $thetas")
    secondary && filter!(t -> t[1] > 1.0e-12, thetas)
    # error("failure finding ray trace, thetas vec is empty")
    length(thetas) == 0 && return NaN, [NaN, NaN]
    sort!(thetas, alg=MergeSort, lt=(a,b)->a[1]<b[1])
    theta = thetas[1]
    return theta[1], theta[2:end]
end

function calc_normal{S<:Union{Seg2, Seg3, NSeg}}(element::Element{S}, xi, time; deformed=true)
    x = element["geometry"](time)
    if deformed && haskey(element, "displacement")
        x += element["displacement"](time)
    end
    Q = [0.0 1.0; -1.0 0.0]
    dN = get_dbasis(element, xi, time)
    n = Q*(dN*x)
    n /= norm(n)
    return n
end

function calc_normal{S<:Union{Quad4, NSurf}}(element::Element{S}, xi, time; deformed=true)
    x = element["geometry"](time)
    if deformed && haskey(element, "displacement")
        x += element["displacement"](time)
    end
    dN = get_dbasis(element, xi, time)
    J = transpose(sum([kron(dN[:,i], x[i]') for i=1:length(x)]))
    n = cross(J[:,1], J[:,2])
    n /= norm(n)
    return n
end

"""
References
----------
[1] http://fp.optics.arizona.edu/optomech/Fall13/Notes/6%20Mirror%20matrices.pdf
"""
function calc_reflection(element::Element, xi, k, time; deformed=true)
    n = calc_normal(element, xi, time; deformed=deformed)
    k2 = k - 2*vecdot(k, n)*n
    return k2
end


text:
import ..Codim1LimitCycle: limitcycles
using ..Codim1LimitCycle: AbstractLimitCycleData

function dim_ds_state(sweep::Codim2LCSweep)
    prob = as(sweep, ContinuationSweep).sol.value.prob
    return size(prob.xs0, 1)
end

struct FoldLimitCycleData{M, R, V} <: AbstractLimitCycleData
    state::M
    period::R
    param1_value::V
    param2_value::V
    # param1_axis::L
    # param2_axis::L
end

function FoldLimitCycleData(u::AbstractVector, n::Integer)
    state_elements, r = divrem(length(u) - 4, 2)
    @assert r == 0
    FoldLimitCycleData(
        reshape((@view u[1:state_elements]), n, :),
        u[state_elements + 1],  # period
        u[end - 1],             # param1_value
        u[end],                 # param2_value
    )
end
# See: [[./problem.jl::get_u0.*::]]

function limitcycles(sweep::Codim2LCSweep)
    n = dim_ds_state(sweep)
    u_list = as(sweep, ContinuationSweep).u
    return FoldLimitCycleData.(u_list, (n,))
end

limitcycles(sol::Codim2LCSolution) = vcat(limitcycles.(sol.sweeps)...)
limitcycles(solver::Codim2LCSolver) = limitcycles(solver.sol)

BifurcationsBase.measure(lc::FoldLimitCycleData, i::Integer) =
    @view lc.state[i, :]

# TODO: make it more type-based
function BifurcationsBase.measure(lc::FoldLimitCycleData, key::Symbol)
    if key == :p1
        return lc.param1_value
    elseif key == :p2
        return lc.param2_value
    elseif key == :period
        return lc.period
    end
    error("Unsupported key: $key")
end

# TODO: make it more type-based
function BifurcationsBase.measure(sweep::Codim2LCSweep, key::Symbol)
    u_list = as(sweep, ContinuationSweep).u
    if isempty(u_list)
        state_elements = 0  # won't be used
    else
        state_elements, r = divrem(length(u_list[1]) - 4, 2)
        @assert r == 0
    end
    if key == :p1
        return [u[end - 1] for u in u_list]
    elseif key == :p2
        return [u[end] for u in u_list]
    elseif key == :period
        return [u[state_elements + 1] for u in u_list]
    end
    error("Unsupported key: $key")
end


text:
module Qualification
using MLStyle.MatchCore
using MLStyle.Pervasives

export get_qualifier
get_qualifier(node, curmod) =
    @match node begin
        :public   => invasive
        :internal => internal
        :(visible in [$(mods...)]) ||
        :(visible in $mod) && Do(mods = [mod]) =>
            share_with(Set(map(curmod.eval, mods)))
    end
end


text:
export DaxKaczmarz

mutable struct DaxKaczmarz{matT,T,U,Tsparse} <: AbstractLinearSolver
  S::matT
  u::Vector{T}
  λ::Float64
  denom::Vector{U}
  rowindex::Vector{Int64}
  sumrowweights::Vector{Float64}
  zk::Vector{T}
  bk::Vector{T}
  xl::Vector{T}
  yl::Vector{T}
  εw::Vector{T}
  τl::T
  αl::T
  weights::Vector{U}
  enforceReal::Bool
  enforcePositive::Bool
  sparseTrafo::Tsparse
  iterations::Int64
  iterationsInner::Int64
end

"""This function solves a unconstrained linear least squaares problem using an algorithm proposed in [1] combined with a randomized version of kaczmarz [2].
Returns an approximate solution to the linear leaast squares problem Sᵀx = u.

[1] Dax, A. On Row Relaxation Methods for Large Constrained Least Squares Problems. SIAM J. Sci. Comput. 14, 570–584 (1993).
[2] Strohmer, T. & Vershynin, R. A Randomized Kaczmarz Algorithm with Exponential Convergence. J. Fourier Anal. Appl. 15, 262–278 (2008).

### Input Arguments
* `S::AbstractMatrix{T}`: Problem Matrix S.
* `u::Vector{T}`: Righthandside of the linear equation.

### Keyword/Optional Arguments

* `iterations::Int`: Number of Iterations of outer dax scheme.
* `iterationsInner::Int`: Number of Iterations of inner dax scheme.
* `λ::Float64`: The regularization parameter ɛ>0 influences the speed of convergence but not the solution.
* `weights::Vector{T}`: Use weights in vector to weight equations. The larger the weight the more one 'trusts' a sqecific equation.
"""
function DaxKaczmarz(S, b=nothing; λ::Real=0.0
              , weights::Vector{R}=ones(Float64,size(S,1))
              , sparseTrafo=nothing
              , enforceReal::Bool=false
              , enforcePositive::Bool=false
              , iterations::Int=3
              , iterationsInner::Int=2
              , kargs...) where R <: Real

  # setup denom and rowindex
  sumrowweights, denom, rowindex = initkaczmarzdax(S,λ,weights)

  T = typeof(real(S[1]))
  M,N = size(S)
  if b != nothing
    u = b
  else
    u = zeros(eltype(S),M)
  end
  zk = zeros(eltype(S),N)
  bk = zeros(eltype(S),M)
  xl = zeros(eltype(S),N)
  yl = zeros(eltype(S),M)
  εw = zeros(eltype(S),length(rowindex))
  τl = zero(eltype(S))
  αl = zero(eltype(S))

  return DaxKaczmarz(S,u,Float64(λ),denom,rowindex,sumrowweights,zk,bk,xl,yl,εw,τl,αl
                  ,T.(weights),enforceReal,enforcePositive
                  ,sparseTrafo,iterations,iterationsInner)
end

function init!(solver::DaxKaczmarz
              ; S::matT=solver.S
              , λ::Real=solver.λ
              , u::Vector{T}=eltype(S)[]
              , zk::Vector{T}=eltype(S)[]
              , weights::Vector{Float64}=solver.weights) where {matT,T}

  if S != solver.S
    solver.sumrowweights, solver.denom, solver.rowindex = initkaczmarzdax(S,solver.λ,solver.weights)
  end
  solver.λ = Float64(λ)

  solver.u[:] .= u
  solver.weights=weights

  # start vector
  if isempty(zk)
    solver.zk[:] .= zeros(T,size(S,2))
  else
    solver.zk[:] .= x
  end

  solver.bk[:] .= zero(T)
  solver.xl[:] .= zero(T)
  solver.yl[:] .= zero(T)
  solver.αl = zero(T)        #temporary scalar
  solver.τl = zero(T)        #temporary scalar

  for i=1:length(solver.rowindex)
    j = solver.rowindex[i]
    solver.ɛw[i] = sqrt(solver.λ)/weights[j]
  end
end

function solve(solver::DaxKaczmarz, u::Vector{T}; λ::Real=solver.λ
                , S::matT=solver.S, startVector::Vector{T}=eltype(S)[]
                , weights::Vector=solver.weights
                , solverInfo=nothing, kargs...) where {T,matT}

  # initialize solver parameters
  init!(solver; S=S, λ=λ, u=u, zk=startVector, weights=weights)

  # log solver information
  solverInfo != nothing && storeInfo(solverInfo,solver.zk,norm(bk))

  # perform CGNR iterations
  for (iteration, item) = enumerate(solver)
    solverInfo != nothing && storeInfo(solverInfo,solver.zk,norm(bk))
  end

  return solver.zk
end

function iterate(solver::DaxKaczmarz, iteration::Int=0)
  if done(solver,iteration)
    applyConstraints(solver.zk, solver.sparseTrafo, solver.enforceReal, solver.enforcePositive)
    return nothing
  end

  copyto!(solver.bk, solver.u)
  gemv!('N',-1.0,solver.S,solver.zk,1.0,solver.bk)

  # solve min ɛ|x|²+|W*A*x-W*bk|² with weightingmatrix W=diag(wᵢ), i=1,...,M.
  for l=1:length(solver.rowindex)*solver.iterationsInner
    i::Int64 = getrandindex(solver.sumrowweights,rand()*solver.sumrowweights[end])  #choose row with propability proportional to its weight.
    j = solver.rowindex[i]
    solver.τl = dot_with_matrix_row(solver.S,solver.xl,j)
    solver.αl = solver.denom[i]*(solver.bk[j]-solver.τl-solver.ɛw[i]*solver.yl[j])
    kaczmarz_update!(solver.S,solver.xl,j,solver.αl)
    solver.yl[j] += solver.αl*solver.ɛw[i]
  end

  BLAS.axpy!(1.0,solver.xl,solver.zk)  # zk += xl
  # reset xl and yl for next Kaczmarz run
  rmul!(solver.xl,0.0)
  rmul!(solver.yl,0.0)

  return solver.bk, iteration+1

end

@inline done(solver::DaxKaczmarz,iteration::Int) = iteration>=solver.iterations

"""This funtion saves the denominators to compute αl in denom and the rowindices,
  which lead to an update of cl in rowindex.
"""
function initkaczmarzdax(S::AbstractMatrix, ɛ, weights::Vector)
  length(weights)==size(S,1) ? nothing : error("number of weights must equal number of equations")
  denom = Float64[]
  sumrowweights = Float64[]
  rowindex = Int64[]

  push!(sumrowweights,0.0)
  for i=1:size(S,1)
    s² = rownorm²(S,i)*weights[i]^2
    if s²>0
      push!(denom,weights[i]^2/(s²+ɛ))
      push!(sumrowweights,s²+ɛ+sumrowweights[end])
      push!(rowindex,i)
    end
  end
  return sumrowweights, denom, rowindex
end

"This function returns the next index for the randomized Kaczmarz algorithm."
function getrandindex(v::AbstractVector, x)
    lo::Int64 = 0
    hi::Int64 = length(v)+1
    while lo < hi-1
        m = (lo+hi)>>>1
        x<v[m] ? hi = m : lo = m
    end
    return lo
end


text:
function get_vecs!((phi,q),M,V,AV,ni)
  F = eigen(Hermitian(M))
  lambda = F.values[1]
  u = F.vectors[:,1]
  phi .= u[1] .* V[1]
  q .= u[1] .* AV[1]
  for n=2:ni
    phi .+= u[n] .* V[n]
    q .+= u[n] .* AV[n]
  end
  q .-= lambda .* phi
  #Fix sign
  if real(u[1]) < 0.0
    phi .*= -1
    q .*= -1
  end
  return lambda
end

function orthogonalize!(q::ITensor,V,ni)
  q0 = copy(q)
  for k=1:ni
    Vq0k = dot(V[k],q0)
    q .+= -Vq0k .* V[k]
  end
  qnrm = norm(q)
  if qnrm < 1E-10 #orthog failure, try randomizing
    randn!(q)
    qnrm = norm(q)
  end
  q .*= 1.0/qnrm
  return
end

function expand_krylov_space(M::Matrix{elT},V,AV,ni) where {elT}
  newM = fill(zero(elT),(ni+1,ni+1))
  newM[1:ni,1:ni] = M
  for k=1:ni+1
    newM[k,ni+1] = dot(V[k],AV[ni+1])
    newM[ni+1,k] = conj(newM[k,ni+1])
  end
  return newM
end

function davidson(A,
                  phi0::ITensorT;
                  kwargs...) where {ITensorT<:ITensor}
  elTA   = eltype(A)
  elTphi = eltype(phi0)

  # if the matrix is complex and the starting vector is real,
  # that's not going to last long
  # and we're going to need phi to be complex
  # down when we reuse the storage in get_vecs!
  if !(elTA <: Real) && (elTphi <: Real)
    phi = complex(phi0)
  else
    phi = copy(phi0)
  end

  maxiter = get(kwargs,:maxiter,2)
  miniter = get(kwargs,:miniter,1)
  errgoal = get(kwargs,:errgoal,1E-14)
  Northo_pass = get(kwargs,:Northo_pass,1)

  approx0 = 1E-12

  nrm = norm(phi)
  if nrm < 1E-18
    phi_ = similar(phi)
    randn!(phi_)
    phi = phi_
    nrm = norm(phi)
  end
  phi .*= 1.0/nrm

  maxsize = size(A)[1]
  actual_maxiter = min(maxiter,maxsize-1)

  if dim(inds(phi)) != maxsize
    error("linear size of A and dimension of phi should match in davidson")
  end

  Aphi = A(phi)

  V = ITensorT[copy(phi)]
  AV = ITensorT[Aphi]

  last_lambda = NaN
  lambda::Float64 = real(dot(V[1],AV[1]))
  q = AV[1] - lambda*V[1];

  M = fill(elTA(lambda),(1,1))

  for ni=1:actual_maxiter

    qnorm = norm(q)

    errgoal_reached = (qnorm < errgoal && abs(lambda-last_lambda) < errgoal)
    small_qnorm = (qnorm < max(approx0,errgoal*1E-3))
    converged = errgoal_reached || small_qnorm

    if (qnorm < 1E-20) || (converged && ni > miniter) #|| (ni >= actual_maxiter)
      #@printf "  done with davidson, ni=%d, qnorm=%.3E\n" ni qnorm
      break
    end

    last_lambda = lambda

    for pass = 1:Northo_pass
      orthogonalize!(q,V,ni)
    end

    Aq = A(q)

    push!(V,copy(q))
    push!(AV,Aq)

    M = expand_krylov_space(M,V,AV,ni)

    lambda = get_vecs!((phi,q),M,V,AV,ni+1)

  end #for ni=1:actual_maxiter+1

  return lambda,phi

end


text:
helloworld = UInt8[0x68, 0x65, 0x6C, 0x6C, 0x6F, 0x20, 0x77, 0x6F, 0x72, 0x6C, 0x64, 0x00]
@test VirtualTestDrive.string_to_zeroed_null_terminated_array("") == UInt8[0x00]
@test VirtualTestDrive.string_to_zeroed_null_terminated_array("hello world") == helloworld
@test VirtualTestDrive.string_to_zeroed_null_terminated_array("", 2) == UInt8[0x00, 0x00]
@test VirtualTestDrive.string_to_zeroed_null_terminated_array("hello world", 4) == UInt8[0x68, 0x65, 0x6C, 0x00]
@test VirtualTestDrive.string_to_zeroed_null_terminated_array("hello world", 14) == append!(helloworld, zeros(UInt8, 2))

@test VirtualTestDrive.string_from_buffer(helloworld) == "hello world"
@test VirtualTestDrive.string_from_buffer(
    VirtualTestDrive.string_to_zeroed_null_terminated_array("hi",10)
    ) == "hi"

e = ETree("Set")
e.attr["entity"] = "player"
push!(e.elements, ETree("Speed"))
e.elements[1].attr["value"] = "29.06"

f = VirtualTestDrive.message_payload_to_etree(
    "<Set entity=\"player\"><Speed value=\"29.06\"/></Set>")

@test f.name == e.name
@test f.attr == e.attr
@test f.elements[1].name == e.elements[1].name
@test f.elements[1].attr == e.elements[1].attr
@test length(f.elements) == length(e.elements)
@test f == e

@test VirtualTestDrive.message_contains_packet_with_element(
    e, "Set", "Speed")

@test isa(find_vtd_directory(false), String)

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
using WaterLily
using LinearAlgebra: norm2
using JLD2

function circle(radius=8;Re=250,n=10,m=6)
    center, ν = radius*m/2, radius/Re
    body = AutoBody((x,t)->norm2(x .- center) - radius)
    Simulation((n*radius+2,m*radius+2), [1.,0.], radius; ν, body)
end

function sim_integrate!(sim;duration=1,step=0.1,verbose=false,
                        remeasure=false,kv...)
    t₀ = round(sim_time(sim))
    t = range(t₀,t₀+duration;step)
    duration = @elapsed for tᵢ in t
    sim_step!(sim,tᵢ;remeasure)
    verbose && println("tU/L=",round(tᵢ,digits=4),
    ", Δt=",round(sim.flow.Δt[end],digits=3))
    end
    return duration
    end

sim = circle(16)
runtime = sim_integrate!(sim;duration=1,step=0.25)
filename_u = string(tempdir(), "/testing_using_threads#",
                    Threads.nthreads(), "_u.jld2")
save(filename_u, "data", sim.flow.u)
# filename_t = string(tempdir(), "\\testing_using_threads#",
#                     Threads.nthreads(), "_t.jld")
# save(filename_t, "data", runtime)
@info("Saved final velocity field to ", filename_u)


text:
import PTChain
import PlotChain

function score_convex(vec)
    return sum(abs.(vec))
end

function min_convex(micro_steps, macro_steps)
    dims = 100
    start = randn(dims)*100
    temps = [128, 64, 32, 16, 8, 4, 2, 1, 0.5, 0.25, 0.125, 0]
    min_val, x = PTChain.sampleTempering(start, score_convex, temps, widescaleRandomNoiseTL, micro_steps, macro_steps)
    return min_val, x
end

function min_convex_SA()
    dims = 100
    start = randn(dims)*100
    reps = 15000
    temps = [128, 64, 32, 16, 8, 4, 2, 1, 0.5, 0.25, 0.125, 0]
    min_val = score_vec(start)
    min_x = copy(start)
    for i = 1
        val, x = PTChain.simulatedAnnealing(start, score_convex, temps, widescaleRandomNoiseTL, fill(reps, length(temps)))
        if min_val > val
            min_val = copy(val)
            min_x = copy(x)
        end
    end
    min_val, min_x
end

function compare_convex(micro_steps, macro_steps)
    dims = 100
    start = randn(dims)*100
    temps = collect(128*(0.85.^(0:40)))
    min_val_st, x_st, p_st = PTChain.sampleTempering(start, score_convex, temps, widescaleRandomNoiseTL, micro_steps, macro_steps)
    min_val_gt, x_gt, p_gt = PTChain.greedyTempering(start, score_convex, reverse(temps), widescaleRandomNoiseTL, micro_steps, macro_steps)
    min_val_sa, x_sa = PTChain.simulatedAnnealing(start, score_convex, temps, widescaleRandomNoiseTL, fill(micro_steps*macro_steps, length(temps)))
    p_g = plot(p_gt', title="Greedy Tempering")
    p_s = plot(p_st', title="Parallel Tempering")
    plot(p_g, p_s, layout=(2,1))
    gui()
    println("Parallel Tempering: $(min_val_st)")
    println("Greedy Tempering: $(min_val_gt)")
    println("Simulated Annealing: $(min_val_sa)")
end

function score_vec(vec)
    score = 0
    for i = 1:length(vec)
        if vec[i] == 0
            score += 2
        elseif vec[i] == -1
            score += 1
        elseif vec[i] == 1
            score += 0
        elseif abs(vec[i]) <= 4
            score += 2*abs(vec[i]) - 1
        elseif abs(vec[i]) <= 6
            score += 15 - 2*abs(vec[i])
        else
            score += abs(vec[i]) - 4
        end
    end
    score
end

function random_walk(x::Vector)
    n = length(x)
    direction = rand(1:n)
    y = copy(x)
    y[direction] = y[direction] + rand(-1:2:1)
    return y
end

function widescaleRandomNoiseTL(x::Vector{Float64})
    n = length(x)
    expon = rand()*6 - 4
    displacement = randn(n)
    new_x = x + displacement*10^expon
end

function annealSolverTL(dims)
    start = fill(-7, dims)
    reps = 100
    temps = 7*(0.975.^(0:99))
    min_val = score_vec(start)
    min_x = copy(start)
    for i = 1
        val, x = PTChain.simulatedAnnealing(start, score_vec, temps, random_walk, fill(reps, length(temps)))
        if min_val > val
            min_val = copy(val)
            min_x = copy(x)
        end
    end
    min_val, min_x
end

function insSolver(dims)
    start = rand(-10:10, dims)
    reps = 100
    temps = 7*(0.975.^(0:99))
    min_val, x = PTChain.sampleTempering(start, score_vec, [0, 0.5, 0.8, 5, 7], random_walk, reps, 100)
    return min_val, x
end

function compare_score_vec(dims,micro_steps, macro_steps, reps)
    min_val_st = zeros(reps)
    min_val_gt = zeros(reps)
    min_val_sa = zeros(reps)
    for i = 1:reps
        start = fill(-7, dims)
        temps = 15 * 0.7.^(0:15)
        min_val_st[i], _, _ = PTChain.sampleTempering(start, score_vec, temps, random_walk, micro_steps, macro_steps)
        min_val_gt[i], _, _ = PTChain.greedyTempering(start, score_vec, temps, random_walk, micro_steps, macro_steps)
        min_val_sa[i], _ = PTChain.simulatedAnnealing(start, score_vec, temps, random_walk, fill(micro_steps*macro_steps, length(temps)))
    end
    println("Parallel Tempering: $(min_val_st), mean: $(mean(min_val_st))")
    println("Greedy Tempering: $(min_val_gt), mean: $(mean(min_val_gt))")
    println("Simulated Annealing: $(min_val_sa), mean: $(mean(min_val_sa))")
end

function plot_run(dims, micro_steps, macro_steps)
    start = fill(-7, dims)
    temps = collect(7*(0.65.^(0:9)))
    min_val_st, _, plot_data_st = PTChain.sampleTempering(start, score_vec, temps, random_walk, micro_steps, macro_steps)
    min_val_gt, _, plot_data_gt = PTChain.greedyTempering(start, score_vec, temps, random_walk, micro_steps, macro_steps)
    #p_gt = plot(plot_data_gt', title="Greedy Tempering")
    #p_st = plot(plot_data_st', title="Parallel Tempering")
    #plot(p_gt, p_st, layout=(2,1))
    #gui()
    min_val_gt, min_val_st
end

#calculates the expected number of steps from a position on the function to a max
#steps up the recurrance that F(x) = p*(F(x+1) + 1) + (1 - p)*(F(x - 1) + 1)
function rand_walk_solver(steps, temp, ΔF=1)
    p = exp(-ΔF/temp)/2
    vec = ones(steps)
    vec[end] = 0
    mat = zeros(steps, steps)
    mat[1, 1] = p
    mat[1, 2] = -p
    for  i = 2:(steps - 1)
        mat[i, i - 1] = -(1 - p)
        mat[i, i] = 1
        mat[i, i + 1] = -p
    end
    mat[steps, steps] = 1
    return inv(mat)*vec
end

compare_score_vec(20, 500, 100, 1)

max_stars_repo_name: jacob-ruth/pt-optimization

max_issues_repo_name: jacob-ruth/pt-optimization

max_forks_repo_name: jacob-ruth/pt-optimization


text:
const TupleType = Type # should be Tuple{T1,T2,...}
const Req = Tuple{Function, TupleType}

abstract type AbstractRequirementSet end

mutable struct Unspecified <: AbstractRequirementSet
    requirer
    parent::Union{Nothing, Any}
end

Unspecified(requirer) = Unspecified(requirer, nothing)

mutable struct RequirementSet <: AbstractRequirementSet
    requirer
    reqs::Vector{Req} # not actually a set - to preserve intuitive ordering
    deps::Vector{AbstractRequirementSet}
    parent::Union{Nothing, Any}
    exception::Union{Nothing, Exception}
end

function RequirementSet(requirer, parent=nothing)
    return RequirementSet(requirer,
                          Vector{Tuple{Function, TupleType}}(),
                          AbstractRequirementSet[],
                          parent,
                          nothing)
end

Base.push!(r::RequirementSet, func::Function, argtypes::TupleType) = push!(r, (func, argtypes))
Base.push!(r::RequirementSet, t::Tuple{Function, TupleType}) = push!(r.reqs, t)

function push_dep!(r::RequirementSet, dep::AbstractRequirementSet)
    dep.parent = r.requirer
    push!(r.deps, dep)
end

"""
Return an expression that creates a RequirementSet using the code in the block. The resulting code will *always* return a RequirementSet, but it may be incomplete if the exception field is not null.
"""
function pomdp_requirements(name::Union{Expr,Symbol,String}, block::Expr)
    block = deepcopy(block)
    req_found = handle_reqs!(block, :reqs)
    if !req_found
        block = esc(block)
        @warn("No @req or @subreq found in @POMDP_requirements block.")
    end

    newblock = quote
        reqs = RequirementSet($(esc(name)))
        try
            $block
        catch exception
            reqs.exception = exception
        end
        reqs
    end
    return newblock
end

const CheckedList = Vector{Tuple{Union{Bool,Missing}, Function, TupleType}}

"""
Return a `(f, Tuple{T1,T2})::Req` expression given a `f( ::T1, ::T2)` expression.
"""
function convert_req(ex::Expr)
    malformed = false
    if ex.head == :call
        func = ex.args[1]
        argtypes = Union{Symbol, Expr}[]
        for a in ex.args[2:end]
            if isa(a, Expr)
                if a.head == :(::)
                    if length(a.args) == 1
                        push!(argtypes, a.args[1])
                    elseif length(a.args) == 2
                        push!(argtypes, a.args[2])
                    else
                        malformed = true
                        break
                    end
                else
                    malformed = true
                    break
                end
            else
                push!(argtypes, :(typeof($a)))
            end
        end
    else
        malformed = true
    end
    if malformed # throw error at parse time so solver writers will have to deal with this
        error("""
              Malformed requirement expression: $ex
              Requirements should be expressed in the form `function_name(::Type1, ::Type2)` or `function_name(arg1, arg2)`.
              """)
    else
        return quote ($func, Tuple{$(argtypes...)}) end
    end
end

function recursively_show(io::IO,
                           r::RequirementSet,
                           analyzed::Set,
                           reported::Set{Req})
    if r.requirer in analyzed
        return true
    end

    push!(analyzed, r.requirer)

    checked = CheckedList()
    allthere = true
    for fp in r.reqs
        if !(fp in reported)
            push!(reported, fp)
            exists = implemented(first(fp), last(fp))
            allthere = exists & allthere
            push!(checked, (exists, first(fp), last(fp)))
        end
    end

    show_requirer(io, r)
    if isempty(checked)
        println(io, "  [No additional requirements]")
    else
        show_checked_list(io, checked)
    end

    if r.exception == nothing # no exception
        first_exception = nothing
    else
        allthere = false
        show_incomplete(io, r)
        first_exception = r
    end

    for dep in r.deps
        depcomplete, depexception = recursively_show(io, dep, analyzed, reported)
        allthere = allthere & depcomplete
        if first_exception == nothing && depexception != nothing
            first_exception = depexception
        end
    end

    return allthere, first_exception
end

function recursively_show(io::IO, r::Unspecified, analyzed::Set, reported::Set{Req})
    if r.requirer in analyzed
        return true, nothing
    else
        push!(analyzed, r.requirer)
        show_requirer(io::IO, r)
        println(io, "  [No requirements specified]")
        return true, nothing
    end
end

function recursively_check(r::RequirementSet, analyzed::Set)
    if r.requirer in analyzed
        return true
    end

    push!(analyzed, r.requirer)

    allthere = r.exception == nothing
    for fp in r.reqs
        allthere = allthere & implemented(first(fp), last(fp))
    end

    for dep in r.deps
        allthere = allthere & recursively_check(dep, analyzed)
    end

    return allthere
end

function recursively_check(r::Unspecified, analyzed::Set)
    push!(analyzed, r.requirer)
    return true
end

"""
Return a tuple (not an Expr) of the function name, arguments, and argument types.

E.g. `f(arg1::T1, arg2::T2)` would be unpacked to (:f, [:arg1, :arg2], [:T1, :T2])
"""
function unpack_typedcall(typedcall::Expr)
    malformed = false
    if typedcall.head != :call
        malformed = true
    end

    args = Union{Symbol,Expr}[]
    types = Union{Symbol,Expr}[]
    for expr in typedcall.args[2:end]
        if isa(expr,Expr) && expr.head == :(::)
            push!(args, expr.args[1])
            push!(types, expr.args[2])
        elseif isa(expr,Symbol)
            push!(args, expr)
            push!(types, :Any)
        else
            malformed = true
        end
    end

    if malformed
        error("""
              Malformed typed funciton call expression: $typedcall
              Expected the form `function_name(arg1::Type1, arg2::Type2)`.
              """)
    end

    return (typedcall.args[1], args, types)
end

"""
Return a `(f, (arg1,arg2))` expression given a `f(arg1, arg2)` expression.
"""
function convert_call(call::Expr)
    malformed = false
    if call.head == :call
        func = call.args[1]
        args = Union{Symbol, Expr}[]
        for a in call.args[2:end]
            if isa(a, Expr) && a.head == :(::)
                @assert length(args) == 2
                push!(args, a.args[1])
            else
                push!(args, a)
            end
        end
    else
        malformed = true
    end
    if malformed # throw error at parse time so solver writers will have to deal with this
        error("""
              Malformed call expression: $call
              Expected the form `funcion_name(arg1, arg2)`
              """)
    else
        return quote ($func, ($(args...),)) end
    end
end

# this is where the freaking magic happens.
"""
    handle_reqs!(block, reqs_name::Symbol)

Replace any @req calls with `push!(\$reqs_name, <requirement>)`

Returns true if there was a requirement in there and so should not be escaped.
"""
function handle_reqs!(node::Expr, reqs_name::Symbol)

    if node.head == :macrocall && node.args[1] == Symbol("@req")
        macro_node = copy(node)
        node.head = :call
        expanded = macroexpand(POMDPLinter, macro_node)
        if isa(expanded, Expr) && expanded.head == :error
            throw(expanded.args[1])
        end
        node.args = [:push!, reqs_name, esc(expanded)]
        return true
    elseif node.head == :macrocall && node.args[1] == Symbol("@subreq")
        macro_node = copy(node)
        node.head = :call
        expanded = macroexpand(POMDPLinter, macro_node)
        if isa(expanded, Expr) && expanded.head == :error
            throw(expanded.args[1])
        end
        node.args = [:push_dep!, reqs_name, esc(macroexpand(POMDPLinter, expanded))]
        return true
    else
        found = falses(length(node.args))
        for (i, arg) in enumerate(node.args)
            found[i] = handle_reqs!(arg, reqs_name)
        end
        if any(found)
            for i in 1:length(node.args)
                if !found[i] # && !(isa(node.args[i], Expr) && node.args[i].head == :line) # this would not escape lines (I don't know what implications that has)
                    node.args[i] = esc(node.args[i])
                end
            end
        end
        return any(found)
    end
end

function handle_reqs!(node::Any, reqs_name::Symbol)
    # for anything that's not an Expr
    return false
end

"""
    @impl_dep reward(::P,::S,::A,::S) where {P<:POMDP,S,A} reward(::P,::S,::A)

Declare an implementation dependency and automatically implement `implemented`.

In the example above, `@implemented reward(::P,::S,::A,::S)` will return true if the user has implemented `reward(::P,::S,::A,::S)` OR `reward(::P,::S,::A)`

THIS IS ONLY INTENDED FOR USE INSIDE POMDPLINTER AND MAY NOT FUNCTION CORRECTLY ELSEWHERE
"""
macro impl_dep(signature, dependency, module_mask=[])
    if signature.head == :where
        sig_req = signature.args[1]
        wheres = signature.args[2:end]
    else
        sig_req = signature
        wheres = ()
    end
    tplex = convert_req(sig_req)
    deptplex = convert_req(dependency)
    impled = quote
        function implemented(f::typeof(first($tplex)), TT::Type{last($tplex)}) where {$(wheres...)}
            m = which(f,TT)
            if m.module in $(module_mask) && !implemented($deptplex...)
                return false
            else # a more specific implementation exists
                return true
            end
            return false
        end
    end
    return esc(impled)
end


text:
mutable struct PReLU{T, W, O} <: AbstractBlock
    weight::W
    output::O
    PReLU{T, W, O}(weight::W) where {T, W, O} = new{T, W, O}(weight)
end

PReLU(::Type{T}, ::Type{O}, weight::W) where {T, W, O} = PReLU{T, W, O}(weight)

PReLU(::Type{T}, weight::W; nbatch=1) where {T, W} =
    PReLU(T, nbatch > 1 ? Matrix{T} : Vector{T}, weight)
PReLU(weight::T;nbatch=1) where T = PReLU(T, weight; nbatch=nbatch)
PReLU(::Type{T};nbatch=1, init=0.25, nparam=1) where T =
    PReLU(nparam > 1 ? fill(T(init), nparam) : T(init);nbatch=nbatch)
PReLU(;nbatch=1, init=0.25, nparam=1) =
    PReLU(Float64; nbatch=nbatch, init=init, nparam=nparam)

_prelu(x::T, w::T) where {T <: Real} = x > 0 ? x : w * x
_prelu_grad(x::T, dx::T, w::T) where {T <: Real} = x > 0 ? dx : w * dx

# single output plane
function forward(op::PReLU{T, T, O}, input::O) where{T <: Real, O}
    return _prelu.(input, op.weight)
end

function backward(op::PReLU{T, T, O}, grad::O) where {T <: Real, O}
    return _prelu_grad.(op.input, grad, op.weight)
end

# TODO: multiple output plane, ref: THNN/generic/PReLU.c
# TODO: complex PReLU


text:
function _write_fields(list::Vector{Tuple{AnticNumberField, fmpz}}, filename::String)
  f=open(filename, "a")
  for L in list
    x=([coeff(L[1].pol, i) for i=0:degree(L[1].pol)], L[2])
    Base.write(f, "$x\n")
  end
  close(f)
end

function __to_univariate(Qx, f)
  z = zero(Qx)
  x = gen(Qx)
  for (c, m) in zip(coefficients(f), monomials(f))
    z = z + FlintQQ(c) * x^total_degree(m)
  end
  return z
end

function _write_fields(list::Vector{Tuple{AnticNumberField, NfRelNS{nf_elem}, fmpz}}, filename::String)
  f=open(filename, "a")
  for L in list
    Qx = parent(L[1].pol)
    unis = [ __to_univariate(Qx, L[2].pol[i]) for i in 1:length(L[2].pol)]
    y = [[coeff(g, i) for i=0:degree(g)] for g in unis]
    x=([coeff(L[1].pol, i) for i=0:degree(L[1].pol)], y, L[3])
    Base.write(f, "$x\n")
  end
  close(f)
end

function _read_fields(filename::String)
  f=open(filename, "r")
  Qx,x=PolynomialRing(FlintQQ,"x")
  pols = []
  for s in eachline(f)
    a=Main.eval(Meta.parse(s))
    if length(a) == 3
      push!(pols,(Qx(a[1]), [ Qx(b) for b in a[2]], a[3]))
    else
      @assert length(a) == 2
      push!(pols,(Qx(a[1]), a[end]))
    end
  end
  close(f)
  return identity.(pols)
end

global small_solvable_groups =
[ (1, 1), (2, 1), (3, 1), (4, 1), (4, 2), (5, 1), (6, 1), (6, 2), (7, 1), (8,
1), (8, 2), (8, 3), (8, 4), (8, 5), (9, 1), (9, 2), (10, 1), (10, 2), (11, 1),
(12, 1), (12, 2), (12, 3), (12, 4), (12, 5), (13, 1), (14, 1), (14, 2), (15, 1),
(16, 1), (16, 2), (16, 3), (16, 4), (16, 5), (16, 6), (16, 7), (16, 8), (16, 9),
(16, 10), (16, 11), (16, 12), (16, 13), (16, 14), (17, 1), (18, 1), (18, 2),
(18, 3), (18, 4), (18, 5), (19, 1), (20, 1), (20, 2), (20, 3), (20, 4), (20, 5),
(21, 1), (21, 2), (22, 1), (22, 2), (23, 1), (24, 1), (24, 2), (24, 3), (24, 4),
(24, 5), (24, 6), (24, 7), (24, 8), (24, 9), (24, 10), (24, 11), (24, 12), (24,
13), (24, 14), (24, 15), (25, 1), (25, 2), (26, 1), (26, 2), (27, 1), (27, 2),
(27, 3), (27, 4), (27, 5), (28, 1), (28, 2), (28, 3), (28, 4), (29, 1), (30, 1),
(30, 2), (30, 3), (30, 4), (31, 1), (32, 1), (32, 2), (32, 3), (32, 4), (32, 5),
(32, 6), (32, 7), (32, 8), (32, 9), (32, 10), (32, 11), (32, 12), (32, 13), (32,
14), (32, 15), (32, 16), (32, 17), (32, 18), (32, 19), (32, 20), (32, 21), (32,
22), (32, 23), (32, 24), (32, 25), (32, 26), (32, 27), (32, 28), (32, 29), (32,
30), (32, 31), (32, 32), (32, 33), (32, 34), (32, 35), (32, 36), (32, 37), (32,
38), (32, 39), (32, 40), (32, 41), (32, 42), (32, 43), (32, 44), (32, 45), (32,
46), (32, 47), (32, 48), (32, 49), (32, 50), (32, 51), (33, 1), (34, 1), (34,
2), (35, 1), (36, 1), (36, 2), (36, 3), (36, 4), (36, 5), (36, 6), (36, 7), (36,
8), (36, 9), (36, 10), (36, 11), (36, 12), (36, 13), (36, 14), (37, 1), (38, 1),
(38, 2), (39, 1), (39, 2), (40, 1), (40, 2), (40, 3), (40, 4), (40, 5), (40, 6),
(40, 7), (40, 8), (40, 9), (40, 10), (40, 11), (40, 12), (40, 13), (40, 14),
(41, 1), (42, 1), (42, 2), (42, 3), (42, 4), (42, 5), (42, 6), (43, 1), (44, 1),
(44, 2), (44, 3), (44, 4), (45, 1), (45, 2), (46, 1), (46, 2), (47, 1), (48, 1),
(48, 2), (48, 3), (48, 4), (48, 5), (48, 6), (48, 7), (48, 8), (48, 9), (48,
10), (48, 11), (48, 12), (48, 13), (48, 14), (48, 15), (48, 16), (48, 17), (48,
18), (48, 19), (48, 20), (48, 21), (48, 22), (48, 23), (48, 24), (48, 25), (48,
26), (48, 27), (48, 28), (48, 29), (48, 30), (48, 31), (48, 32), (48, 33), (48,
34), (48, 35), (48, 36), (48, 37), (48, 38), (48, 39), (48, 40), (48, 41), (48,
42), (48, 43), (48, 44), (48, 45), (48, 46), (48, 47), (48, 48), (48, 49), (48,
50), (48, 51), (48, 52), (49, 1), (49, 2), (50, 1), (50, 2), (50, 3), (50, 4),
(50, 5), (51, 1), (52, 1), (52, 2), (52, 3), (52, 4), (52, 5), (53, 1), (54, 1),
(54, 2), (54, 3), (54, 4), (54, 5), (54, 6), (54, 7), (54, 8), (54, 9), (54,
10), (54, 11), (54, 12), (54, 13), (54, 14), (54, 15), (55, 1), (55, 2), (56,
1), (56, 2), (56, 3), (56, 4), (56, 5), (56, 6), (56, 7), (56, 8), (56, 9), (56,
10), (56, 11), (56, 12), (56, 13), (57, 1), (57, 2), (58, 1), (58, 2), (59, 1),
(60, 1), (60, 2), (60, 3), (60, 4), (60, 6), (60, 7), (60, 8), (60, 9), (60,
10), (60, 11), (60, 12), (60, 13), (61, 1), (62, 1), (62, 2), (63, 1), (63, 2),
(63, 3), (63, 4), (64, 1), (64, 2), (64, 3), (64, 4), (64, 5), (64, 6), (64, 7),
(64, 8), (64, 9), (64, 10), (64, 11), (64, 12), (64, 13), (64, 14), (64, 15),
(64, 16), (64, 17), (64, 18), (64, 19), (64, 20), (64, 21), (64, 22), (64, 23),
(64, 24), (64, 25), (64, 26), (64, 27), (64, 28), (64, 29), (64, 30), (64, 31),
(64, 32), (64, 33), (64, 34), (64, 35), (64, 36), (64, 37), (64, 38), (64, 39),
(64, 40), (64, 41), (64, 42), (64, 43), (64, 44), (64, 45), (64, 46), (64, 47),
(64, 48), (64, 49), (64, 50), (64, 51), (64, 52), (64, 53), (64, 54), (64, 55),
(64, 56), (64, 57), (64, 58), (64, 59), (64, 60), (64, 61), (64, 62), (64, 63),
(64, 64), (64, 65), (64, 66), (64, 67), (64, 68), (64, 69), (64, 70), (64, 71),
(64, 72), (64, 73), (64, 74), (64, 75), (64, 76), (64, 77), (64, 78), (64, 79),
(64, 80), (64, 81), (64, 82), (64, 83), (64, 84), (64, 85), (64, 86), (64, 87),
(64, 88), (64, 89), (64, 90), (64, 91), (64, 92), (64, 93), (64, 94), (64, 95),
(64, 96), (64, 97), (64, 98), (64, 99), (64, 100), (64, 101), (64, 102), (64,
103), (64, 104), (64, 105), (64, 106), (64, 107), (64, 108), (64, 109), (64,
110), (64, 111), (64, 112), (64, 113), (64, 114), (64, 115), (64, 116), (64,
117), (64, 118), (64, 119), (64, 120), (64, 121), (64, 122), (64, 123), (64,
124), (64, 125), (64, 126), (64, 127), (64, 128), (64, 129), (64, 130), (64,
131), (64, 132), (64, 133), (64, 134), (64, 135), (64, 136), (64, 137), (64,
138), (64, 139), (64, 140), (64, 141), (64, 142), (64, 143), (64, 144), (64,
145), (64, 146), (64, 147), (64, 148), (64, 149), (64, 150), (64, 151), (64,
152), (64, 153), (64, 154), (64, 155), (64, 156), (64, 157), (64, 158), (64,
159), (64, 160), (64, 161), (64, 162), (64, 163), (64, 164), (64, 165), (64,
166), (64, 167), (64, 168), (64, 169), (64, 170), (64, 171), (64, 172), (64,
173), (64, 174), (64, 175), (64, 176), (64, 177), (64, 178), (64, 179), (64,
180), (64, 181), (64, 182), (64, 183), (64, 184), (64, 185), (64, 186), (64,
187), (64, 188), (64, 189), (64, 190), (64, 191), (64, 192), (64, 193), (64,
194), (64, 195), (64, 196), (64, 197), (64, 198), (64, 199), (64, 200), (64,
201), (64, 202), (64, 203), (64, 204), (64, 205), (64, 206), (64, 207), (64,
208), (64, 209), (64, 210), (64, 211), (64, 212), (64, 213), (64, 214), (64,
215), (64, 216), (64, 217), (64, 218), (64, 219), (64, 220), (64, 221), (64,
222), (64, 223), (64, 224), (64, 225), (64, 226), (64, 227), (64, 228), (64,
229), (64, 230), (64, 231), (64, 232), (64, 233), (64, 234), (64, 235), (64,
236), (64, 237), (64, 238), (64, 239), (64, 240), (64, 241), (64, 242), (64,
243), (64, 244), (64, 245), (64, 246), (64, 247), (64, 248), (64, 249), (64,
250), (64, 251), (64, 252), (64, 253), (64, 254), (64, 255), (64, 256), (64,
257), (64, 258), (64, 259), (64, 260), (64, 261), (64, 262), (64, 263), (64,
264), (64, 265), (64, 266), (64, 267), (65, 1), (66, 1), (66, 2), (66, 3), (66,
4), (67, 1), (68, 1), (68, 2), (68, 3), (68, 4), (68, 5), (69, 1), (70, 1), (70,
2), (70, 3), (70, 4), (71, 1), (72, 1), (72, 2), (72, 3), (72, 4), (72, 5), (72,
6), (72, 7), (72, 8), (72, 9), (72, 10), (72, 11), (72, 12), (72, 13), (72, 14),
(72, 15), (72, 16), (72, 17), (72, 18), (72, 19), (72, 20), (72, 21), (72, 22),
(72, 23), (72, 24), (72, 25), (72, 26), (72, 27), (72, 28), (72, 29), (72, 30),
(72, 31), (72, 32), (72, 33), (72, 34), (72, 35), (72, 36), (72, 37), (72, 38),
(72, 39), (72, 40), (72, 41), (72, 42), (72, 43), (72, 44), (72, 45), (72, 46),
(72, 47), (72, 48), (72, 49), (72, 50), (73, 1), (74, 1), (74, 2), (75, 1), (75,
2), (75, 3), (76, 1), (76, 2), (76, 3), (76, 4), (77, 1), (78, 1), (78, 2), (78,
3), (78, 4), (78, 5), (78, 6), (79, 1), (80, 1), (80, 2), (80, 3), (80, 4), (80,
5), (80, 6), (80, 7), (80, 8), (80, 9), (80, 10), (80, 11), (80, 12), (80, 13),
(80, 14), (80, 15), (80, 16), (80, 17), (80, 18), (80, 19), (80, 20), (80, 21),
(80, 22), (80, 23), (80, 24), (80, 25), (80, 26), (80, 27), (80, 28), (80, 29),
(80, 30), (80, 31), (80, 32), (80, 33), (80, 34), (80, 35), (80, 36), (80, 37),
(80, 38), (80, 39), (80, 40), (80, 41), (80, 42), (80, 43), (80, 44), (80, 45),
(80, 46), (80, 47), (80, 48), (80, 49), (80, 50), (80, 51), (80, 52), (81, 1),
(81, 2), (81, 3), (81, 4), (81, 5), (81, 6), (81, 7), (81, 8), (81, 9), (81,
10), (81, 11), (81, 12), (81, 13), (81, 14), (81, 15), (82, 1), (82, 2), (83,
1), (84, 1), (84, 2), (84, 3), (84, 4), (84, 5), (84, 6), (84, 7), (84, 8), (84,
9), (84, 10), (84, 11), (84, 12), (84, 13), (84, 14), (84, 15), (85, 1), (86,
1), (86, 2), (87, 1), (88, 1), (88, 2), (88, 3), (88, 4), (88, 5), (88, 6), (88,
7), (88, 8), (88, 9), (88, 10), (88, 11), (88, 12), (89, 1), (90, 1), (90, 2),
(90, 3), (90, 4), (90, 5), (90, 6), (90, 7), (90, 8), (90, 9), (90, 10), (91,
1), (92, 1), (92, 2), (92, 3), (92, 4), (93, 1), (93, 2), (94, 1), (94, 2), (95,
1), (96, 1), (96, 2), (96, 3), (96, 4), (96, 5), (96, 6), (96, 7), (96, 8), (96,
9), (96, 10), (96, 11), (96, 12), (96, 13), (96, 14), (96, 15), (96, 16), (96,
17), (96, 18), (96, 19), (96, 20), (96, 21), (96, 22), (96, 23), (96, 24), (96,
25), (96, 26), (96, 27), (96, 28), (96, 29), (96, 30), (96, 31), (96, 32), (96,
33), (96, 34), (96, 35), (96, 36), (96, 37), (96, 38), (96, 39), (96, 40), (96,
41), (96, 42), (96, 43), (96, 44), (96, 45), (96, 46), (96, 47), (96, 48), (96,
49), (96, 50), (96, 51), (96, 52), (96, 53), (96, 54), (96, 55), (96, 56), (96,
57), (96, 58), (96, 59), (96, 60), (96, 61), (96, 62), (96, 63), (96, 64), (96,
65), (96, 66), (96, 67), (96, 68), (96, 69), (96, 70), (96, 71), (96, 72), (96,
73), (96, 74), (96, 75), (96, 76), (96, 77), (96, 78), (96, 79), (96, 80), (96,
81), (96, 82), (96, 83), (96, 84), (96, 85), (96, 86), (96, 87), (96, 88), (96,
89), (96, 90), (96, 91), (96, 92), (96, 93), (96, 94), (96, 95), (96, 96), (96,
97), (96, 98), (96, 99), (96, 100), (96, 101), (96, 102), (96, 103), (96, 104),
(96, 105), (96, 106), (96, 107), (96, 108), (96, 109), (96, 110), (96, 111),
(96, 112), (96, 113), (96, 114), (96, 115), (96, 116), (96, 117), (96, 118),
(96, 119), (96, 120), (96, 121), (96, 122), (96, 123), (96, 124), (96, 125),
(96, 126), (96, 127), (96, 128), (96, 129), (96, 130), (96, 131), (96, 132),
(96, 133), (96, 134), (96, 135), (96, 136), (96, 137), (96, 138), (96, 139),
(96, 140), (96, 141), (96, 142), (96, 143), (96, 144), (96, 145), (96, 146),
(96, 147), (96, 148), (96, 149), (96, 150), (96, 151), (96, 152), (96, 153),
(96, 154), (96, 155), (96, 156), (96, 157), (96, 158), (96, 159), (96, 160),
(96, 161), (96, 162), (96, 163), (96, 164), (96, 165), (96, 166), (96, 167),
(96, 168), (96, 169), (96, 170), (96, 171), (96, 172), (96, 173), (96, 174),
(96, 175), (96, 176), (96, 177), (96, 178), (96, 179), (96, 180), (96, 181),
(96, 182), (96, 183), (96, 184), (96, 185), (96, 186), (96, 187), (96, 188),
(96, 189), (96, 190), (96, 191), (96, 192), (96, 193), (96, 194), (96, 195),
(96, 196), (96, 197), (96, 198), (96, 199), (96, 200), (96, 201), (96, 202),
(96, 203), (96, 204), (96, 205), (96, 206), (96, 207), (96, 208), (96, 209),
(96, 210), (96, 211), (96, 212), (96, 213), (96, 214), (96, 215), (96, 216),
(96, 217), (96, 218), (96, 219), (96, 220), (96, 221), (96, 222), (96, 223),
(96, 224), (96, 225), (96, 226), (96, 227), (96, 228), (96, 229), (96, 230),
(96, 231), (97, 1), (98, 1), (98, 2), (98, 3), (98, 4), (98, 5), (99, 1), (99,
2), (100, 1), (100, 2), (100, 3), (100, 4), (100, 5), (100, 6), (100, 7), (100,
8), (100, 9), (100, 10), (100, 11), (100, 12), (100, 13), (100, 14), (100, 15),
(100, 16) ]


text:
# TODO: Optimize cache size
mutable struct AN5ConstantCache{zType,lType,dtsType,dType,tsit5Type} <: OrdinaryDiffEqConstantCache
  # `z` is the Nordsieck vector
  z::zType
  # `l` is used for the corrector iteration
  l::Vector{lType}
  # `m` is a tmp vector that is used for calculating `l`
  m::Vector{lType}
  # `c_LTE` is used for the error estimation for the current order
  c_LTE::lType
  c_conv::lType
  # `dts` stores `dt`s
  dts::dtsType
  # `Δ` is the difference between the predictor `uₙ₀` and `uₙ`
  Δ::dType
  # `Tsit5` for the first step
  tsit5tab::tsit5Type
  order::Int
end

function alg_cache(alg::AN5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  N = 5
  z = [zero(rate_prototype) for i in 1:N+1]
  Δ = u
  l = fill(zero(tTypeNoUnits),N+1); m = zero(l)
  c_LTE = c_conv = zero(tTypeNoUnits)
  dts = fill(zero(dt), 6)
  tsit5tab = Tsit5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  AN5ConstantCache(z,l,m,c_LTE,c_conv,dts,Δ,tsit5tab,1)
end

mutable struct AN5Cache{uType,dType,rateType,zType,lType,dtsType,tsit5Type} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  Δ::dType
  # Error estimation
  atmp::dType
  fsalfirst::rateType
  ratetmp::rateType
  # `z` is the Nordsieck vector
  z::zType
  # `l` is used for the corrector iteration
  l::Vector{lType}
  # `m` is a tmp vector that is used for calculating `l`
  m::Vector{lType}
  # `c_LTE` is used for the error estimation for the current order
  c_LTE::lType
  c_conv::lType
  # `dts` stores `dt`s
  dts::dtsType
  # `Tsit5` for the first step
  tsit5cache::tsit5Type
  order::Int
end

function alg_cache(alg::AN5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  #################################################
  # Tsit5
  tab = Tsit5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  # Cannot alias pointers, since we have to use `k`s to start the Nordsieck vector
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)
  k3 = zero(rate_prototype); k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype); k7 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits); tmp = similar(u)
  tsit5cache = Tsit5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
  #################################################
  N = 5
  Δ = similar(atmp)
  l = fill(zero(tTypeNoUnits),N+1); m = zero(l)
  c_LTE = c_conv = zero(tTypeNoUnits)
  dts = fill(zero(dt), 6)
  fsalfirst = zero(rate_prototype)
  z = [zero(rate_prototype) for i in 1:N+1]
  for i in 1:N+1
    z[i] = zero(rate_prototype)
  end
  ratetmp = zero(rate_prototype)

  AN5Cache(u,uprev,tmp,Δ,atmp,fsalfirst,ratetmp,
           z,l,m,c_LTE,c_conv,dts,
           tsit5cache, 1)
end

mutable struct JVODEConstantCache{zType,lType,dtsType,dType,tsit5Type,etaType} <: OrdinaryDiffEqConstantCache
  # `z` is the Nordsieck vector
  z::zType
  # `l` is used for the corrector iteration
  l::Vector{lType}
  # `m` is a tmp vector that is used for calculating `l`
  m::Vector{lType}
  # `c_LTE₊₁` is used for the error estimation for the current order + 1
  c_LTE₊₁::lType
  # `c_LTE` is used for the error estimation for the current order
  c_LTE::lType
  # `c_LTE₋₁` is used for the error estimation for the current order - 1
  c_LTE₋₁::lType
  # `c_conv` is used in convergence test
  c_conv::lType
  # `c_𝒟` is used to get the order q+2 derivative vector
  c_𝒟::lType
  prev_𝒟::lType
  # `dts` stores `dt`s
  dts::dtsType
  # `Δ` is the difference between the predictor `uₙ₀` and `uₙ`
  Δ::dType
  # `Tsit5` for the first step
  tsit5tab::tsit5Type
  L::Int
  # same with `order` or `q`
  order::Int
  nextorder::Int
  # number of steps to take before considering to change order
  n_wait::Int
  # `η` is `dtₙ₊₁/dtₙ`
  η  ::etaType
  ηq ::etaType
  η₊₁::etaType
  η₋₁::etaType
  maxη::etaType
end

function alg_cache(alg::JVODE,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  N = 12
  z = [rate_prototype for i in 1:N+1]
  Δ = u
  l = fill(zero(tTypeNoUnits), N+1); m = zero(l)
  c_LTE₊₁ = c_LTE = c_LTE₋₁ = c_conv = c_𝒟 = prev_𝒟 = zero(tTypeNoUnits)
  dts = fill(zero(dt),N+1)
  tsit5tab = Tsit5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  η = zero(dt/dt)
  JVODEConstantCache(z,l,m,
                     c_LTE₊₁,c_LTE,c_LTE₋₁,c_conv,c_𝒟 ,prev_𝒟,
                     dts,Δ,tsit5tab,2,1,1,2,η,η,η,η,η)
end

mutable struct JVODECache{uType,rateType,zType,lType,dtsType,dType,etaType,tsit5Type} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  fsalfirst::rateType
  ratetmp::rateType
  # `z` is the Nordsieck vector
  z::zType
  # `l` is used for the corrector iteration
  l::Vector{lType}
  # `m` is a tmp vector that is used for calculating `l`
  m::Vector{lType}
  # `c_LTE₊₁` is used for the error estimation for the current order + 1
  c_LTE₊₁::lType
  # `c_LTE` is used for the error estimation for the current order
  c_LTE::lType
  # `c_LTE₋₁` is used for the error estimation for the current order - 1
  c_LTE₋₁::lType
  # `c_conv` is used in convergence test
  c_conv::lType
  # `c_𝒟` is used to get the order q+2 derivative vector
  c_𝒟::lType
  prev_𝒟::lType
  # `dts` stores `dt`s
  dts::dtsType
  # `Δ` is the difference between the predictor `uₙ₀` and `uₙ`
  Δ::dType
  # Error estimation
  atmp::dType
  # `Tsit5` for the first step
  tsit5cache::tsit5Type
  L::Int
  # same with `order` or `q`
  order::Int
  nextorder::Int
  # number of steps to take before considering to change order
  n_wait::Int
  # `η` is `dtₙ₊₁/dtₙ`
  η  ::etaType
  ηq ::etaType
  η₊₁::etaType
  η₋₁::etaType
  maxη::etaType
end

function alg_cache(alg::JVODE,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  #################################################
  # Tsit5
  # Cannot alias pointers, since we have to use `k`s to start the Nordsieck vector
  tab = Tsit5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype); k4 = zero(rate_prototype)
  k5 = zero(rate_prototype); k6 = zero(rate_prototype); k7 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits); tmp = similar(u)
  tsit5cache = Tsit5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
  #################################################
  fsalfirst = zero(rate_prototype)
  N = 12
  z = [zero(rate_prototype) for i in 1:N+1]
  Δ = similar(u,uEltypeNoUnits)
  l = fill(zero(tTypeNoUnits), N+1); m = zero(l)
  c_LTE₊₁ = c_LTE = c_LTE₋₁ = c_conv = c_𝒟 = prev_𝒟 = zero(tTypeNoUnits)
  dts = fill(zero(dt),N+1)
  η = zero(dt/dt)
  #################################################
  # Nordsieck Vector
  z[1] = zero(rate_prototype); z[2] = zero(rate_prototype); z[3] = zero(rate_prototype);
  z[4] = zero(rate_prototype); z[5] = zero(rate_prototype); z[6] = zero(rate_prototype);
  z[7] = zero(rate_prototype); z[8] = zero(rate_prototype); z[9] = zero(rate_prototype);
  z[10] = zero(rate_prototype); z[11] = zero(rate_prototype); z[12] = zero(rate_prototype);
  z[13] = zero(rate_prototype)
  ratetmp = zero(rate_prototype)
  #################################################
  JVODECache(u,uprev,tmp,fsalfirst,ratetmp,
             z, l, m,
             c_LTE₊₁, c_LTE, c_LTE₋₁, c_conv, c_𝒟 , prev_𝒟,
             dts, Δ, atmp, tsit5cache, 2, 1, 1, 2, η, η, η, η, η)
end