text:
"""
format(str::AbstractString, args...) -> String

format `str` with `args`
"""
function format(str::AbstractString, args...; kwargs...)
    inbetweens = String[]
    fields = String[]
    field = ""
    nextidx = 1
    while field !== nothing
        nextidx, inbetween, field = iterate_braces_pair(str, nextidx)

        push!(inbetweens, inbetween)
        if field !== nothing
            push!(fields, field)
        end
    end

    #TODO: do replacements in fmtstr
    #getattr and getindex

    finalstr = ""

    for (i, field) in enumerate(fields)
        idx = findfirst(':', field)
        if idx === nothing
            field_name = field
            specifier = ""
        else
            field_name = field[1:(idx-1)]
            specifier = field[(idx+1):end]
        end
        formated =
            if isempty(field_name)
                fmt(specifier, args[i])
            elseif isdigit(field_name[1])
                fmt(specifier, args[parse(Int,field_name)])
            else
                fmt(specifier, kwargs[Symbol(field_name)])
            end
        finalstr *= inbetweens[i] * formated
    end
    finalstr * inbetweens[end]
end

"""
iterate_braces_pair(str, start)

return end of analysed string, in-between string and the format string
"""
function iterate_braces_pair(str, start)
    i = start
    while i <= lastindex(str)
        if str[i] == '}'
            if i == lastindex(str)
                error("single '}' is not allowed at idx $i in string \"$str\"")
            elseif str[i+1] == '}'
                i = nextind(str, i+1)
                continue
            else
                error("single '}' is not allowed at idx $i in string \"$str\"")
            end
        end

        if str[i] == '{'
            if i == lastindex(str)
                error("single '{' is not allowed at idx $i in string \"$str\"")
            elseif str[i+1] == '{'
                i = nextind(str, i+1)
                continue
            else
                break
            end
        end

        i = nextind(str, i)
    end

    count = 1
    j = i+1
    while j <= lastindex(str)
        if str[j] == '{'
            count += 1
        elseif str[j] == '}'
            count -= 1
        end
        if count == 0
            break
        end
        j = nextind(str, j)
    end

    found = i <= lastindex(str) && str[i] == '{'

    if count != 0 && found
        error("braces missmatch: \"$(str[i:end])\" at idx $i")
    end
    (
        j+1,
        replace(replace(str[start:(i-1)], "{{"=>"{"), "}}" => "}"),
        found ? str[i+1:j-1] : nothing
    )
end


text:
# dependencies
using LFAToolkit
using LinearAlgebra
using CSV
using DataFrames

# setup
dimension = 1
numbercomponents = 1
mesh = []
if dimension == 1
    mesh = Mesh1D(1.0)
elseif dimension == 2
    mesh = Mesh2D(1.0, 1.0)
elseif dimension == 3
    mesh = Mesh3D(1.0, 1.0, 1.0)
end
smoothingparameters = DataFrame()

# test range
for P = 1:4
    println("p = ", 2^P)
    # setup
    # -- bases
    p = 2^P

    # -- diffusion operator
    diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

    # -- Jacobi smoother
    jacobi = Jacobi(diffusion)

    # compute smoothing factor
    # -- setup
    numbersteps = 100
    maxeigenvalue = 0
    θ_min = -π / 2
    θ_min_high = π / 2
    θ_max = 3π / 2
    θ_step = 2π / numbersteps
    θ_range = θ_min:θ_step:(θ_max-θ_step)
    ω = 1.0

    # -- compute
    λ_minhigh = 100
    λ_maxhigh = -100
    λ_max = -100
    # -- 1D --
    if dimension == 1
        for i = 1:numbersteps
            θ = [θ_range[i]]
            if abs(θ[1]) > π / 128
                A = computesymbols(jacobi, [ω], θ)
                eigenvalues = [abs(val) for val in eigvals(I - A)]
                if (θ[1] > θ_min_high)
                    λ_minhigh = min(λ_minhigh, eigenvalues...)
                    λ_maxhigh = max(λ_maxhigh, eigenvalues...)
                end
                λ_max = max(λ_max, eigenvalues...)
            end
        end
        # -- 2D --
    elseif dimension == 2
        for i = 1:numbersteps, j = 1:numbersteps
            θ = [θ_range[i], θ_range[j]]
            if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π / 128
                A = computesymbols(jacobi, [ω], θ)
                eigenvalues = [abs(val) for val in eigvals(I - A)]
                if (θ[1] > θ_min_high || θ[2] > θ_min_high)
                    λ_minhigh = min(λ_minhigh, eigenvalues...)
                    λ_maxhigh = max(λ_maxhigh, eigenvalues...)
                end
                λ_max = max(λ_max, eigenvalues...)
            end
        end
        # -- 3D --
    elseif dimension == 3
        for i = 1:numbersteps, j = 1:numbersteps, k = 1:numbersteps
            θ = [θ_range[i], θ_range[j], θ_range[k]]
            if sqrt(abs(θ[1])^2 + abs(θ[2])^2 + abs(θ[3])^2) > π / 128
                A = computesymbols(jacobi, [ω], θ)
                eigenvalues = [abs(val) for val in eigvals(I - A)]
                if (θ[1] > θ_min_high || θ[2] > θ_min_high)
                    λ_minhigh = min(λ_minhigh, eigenvalues...)
                    λ_maxhigh = max(λ_maxhigh, eigenvalues...)
                end
                λ_max = max(λ_max, eigenvalues...)
            end
        end
    end
    ω_classical = 2 / (λ_minhigh + λ_maxhigh)
    ω_highorder = 2 / (λ_minhigh + λ_max)
    append!(
        smoothingparameters,
        DataFrame(p = p, ω_classical = ω_classical, ω_highorder = ω_highorder),
    )
end

CSV.write("jacobi-smoothing-factor-estimate.csv", smoothingparameters)

ext: jl


text:
# based on examples like https://gitlab.eeecs.qub.ac.uk/40126401/csc4006-EdgeBenchmarking/blob/f56338fe1c6ae5b87a8ba59ea17ac64387e8213b/Experiments/YOLO/yolo/darknet.py
function detect(net, meta, img; thresh=0.5, hier_thresh=0.5, nms=0)
    num_dets = Ref(Cint(0))
    network_predict_image(net, img)
    dets = get_network_boxes(net, img.w, img.h, thresh, hier_thresh, C_NULL, 0, num_dets, 0)
    if nms > 0
        do_nms_sort(dets, num_dets[], meta.classes, nms)
    end
    res = []
    detections = unsafe_wrap(Array, dets, num_dets[])
    names = unsafe_wrap(Array, meta.names, meta.classes)
    for j in 1:num_dets[]
        prob = unsafe_wrap(Array, detections[j].prob, meta.classes)
        for i in 1:meta.classes
            if prob[i] > 0
                b = detections[j].bbox
                nameTag = names[i]
                push!(res,(nameTag, prob[i], (b.x, b.y, b.w, b.h)))
            end
        end
    end
    free_detections(dets, num_dets[])
    return res
end


text:
#
println("At) Test of parallel computing in Julia.")
using Distributed
@everywhere using LinearAlgebra

println("Number of cores = $(nprocs()),  number of workers = $(nworkers())")

@everywhere function inside(n)
    ni = 0;     for  i in 1:n   x,y = rand(), rand();   if  x^2 + y^2 <= 1      ni = ni + 1     end     end
    return( ni )
end
#
function piParallel(n)
    p = nworkers();     ni = @distributed (+) for i = 1:p   inside(n/p)     end
    return( 4*ni/n )
end
#
@everywhere function myAction(A,B)
    C = A*B;    return( C[1,1] )
end

nx = 5000
# using Profile
# Profile.clear()
# @profile if  false
@time if false
# if false
    # Sequential processing
    for  i in 1:10
        wa = rand(nx,nx)
        wb = LinearAlgebra.eigvals(wa)
        println("wb = $(wb[1:3])")
    end

elseif true
    # Broadcasting of values and functions
    A = rand(nx,nx);    B = rand(nx,nx);    Cf = Future[]
    @everywhere A, B = $A, $B
    for  i in 1:15     ii = rem(i,nworkers())+1;     push!(Cf, @spawnat ii (A.^17.3)[1,1])       end
    ## for  i in 1:15     rf = remotecall(*, WorkerPool(workers()), A,B);   push!(Cf, rf)   end
    println("first loop done")
    for  i in 1:15     C = fetch(Cf[i]);   println("$i  $(C[1,1])")    end

elseif false
    # Producer-consumer model
    @async foo()

elseif true
    # Distributed calculations
    nx = 1000000000
    @time wc = inside(nx);         println("wc = $wc")
    @time wd = piParallel(nx);     println("wd = $wd")

elseif false

end


text:
module seq2seq

using Flux, ArgCheck, StatsBase, YAML

using ..modelutils
import ..modelutils: load!, randn_repar
export load!

const D_IM = 32

"""
    Seq2SeqModel(encoder, x0_post_dist, x0decoder, generator)

A sequence-to-sequence model 𝒴→𝒴 which imbibes `T_enc` steps of the sequence,
and predicts `T_steps` (including the original `T_enc`).

This type has 5 slots:
1. `encoder`
2. `x0_post_dist`
3. `x0decoder`
4. `generator`
5. `d`

The posterior over the initial state of the generator is learned variationally
via the `encoder` RNN. The final state of this encoder is passed to the
`x0_post_dist` which is a `MultiDense` objects which spits out two values,
i.e. the mean and variance of the latent posterior. Typically the posterior is
lower dimensional than the state of the generator network, so the latent
variable is expanded via the `x0decoder` network into the state of the
`generator`. Finally the `generator` is iterated in an open-loop fashion for
`T_steps` with no inputs to provide the prediction. The final slot `d` is just
a record of the dimension of 𝒴, used currently just for `show`.

--------------------
This type can be called via
```julia
    m = Seq2SeqModel(...)
    m(y::Vector)
```
where `y` is a vector of slices of the sequence for the encoder. This means each
`y[i]` is an `AbstractMatrix` of size ``d \\times n_{\\text{batch}}``.

Options include `T_enc`, `T_steps` (already discussed), `reset` (which performs
a state reset before execution; true by default), and `stoch`, whether to
randomly sample the latent state, or use the state mean (default true).
"""
struct Seq2SeqModel{S,U,V,W}
    encoder::S
    x0_post_dist::U
    x0decoder::V
    generator::W
    d::Int
end

Flux.@treelike Seq2SeqModel
unpack(m::Seq2SeqModel) = Flux.children(m)[1:4] # unpack struct; exclude `d`.

modelutils.load!(m::Seq2SeqModel, fname::String) = load!(Flux.params(m), fname)

function Base.show(io::IO, l::Seq2SeqModel)
    has_cnn = length(l.encoder) != 1
    out_type = has_cnn ? "CNN" : l.generator.layers[end][end] isa MultiDense ? "Probabilistic" : "Deterministic"
    print(io, "Seq2SeqModel(dim=", l.d, ", ", out_type, ")")
end

# ≡ __call__
function (m::Seq2SeqModel)(y::AbstractVector; T_steps=70, T_enc=10, reset=true, stoch=true)
    encoder, x0_post_dist, x0decoder, generator = unpack(m)
    reset && Flux.reset!(encoder)

    [encoder(y[tt]) for tt in 1:T_enc];
    enc_state = encoder.layers[end].state # d × d × nbatch → d_enc_state
    enc_state = enc_state isa Tuple ? enc_state[1] : enc_state  # handle LSTM
    μ_x0, σ_x0 = x0_post_dist(enc_state)  # 2×Dense: d_enc_state × nbatch → d_x0 × nbatch
    d_x0, nbatch = size(μ_x0)
    x0_smp = μ_x0 + randn_repar(σ_x0, d_x0, nbatch, stoch)

    # Set state (no need to `reset'!)
    generator.layers[1].state = x0decoder(x0_smp)  # RHS: d_x0 × nbatch → d_x × nbatch
    gmove = Flux.has_cuarrays() ? gpu : identity   # isn't super important, can probably make `gpu`
    yhats = [generator(gmove(zeros(1, nbatch))) for t in 1:T_steps]
end

function nllh(m::Seq2SeqModel, y::AbstractVector; T_steps=70, T_enc=10, reset=true, stoch=true)
    yhats = m(y; T_steps=T_steps, T_enc=T_enc, reset=reset, stoch=stoch)
    sum([sum(Flux.logitbinarycrossentropy.(ŷ, y)) for (ŷ, y) in zip(yhats, y)])
end

StatsBase.loglikelihood(m::Seq2SeqModel, y::AbstractVector; T_steps=70, T_enc=10, reset=true,
    stoch=true) = -nllh(m, y; T_steps=T_steps, T_enc=T_enc, reset=reset, stoch=stoch)

"""
    create_model(d_x, d_x0, d_y, d_enc_state; encoder=:LSTM,
    cnn=false, out_heads=1, d_hidden=d_x)

Create `Seq2SeqModel` with the number of hidden units in the recurrent
generative model as `d_x`, the size of the latent variable for the initial state
as `d_x0` (typically ``d_x0 \\ll d_x``), `d_y` the dimension of the observations
and `d_enc_state` the number of hidden units in the recurrent encoder. This
encoder can be specified as `:LSTM`, `:GRU` or `:Bidirectional`. If the
observations are video data and one wishes to use a CNN to encode and decode,
specify `cnn=true`. By default the `Seq2SeqModel` looks to make a point estimate
of the future series, but specifying `out_heads=2` will result in uncertainty
estimation too. Finally `d_hidden` specifies the number of hidden units in non
CNN decoders.
"""
function create_model(d_x, d_x0, d_y, d_enc_state; encoder=:LSTM,
    cnn=false, out_heads=1, d_hidden=d_x)

    @assert !(out_heads > 1 && cnn) "cannot have multiple output heads and CNN."
    @argcheck out_heads in 1:2

    # for trained models cnn_enc_dec=(32,32,32)
    if cnn
        N_FILT=32
        d_conv_result = 4*4*N_FILT
        tform_enc = Chain(
            Conv((3, 3), 2=>N_FILT, relu, stride=2, pad=1),        # out: 16 × 16 × 32 × nbatch
            Conv((3, 3), N_FILT=>N_FILT, relu, stride=2, pad=1),   # out: 8 × 8 × 32 × nbatch
            Conv((3, 3), N_FILT=>N_FILT, relu, stride=2, pad=1),   # out: 4 × 4 × 32 × nbatch
            x->reshape(x, d_conv_result, :)
        )
        d_rnn_in = d_conv_result
    else
        tform_enc = Chain()
        d_rnn_in = d_y
    end

    if encoder == :LSTM
        init_enc = LSTM(d_rnn_in, d_enc_state)
    elseif encoder == :GRU
        init_enc = GRU(d_rnn_in, d_enc_state)
    elseif encoder == :Bidirectional
        init_enc = BRNNenc(d_rnn_in, d_enc_state)
        d_enc_state *= 2
    else
        error("encoder must be specified as :LSTM, :GRU, :Bidirectional")
    end

    enc_rnn = Chain(tform_enc.layers..., init_enc)

    post_x0 = MultiDense(Dense(d_enc_state, d_x0, identity), Dense(d_enc_state, d_x0, σ))
    post_x0.Dense2.b.data .= -2
    x0decoder = Dense(d_x0, d_x, identity)

    # (tform_enc, init_enc, post_x0, x0decoder)

    gen_rnn = GRU(1, d_x)

    if out_heads == 1 && !cnn
        decoder = Chain(Dense(d_x, d_hidden, relu), Dense(d_hidden, d_y, identity))
    elseif out_heads == 2 && !cnn
        decoder = Chain(Dense(d_x, d_hidden, relu), MultiDense(Dense(d_hidden, d_y, identity), Dense(d_hidden, d_y, identity)))
    else    # cnn
        decoder = Chain(
            Dense(d_x, d_conv_result, identity),
            x->reshape(x, 4, 4, N_FILT, :),                               # out: 4 × 4 × 32 × nbatch
            ConvTranspose((3,3), N_FILT=>N_FILT, relu, stride=2),         # out: 9 × 9 × 16 × nbatch
            ConvTranspose((3,3), N_FILT=>N_FILT, relu, stride=2, pad=1),  # out: 17 × 17 × 8 × nbatch
            ConvTranspose((3,3), N_FILT=>2, identity, stride=2, pad=1),   # out: 33 × 33 × 1 × nbatch
            x->x[1:D_IM, 1:D_IM, :, :]   # out: 32 × 32 × 2 × nbatch
        )
    end
    gen_model = Chain(gen_rnn, decoder)

    return Seq2SeqModel(enc_rnn, post_x0, x0decoder, gen_model, d_y)
end

"""
    create_model_opt_dict(...)
Takes same arguments as `create_model` but just wraps up arguments in a Dict for
saving out as JSON/YAML/BSON etc.
"""
function create_model_opt_dict(d_x, d_x0, d_y, d_enc_state; encoder=:LSTM,
    cnn=false, out_heads=1, d_hidden=d_x)

    Dict("d_x"=>d_x, "d_x0"=>d_x0, "d_y"=>d_y, "d_enc_state"=>d_enc_state,
    "encoder"=>encoder, "cnn"=>cnn, "out_heads"=>out_heads, "d_hidden"=>d_hidden)
end

"""
    load_model_from_def(ymlfile)
We have YAML files saved with metadata associated with the model parameters
in the /data folder. These YAML files point to the relevant parameters too, and
so suffices to construct the model and load the parameters. That is the purpose
of this function.
"""
function load_model_from_def(ymlfile::String)
    model_details = YAML.load_file(ymlfile)
    filename = model_details["filename"]
    constructor = model_details["constructor"]
    @argcheck constructor == "seq2seq.create_model"
    D = model_details["model_def"]
    D["encoder"] = Symbol(D["encoder"])
    m = create_model(D["d_x"], D["d_x0"], D["d_y"], D["d_enc_state"];
        encoder=D["encoder"], cnn=D["cnn"], out_heads=D["out_heads"], d_hidden=D["d_hidden"])
    modelutils.load!(m, filename)
    return m
end
# # Example saving model info
# YAML.write_file("saved_models/video_rnn_128_8_40_T20.yml", Dict(
#     "filename"=>"data/video_rnn_128_8_40_T20.bson",
#     "description"=>"Standard GRU model trained in open loop mode. T=80, Tenc=20.",
#     "constructor"=>"seq2seq.create_model",
#     "model_def"=>create_model_opt_dict(128, 8, 4, 40, cnn=true)
# ))
end

max_stars_repo_name: ornithos/mtds-dblpend

max_issues_repo_name: ornithos/mtds-dblpend

max_forks_repo_name: ornithos/mtds-dblpend


text:
###############################################################################
#
# Choose text font and rendering method
#
###############################################################################
"""
    use_sans_tex()

Set default font to sans and use latex wrapper for texts
"""
function use_sans_tex()
    use_sans();
    rc("text"      , usetex=true);
    rc("text.latex",
       preamble = "\\usepackage{amsmath,DejaVuSans,sansmath,upgreek}" *
                  "\\sansmath");

    return nothing
end

"""
    use_sans_text()

Set default font to sans and use normal text wrapper for texts
"""
function use_sans_text()
    use_sans();
    rc("text", usetex=false);

    return nothing
end

"""
    use_serif_tex()

Set default font to serif and use latex wrapper for texts
"""
function use_serif_tex()
    use_serif();
    rc("text"      , usetex=true);
    rc("text.latex", preamble="\\usepackage{amsmath,newpxmath,newpxtext}");

    return nothing
end

"""
    use_serif_text()

Set default font to serif and use normal text wrapper for texts
"""
function use_serif_text()
    use_serif();
    rc("text", usetex=false);

    return nothing
end


text:
### Convert CO2 concentration percentage to compressional wavespeed by Patchy saturation model

## Ziyi Yin followed the example in FwiFlow software, developed by Dongzhuo Li and Kailai Xu, from the publication https://doi.org/10.1029/2019WR027032

using DrWatson
@quickactivate "CompassTimeLapseCCS"

using PyPlot
using Random, Images, JLD2, LinearAlgebra
using JOLI, Statistics, FFTW
using JUDI

JLD2.@load "../model/CompassSlice.jld2" v rho d
JLD2.@load "../data/CompassCO2.jld2" conc

n = (1590, 205)
o = (0f0, 0f0)
d = (10f0, 10f0)

function down_sample(temp1)
    ## downsample the model to 10m grid spacing w/ low-pass filter
    fft1 = fftshift(fft(temp1))
    mask = zeros(size(fft1))
    taper_length = 8
    half_taper_length = Int(taper_length/2)
    taper = (1f0 .+ sin.((Float32(pi)*(0:taper_length-1))/(taper_length - 1) .- Float32(pi)/2f0))/2f0
    mask[:,171-102:171+102] .= 1f0
    mask[:,171-102-taper_length+1+half_taper_length:171-102+half_taper_length] .= reshape(taper,1,taper_length)
    mask[:,171+102-half_taper_length:171+102+taper_length-1-half_taper_length] .= reshape(reverse(taper),1,taper_length)
    fft2 = mask .* fft1
    temp2 = real.(ifft(ifftshift(fft2)))
    temp2[:,1:35] .= temp1[:,1:35]  # fix water layer
    v = imresize(temp2,n)
    v[v.<=minimum(temp1)] .= minimum(temp1)
    v[v.>=maximum(temp1)] .= maximum(temp1)
    return v
end

v = down_sample(deepcopy(v))
rho = down_sample(deepcopy(rho))

idx_wb = find_water_bottom(v.-minimum(v))[1]

extentx = (n[1]-1)*d[1]
extentz = (n[2]-1)*d[2]

perm = 5f0*(v.-minimum(v))/(maximum(v)-minimum(v))

for i = 1:n[1]
  perm[i,1:idx_wb] .= 1f-4
  unconformaty = findall(v[i,:].>=3500f0)[1]
  perm[i,idx_wb+1:unconformaty] = perm[i,idx_wb+1:unconformaty] .+ 15f0
  perm[i,unconformaty+1:unconformaty+5] = perm[i,unconformaty+1:unconformaty+5]*1f-3
  perm[i,unconformaty+6:end] = perm[i,unconformaty+6:end] .+ 170f0
end

function poro_to_perm(phi)
    # Kozeny-carman equation
  c = (1.527/0.0314)^2f0
  perm = c*phi.^3f0./(1f0.-phi).^2f0
  return perm
end

phi = (perm*(1-0.22)^2/(1.527/0.0314)^2f0).^(1/3f0) # a proxy

vp = convert(Array{Float32},v)
vs = vp ./ sqrt(3f0)

num_vintage = 5
sw = zeros(Float32,num_vintage,n[1],n[2])
for i = 1:num_vintage
	sw[i,:,:] = imresize(conc[i,:,:]',n)
end

function Patchy(sw, vp, vs, rho, phi; bulk_min = 36.6f9, bulk_fl1 = 2.735f9, bulk_fl2 = 0.125f9, ρw = 7.766f2, ρo = 1.053f3)
    bulk_sat1 = rho .* (vp.^2f0 - 4f0/3f0 .* vs.^2f0) * 1f3
    shear_sat1 = rho .* (vs.^2f0) * 1f3

    bulk_min = zeros(Float32,size(bulk_sat1))

    bulk_min[findall(v.>=3350)] .= 50f9   # mineral bulk moduli
    bulk_min[findall(v.<3350)] .= 2.5f9   # mineral bulk moduli
    patch_temp = bulk_sat1 ./(bulk_min .- bulk_sat1) -
    bulk_fl1 ./ phi ./ (bulk_min .- bulk_fl1) +
    bulk_fl2 ./ phi ./ (bulk_min .- bulk_fl2)

    bulk_sat2 = bulk_min./(1f0./patch_temp .+ 1f0)

    bulk_new = 1f0./( (1f0.-sw)./(bulk_sat1+4f0/3f0*shear_sat1)
    + sw./(bulk_sat2+4f0/3f0*shear_sat1) ) - 4f0/3f0*shear_sat1

	bulk_new[:,1:idx_wb] = bulk_sat1[:,1:idx_wb]

    rho_new = rho + phi .* sw * (ρw - ρo) / 1f3

    Vp_new = sqrt.((bulk_new+4f0/3f0*shear_sat1)./rho_new/1f3)
    Vs_new = sqrt.((shear_sat1)./rho_new/1f3)

    return Vp_new, Vs_new, rho_new
end

vp_stack = [Patchy(sw[i,:,:], vp, vs, rho, phi)[1] for i = 1:num_vintage]
m_stack = [(1f3./vp_stack[i]).^2f0 for i = 1:num_vintage]

m0_stack = [convert(Array{Float32,2},imfilter(m_stack[i], Kernel.gaussian(7))) for i=1:num_vintage]

for i = 1:num_vintage
    m0_stack[i][:,1:idx_wb] = m_stack[i][:,1:idx_wb]
end

dm_stack = [vec(m_stack[i]-m0_stack[i]) for i = 1:num_vintage]

JLD2.@save "../model/CompassTimeLapse.jld2" num_vintage n d o idx_wb m0_stack dm_stack m_stack


text:
# Color mapping

# import .ImPlot_GetColormapSize

# function SetColormap(colormap, samples::Integer = 1)
#     ImPlot_SetColormapPlotColormap(colormap, Cint(samples))
# end

# function SetColormap(colors::ImVec4, num_colors::Integer)
#     ImPlot_SetColormapVec4Ptr(colors, Cint(num_colors))
# end

# function GetColormapColor(index::Integer)
#     out = Ref{ImVec4}()
#     ImPlot_GetColormapColor(out, Cint(index))
#     return out[]
# end

# function LerpColormap(t::Integer)
#     out = Ref{ImVec4}()
#     ImPlot_LerpColormap(out, Cint(t))
#     return out[]
# end

# function PushColormap(colormap)
#     ImPlot_PushColormapPlotColormap(colormap)
# end

# function PushColormap(colormap::Vector{ImVec4}, size::Integer)
#     ImPlot_PushColormapVec4Ptr(colormap, size)
# end

function PopColormap(count::Integer = 1)
    ImPlot_PopColormap(count)
end

function PushColormap(colormap::ImPlot.ImPlotColormap_)
    PushColormap_PlotColormap(colormap)
end
function PushColormap(colormap::AbstractString)
    PushColormap_Str(colormap)
end


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

""" The semantic enrichment algorithm for dataflow graphs.
"""
module SemanticEnrichment
export to_semantic_graph

using Base.Iterators: product
using Nullables

using LightGraphs, MetaGraphs

using Catlab.WiringDiagrams
using ..Doctrine
using ..Ontology
using ..RawFlowGraphs

""" Convert a raw flow graph into a semantic flow graph.
"""
function to_semantic_graph(db::OntologyDB, raw::WiringDiagram)::WiringDiagram
  sem = WiringDiagram(to_semantic_ports(db, input_ports(raw)),
                      to_semantic_ports(db, output_ports(raw)))

  # Add boxes.
  to_substitute = Int[]
  for v in box_ids(raw)
    raw_box = box(raw, v)
    sem_box = to_semantic_graph(db, raw_box)
    @assert add_box!(sem, sem_box) == v
    if isa(raw_box, Box) && isa(sem_box, WiringDiagram)
      # If the raw box is atomic but the semantic box is a wiring diagram,
      # we should expand it by substitution. We defer the substitution until
      # later so that the box IDs on the wires remain valid.
      push!(to_substitute, v)
    end
  end

  # Add wires.
  # FIXME: If a raw box expands to a semantic box that is missing a port with
  # incoming or outoging wires in the raw graph (e.g., when an unannotated
  # keyword argument is passed), this logic will fail. We should detect this
  # situation and either raise an informative error or discard the wire.
  add_wires!(sem, wires(raw))

  # Perform deferred substitutions (see above).
  substitute!(sem, to_substitute)

  # Simplify the diagram by collapsing adjacent unannotated boxes.
  collapse_unannotated_boxes!(sem)

  return sem
end

function to_semantic_graph(db::OntologyDB, raw_box::Box{RawNode})::AbstractBox
  if isnull(raw_box.value.annotation)
    Box(nothing,
        to_semantic_ports(db, input_ports(raw_box)),
        to_semantic_ports(db, output_ports(raw_box)))
  else
    expand_annotated_box(db, raw_box)
  end
end

function to_semantic_ports(db::OntologyDB, ports::Vector{RawPort})
  [ isnull(port.annotation) ? nothing : expand_annotated_port(db, port)
    for port in ports ]
end

""" Expand a single annotated box from a raw flow graph.
"""
function expand_annotated_box(db::OntologyDB, raw_box::Box{RawNode})::WiringDiagram
  expand_annotated_box(db, raw_box, Val{raw_box.value.annotation_kind})
end

function expand_annotated_box(db::OntologyDB, raw_box::Box{RawNode},
                              ::Type{Val{FunctionAnnotation}})
  inputs = input_ports(raw_box)
  outputs = output_ports(raw_box)
  note = load_annotation(db, get(raw_box.value.annotation))::HomAnnotation
  f = WiringDiagram(inputs, outputs)
  v = add_box!(f, to_wiring_diagram(note.definition))
  add_wires!(f, ((input_id(f), i) => (v, get(port.annotation_index))
                 for (i, port) in enumerate(inputs)
                 if !isnull(port.annotation_index)))
  add_wires!(f, ((v, get(port.annotation_index)) => (output_id(f), i)
                 for (i, port) in enumerate(outputs)
                 if !isnull(port.annotation_index)))
  substitute!(f, v)
  return f
end

function expand_annotated_box(db::OntologyDB, raw_box::Box{RawNode},
                              ::Type{Val{ConstructAnnotation}})
  note = load_annotation(db, get(raw_box.value.annotation))::ObAnnotation
  to_wiring_diagram(construct(note.definition))
end

function expand_annotated_box(db::OntologyDB, raw_box::Box{RawNode},
                              ::Type{Val{SlotAnnotation}})
  note = load_annotation(db, get(raw_box.value.annotation))::ObAnnotation
  index = get(raw_box.value.annotation_index)
  to_wiring_diagram(note.slots[index])
end

""" Expand a single annotated port from a raw flow graph.
"""
function expand_annotated_port(db::OntologyDB, raw_port::RawPort)::Monocl.Ob
  note = load_annotation(db, get(raw_port.annotation))::ObAnnotation
  note.definition
end

""" Collapse adjacent unannotated boxes into single boxes.
"""
function collapse_unannotated_boxes!(diagram::WiringDiagram)
  # Find maximal groups of unannotated boxes.
  unannotated = filter(box_ids(diagram)) do v
    isnothing(box(diagram,v).value)
  end
  groups = group_blank_vertices(graph(diagram), unannotated)

  # Encapsulate the groups, including groups of size 1 because encapsulation
  # will simplify the ports.
  encapsulate!(diagram, groups, nothing)
  return diagram
end

""" Group adjacent blank vertices of a directed graph.
"""
function group_blank_vertices(graph::SimpleDiGraph, blank::Vector{Int})::Vector{Vector{Int}}
  # Create transitive closure of graph.
  closure = transitiveclosure(graph)
  has_path(u::Int, v::Int) = has_edge(closure, u, v)
  ancestors(v::Int) = inneighbors(closure, v)
  descendants(v::Int) = outneighbors(closure, v)

  # Initialize groups as singletons.
  graph = MetaDiGraph(graph)
  for v in blank
    set_prop!(graph, v, :vertices, [v])
  end
  get_group(v::Int) = get_prop(graph, v, :vertices)
  is_blank(v::Int) = has_prop(graph, v, :vertices)
  not_blank(v::Int) = !is_blank(v)

  # Definition: Two adjacent blank vertices are mergeable if their merger
  # does not introduce any new dependencies between non-blank vertices.
  # That is, if there is no directed path from non-blank vertex `v1` to
  # non-blank vertex `v2` before merging, there will not be one afterwards.
  #
  # This criterion holds iff it's *not* the case that
  #   there exists a non-blank ancestor of the child that isn't an ancestor of
  #   the parent and a non-blank descendant of the parent isn't a descendant of
  #   the child
  # (because if there was, merging would create a new dependency) iff
  #   each non-blank ancestor of the child is also an ancestor of the parent, or
  #   each non-blank descendant of the parent is also a descendent of the child.
  function is_mergable(u::Int, v::Int)
    (is_blank(u) && is_blank(v)) || return false
    parent, child = if has_edge(graph, u, v); (u, v)
      elseif has_edge(graph, v, u); (v, u)
      else return false end
    all(has_path(v, parent) for v in filter(not_blank, ancestors(child))) ||
      all(has_path(child, v) for v in filter(not_blank, descendants(parent)))
  end
  function merge_blank!(u::Int, v::Int)
    append!(get_group(min(u,v)), get_group(max(u,v)))
    merge_vertices_directed!(graph, [u,v])
    merge_vertices_directed!(closure, [u,v])
  end

  # Merge pairs of mergable vertices until there are no more mergable pairs.
  # The loop maintains the invariant that no two vertices less than the current
  # one are mergable.
  v = 1
  while v <= nv(graph)
    merged = false
    for u in all_neighbors(graph, v)
      if u > v && is_mergable(u, v)
        merge_blank!(u, v)
        merged = true
        break
      end
    end
    v += !merged
  end

  Vector{Int}[ get_group(v) for v in vertices(graph) if is_blank(v) ]
end

""" Merge the vertices into a single vertex, preserving edges.

Note: LightGraphs.merge_vertices! only supports undirected graphs.
"""
function merge_vertices_directed!(graph::AbstractGraph, vs::Vector{Int})
  @assert is_directed(graph)
  vs = sort(vs, rev=true)
  v0 = vs[end]
  for v in vs[1:end-1]
    for u in inneighbors(graph, v)
      if !(u in vs)
        add_edge!(graph, u, v0)
      end
    end
    for u in outneighbors(graph, v)
      if !(u in vs)
        add_edge!(graph, v0, u)
      end
    end
    rem_vertex!(graph, v)
  end
  v0
end

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
# example1.jl
#
# Vin = 20V
# Iout = 5A
# 1:1 turns ratio
# 4 layer pcb, 1oz copper
# Use PlanerTrnasformer.jl to find an efficient design.

using PlanerTransformer

# define pcb
trace_edge_gap = 0.16e-3
trace_trace_gap = 0.13e-3
outer = copper_weight_to_meters(1.0)
inner = copper_weight_to_meters(1.0)
number_of_layers = 4
dielectric = (1.6e-3-2*outer-2*inner)/(number_of_layers-1)

stackup = Stackup([copper,fr408,copper,fr408,copper,fr408,copper],
                  [outer,dielectric,inner,dielectric,inner,dielectric,outer])
pcb = PCB_Specification(trace_edge_gap,
                        trace_trace_gap,
                        stackup)

# core geometries to analyze
my_geometry_list = ["e14_plt","e18_plt","e22_plt","e32_plt","e58_plt"]
my_cores = getindex.(core_geometry_dict,my_geometry_list)

# ferrites to analyze
my_ferrite_list = ["3f4","3f45","3f5","4f1"]
my_ferrites = getindex.(ferrite_dict,my_ferrite_list)

# 1 to 10 turns per layer
turns_per_layer = 1:10
my_windings = [windings(pcb,x,y,y,prisec"P-S-S-P",false,false) for x in my_cores, y in turns_per_layer]

# create transformers for all conbinations
my_transformers = [Transformer(f,w) for f in my_ferrites, w in my_windings]

v_in = 20.0
i_out = -5.0 # power flows out of secondary
frequency = 10e6
my_analysis = transformer_power_analysis.(my_transformers, PushPull(v_in, i_out, frequency))

(pmin,index) = findmin(total_power.(my_analysis))
best_analysis = my_analysis[index]
(ferrite_name(best_analysis), core_name(best_analysis), turns(best_analysis))
e = efficiency(best_analysis)
l = leakage_inductance(best_analysis)
c = capacitance(best_analysis)

# what is the efficiency of the smallest core
my_analysis_e14plt = transformer_power_analysis.(my_transformers[:,1,:], PushPull(v_in, i_out, frequency))

(pmin,index) = findmin(total_power.(my_analysis_e14plt))
best_analysis_e14plt = my_analysis_e14plt[index]
(ferrite_name(best_analysis_e14plt), core_name(best_analysis_e14plt), turns(best_analysis_e14plt))
e = efficiency(best_analysis_e14plt)
l = leakage_inductance(best_analysis_e14plt)
c = capacitance(best_analysis_e14plt)


text:
as(d::EqualMix, _data) = as(d.components[1], _data)

# StudentT(ν, μ = 0.0, σ = 1.0) = LocationScale(μ, σ, TDist(ν))

as(d::Dists.Dirichlet, _data) = UnitSimplex(length(d.alpha))


text:
# Author:
# Jürgen Vorloeper, Hochschule Ruhr West, 2018

import GLPK

"""
 simple interface to glpk
 solve f*x->max under constraints Ax<=b, x>=0,
 with A \in R^(n\times m), b\in R^n, f,lb,ub \in R^m

 #### Fields
 - 'f::Vector{Real}': xx
 - 'A::Matrix{Real}': xx
 - 'b::Vector{Real}': xx
"""
function glpk(f::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64})
#function glpk(f::Array{Real,1}, A::Array{Real,2}, b::Array{Real,1})
#function glpk{T}(f::AbstractArray{T,1}, A::AbstractArray{T,2}, b::AbstractArray{T,1})
#function glpk(f::Array{Float64,1}, A::Array{Float64,2}, b::Array{Float64,1})

	eflag = 0

	n = size(A, 1)
	m = size(A, 2)

	lp = GLPK.Prob()
	GLPK.set_prob_name(lp, "sample")
	GLPK.set_obj_dir(lp, GLPK.MAX)

	GLPK.add_rows(lp, n)

	for i=1:n
		GLPK.set_row_bnds(lp, i, GLPK.UP, 0.0, b[i])
	end

	GLPK.add_cols(lp, m)

	for j=1:m
		GLPK.set_col_bnds(lp, j, GLPK.LO, 0.0, 0.0)
		GLPK.set_obj_coef(lp, j, f[j])
	end

	GLPK.load_matrix(lp, sparse(A))
	GLPK.simplex(lp)

	z = GLPK.get_obj_val(lp)

	x = Array{Float64}(m)

	for j=1:m
		x[j] = GLPK.get_col_prim(lp, j)
    end

	return x, z, eflag

end

# simple interface to glpk
# solve LP Ax (<=,==,>=) b, lb<= x <= ub, f*x->max
# with A \in R^(n\times m), b\in R^n, f,lb,ub \in R^m
# with A[i,:] <= b[i], if ctype[i] == "U"
# with A[i,:] == b[i], if ctype[i] == "S",
# with A[i,:] >= b[i], if ctype[i] == "L"
#
# with x[i] integer,    if vartype[i] == "I"
# with x[i] continuous, if vartype[i] == "C"
# with x[i] binary,     if vartype[i] == "B"
function glpk(f::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}, lb::Vector{Float64}, ub::Vector{Float64}, ctype::String, vartype::String)

	eflag = 0

	mip = false

	n = size(A, 1)
	m = size(A, 2)

	lp = GLPK.Prob()
	GLPK.set_prob_name(lp, "sample")
	GLPK.set_obj_dir(lp, GLPK.MAX)

	GLPK.add_rows(lp, n)

	for i = 1:n
		if ctype[i] == 'U'
			GLPK.set_row_bnds(lp, i, GLPK.UP, 0.0, b[i])
		end

		if ctype[i] == 'L'
			GLPK.set_row_bnds(lp, i, GLPK.LO, b[i], 0.0)
		end

		if ctype[i] == 'S'
			GLPK.set_row_bnds(lp, i, GLPK.FX, b[i], b[i])
		end
	end

	GLPK.add_cols(lp, m)

	for j = 1:m
		if abs(lb[j]-ub[j])<1e-12
			GLPK.set_col_bnds(lp, j, GLPK.FX, lb[j], ub[j])
		else
			GLPK.set_col_bnds(lp, j, GLPK.DB, lb[j], ub[j])
		end
		GLPK.set_obj_coef(lp, j, f[j])

		if vartype[j] == 'C'
			GLPK.set_col_kind(lp,j, GLPK.CV)
		end

		if vartype[j] == 'I'
			mip = true
			GLPK.set_col_kind(lp, j, GLPK.IV)
		end

		if vartype[j] == 'B'
			mip = true
			GLPK.set_col_kind(lp, j, GLPK.BV)
		end
	end

	GLPK.load_matrix(lp, sparse(A))
	v1 = GLPK.simplex(lp)
	if v1>0
		println("*** glpk: No solution found!")
		eflag = -1
	end

	if mip==true
		v2 = GLPK.intopt(lp)

		if v2>0
			println("*** glpk: No solution found!")
			eflag = -1
		end
	end

	if mip==true
		z = GLPK.mip_obj_val(lp)
	else
		z = GLPK.get_obj_val(lp)
	end

	x = Array{Float64}(m)

	if mip==true
		for j=1:m
			x[j] = GLPK.mip_col_val(lp, j)
		end
	else
		for j=1:m
			x[j] = GLPK.get_col_prim(lp, j)
		end
    end

	return x, z, eflag
end

# simple interface to glpk
# solve LP Ax <= b, lb<= x <= ub, f*x->max
# with A \in R^(n\times m), b\in R^n, f,lb,ub \in R^m
function glpk(f::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}, lb::Vector{Float64}, ub::Vector{Float64})

	n = size(A, 1)
	m = size(A, 2)

	ctype = repeat("U", n)
	vartype = repeat("C",m)

	return glpk(f, A, b, lb, ub, ctype, vartype)
end

max_stars_repo_name: j-vorloeper/HRW

max_issues_repo_name: j-vorloeper/HRW

max_forks_repo_name: j-vorloeper/HRW


text:
# wraps variables found in EIGEN.INC for use in Julia

const NEIGX = 300
const NFREQX = 120
const NZX = 18*IIX + NGVX + NGPX + NGFX
const NDFFX = NFLPX + NENGX + 2*NYGUSX-1

struct EIGEN_s
    LROM::Array{Bool,1}
    NVTOT::Array{Int32,1}
    NEIGEN::Array{Int32,1}
    NEIDAT::Array{Int32,1}
    IPEIG::OffsetArrays.OffsetArray{Int32,1,Array{Int32,1}}
    IPDAT::OffsetArrays.OffsetArray{Int32,1,Array{Int32,1}}
    IPFRQ::Array{Int32,1}
    NFREQ::Array{Int32,1}
    NFDAT::Array{Int32,1}
    NROM::Array{Int32,1}
    NROMC::Array{Int32,1}
    KROM::Array{Int32,1}
    ZESEED::Array{ComplexF64,1}
    ZLAMDA::OffsetArrays.OffsetArray{ComplexF64,1,Array{ComplexF64,1}}
    ZLADAT::OffsetArrays.OffsetArray{ComplexF64,1,Array{ComplexF64,1}}
    ZQEIG::OffsetArrays.OffsetArray{ComplexF64,2,Array{ComplexF64,2}}
    ZLEIG::OffsetArrays.OffsetArray{ComplexF64,2,Array{ComplexF64,2}}
    ZWB::Array{ComplexF64,2}
    ZWBC::Array{ComplexF64,2}
    ZROML::Array{ComplexF64,2}
    ZROMCL::Array{ComplexF64,2}
    ZROM::Array{ComplexF64,1}
    ZROMC::Array{ComplexF64,1}
    ZQ::Array{ComplexF64,2}
    ZPARAM::Array{ComplexF64,1}
    ZQSENS::Array{ComplexF64,2}
    ZQSENSF::Array{ComplexF64,3}
    ZAN::Array{ComplexF64,1}
    ZQJOIN::Array{ComplexF64,2}
    ZAFORCE::Array{ComplexF64,1}
    ZRFORCE::Array{ComplexF64,1}
    ZGFORCE::Array{ComplexF64,1}
    ZTFORCE::Array{ComplexF64,1}
    ZAMOMNT::Array{ComplexF64,1}
    ZRMOMNT::Array{ComplexF64,1}
    ZGMOMNT::Array{ComplexF64,1}
    ZTMOMNT::Array{ComplexF64,1}
    ZFLIFT::Array{ComplexF64,1}
    ZFDRAG::Array{ComplexF64,1}
    ZFUNS::OffsetArrays.OffsetArray{ComplexF64,3,Array{ComplexF64,3}}
    ZDATS::OffsetArrays.OffsetArray{ComplexF64,3,Array{ComplexF64,3}}
    ZQSFUNS::OffsetArrays.OffsetArray{ComplexF64,4,Array{ComplexF64,4}}
    ZQSDATS::OffsetArrays.OffsetArray{ComplexF64,4,Array{ComplexF64,4}}
    AQEIG::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}
    EPHASE::Array{Float64,1}
    FREQ::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}
    FDAT::OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}}
    function EIGEN_s()
        as_eigl = cglobal((:as_eigl_, libaswing), Int32); as_eigl_idx = as_eigl
        LROM = unsafe_wrap(Array, as_eigl_idx, NEIGX)

        as_eigi = cglobal((:as_eigi_, libaswing), Int32); as_eigi_idx = as_eigi
        NVTOT = unsafe_wrap(Array, as_eigi_idx, 1); as_eigi_idx += sizeof(Int32)
        NEIGEN = unsafe_wrap(Array, as_eigi_idx, 1); as_eigi_idx += sizeof(Int32)
        NEIDAT = unsafe_wrap(Array, as_eigi_idx, 1); as_eigi_idx += sizeof(Int32)
        ipeig = unsafe_wrap(Array, as_eigi_idx, NEIGX+1); as_eigi_idx += sizeof(Int32)*(NEIGX+1)
        IPEIG = OffsetArrays.OffsetArray(ipeig, 0:NEIGX)
        ipdat = unsafe_wrap(Array, as_eigi_idx, NEIGX+1); as_eigi_idx += sizeof(Int32)*(NEIGX+1)
        IPDAT = OffsetArrays.OffsetArray(ipdat, 0:NEIGX)
        IPFRQ = unsafe_wrap(Array, as_eigi_idx, 1); as_eigi_idx += sizeof(Int32)
        NFREQ = unsafe_wrap(Array, as_eigi_idx, 1); as_eigi_idx += sizeof(Int32)
        NFDAT = unsafe_wrap(Array, as_eigi_idx, NDFFX); as_eigi_idx += sizeof(Int32)*NDFFX
        NROM = unsafe_wrap(Array, as_eigi_idx, 1); as_eigi_idx += sizeof(Int32)
        NROMC = unsafe_wrap(Array, as_eigi_idx, 1); as_eigi_idx += sizeof(Int32)
        KROM = unsafe_wrap(Array, as_eigi_idx, NEIGX)

        as_eigz = cglobal((:as_eigz_, libaswing), ComplexF64); as_eigz_idx = as_eigz
        ZESEED = unsafe_wrap(Array, as_eigz_idx, 1); as_eigz_idx += sizeof(ComplexF64)
        zlamda = unsafe_wrap(Array, as_eigz_idx, NEIGX+1); as_eigz_idx += sizeof(ComplexF64)*(NEIGX+1)
        ZLAMDA = OffsetArrays.OffsetArray(zlamda, 0:NEIGX)
        zladat = unsafe_wrap(Array, as_eigz_idx, NEIGX+1); as_eigz_idx += sizeof(ComplexF64)*(NEIGX+1)
        ZLADAT = OffsetArrays.OffsetArray(zladat, 0:NEIGX)
        zqeig = unsafe_wrap(Array, as_eigz_idx, (NZX, NEIGX+1)); as_eigz_idx += sizeof(ComplexF64)*NZX*(NEIGX+1)
        ZQEIG = OffsetArrays.OffsetArray(zqeig, 1:NZX, 0:NEIGX)
        zleig = unsafe_wrap(Array, as_eigz_idx, (NZX, NEIGX+1)); as_eigz_idx += sizeof(ComplexF64)*NZX*(NEIGX+1)
        ZLEIG = OffsetArrays.OffsetArray(zleig, 1:NZX, 0:NEIGX)
        ZWB = unsafe_wrap(Array, as_eigz_idx, (NEIGX, NFRPX)); as_eigz_idx += sizeof(ComplexF64)*NEIGX*NFRPX
        ZWBC = unsafe_wrap(Array, as_eigz_idx, (NEIGX, NFRPX)); as_eigz_idx += sizeof(ComplexF64)*NEIGX*NFRPX
        ZROML = unsafe_wrap(Array, as_eigz_idx, (NEIGX, NFRPX)); as_eigz_idx += sizeof(ComplexF64)*NEIGX*NFRPX
        ZROMCL = unsafe_wrap(Array, as_eigz_idx, (NEIGX, NFRPX)); as_eigz_idx += sizeof(ComplexF64)*NEIGX*NFRPX
        ZROM = unsafe_wrap(Array, as_eigz_idx, NEIGX); as_eigz_idx += sizeof(ComplexF64)*NEIGX
        ZROMC = unsafe_wrap(Array, as_eigz_idx, NEIGX); as_eigz_idx += sizeof(ComplexF64)*NEIGX
        ZQ = unsafe_wrap(Array, as_eigz_idx, (18, IIX)); as_eigz_idx += sizeof(ComplexF64)*18*IIX
        ZPARAM = unsafe_wrap(Array, as_eigz_idx, KPTOT); as_eigz_idx += sizeof(ComplexF64)*KPTOT
        ZQSENS = unsafe_wrap(Array, as_eigz_idx, (KSTOT, NSENX)); as_eigz_idx += sizeof(ComplexF64)*KSTOT*NSENX
        ZQSENSF = unsafe_wrap(Array, as_eigz_idx, (KSTOT, NSENX, NFRPX)); as_eigz_idx += sizeof(ComplexF64)*KSTOT*NSENX*NFRPX
        ZAN = unsafe_wrap(Array, as_eigz_idx, NNX); as_eigz_idx += sizeof(ComplexF64)*NNX
        ZQJOIN = unsafe_wrap(Array, as_eigz_idx, (12, NJX)); as_eigz_idx += sizeof(ComplexF64)*12*NJX
        ZAFORCE = unsafe_wrap(Array, as_eigz_idx, 3); as_eigz_idx += sizeof(ComplexF64)*3
        ZRFORCE = unsafe_wrap(Array, as_eigz_idx, 3); as_eigz_idx += sizeof(ComplexF64)*3
        ZGFORCE = unsafe_wrap(Array, as_eigz_idx, 3); as_eigz_idx += sizeof(ComplexF64)*3
        ZTFORCE = unsafe_wrap(Array, as_eigz_idx, 3); as_eigz_idx += sizeof(ComplexF64)*3
        ZAMOMNT = unsafe_wrap(Array, as_eigz_idx, 3); as_eigz_idx += sizeof(ComplexF64)*3
        ZRMOMNT = unsafe_wrap(Array, as_eigz_idx, 3); as_eigz_idx += sizeof(ComplexF64)*3
        ZGMOMNT = unsafe_wrap(Array, as_eigz_idx, 3); as_eigz_idx += sizeof(ComplexF64)*3
        ZTMOMNT = unsafe_wrap(Array, as_eigz_idx, 3); as_eigz_idx += sizeof(ComplexF64)*3
        ZFLIFT = unsafe_wrap(Array, as_eigz_idx, 1); as_eigz_idx += sizeof(ComplexF64)
        ZFDRAG = unsafe_wrap(Array, as_eigz_idx, 1); as_eigz_idx += sizeof(ComplexF64)
        zfuns = unsafe_wrap(Array, as_eigz_idx, (NFREQX+1, KPTOT, NDFFX)); as_eigz_idx += sizeof(ComplexF64)*(NFREQX+1)*KPTOT*NDFFX
        ZFUNS = OffsetArrays.OffsetArray(zfuns, 0:NFREQX, 1:KPTOT, 1:NDFFX)
        zdats = unsafe_wrap(Array, as_eigz_idx, (NFREQX+1, KPTOT, NDFFX)); as_eigz_idx += sizeof(ComplexF64)*(NFREQX+1)*KPTOT*NDFFX
        ZDATS = OffsetArrays.OffsetArray(zdats, 0:NFREQX, 1:KPTOT, 1:NDFFX)
        zqsfuns = unsafe_wrap(Array, as_eigz_idx, (NFREQX+1, KSTOT, NSENX, NDFFX)); as_eigz_idx += sizeof(ComplexF64)*(NFREQX+1)*KSTOT*NSENX*NDFFX
        ZQSFUNS = OffsetArrays.OffsetArray(zqsfuns, 0:NFREQX, 1:KSTOT, 1:NSENX, 1:NDFFX)
        zqsdats = unsafe_wrap(Array, as_eigz_idx, (NFREQX+1, KSTOT, NSENX, NDFFX)); as_eigz_idx += sizeof(ComplexF64)*(NFREQX+1)*KSTOT*NSENX*NDFFX
        ZQSDATS = OffsetArrays.OffsetArray(zqsdats, 0:NFREQX, 1:KSTOT, 1:NSENX, 1:NDFFX)

        as_eigr = cglobal((:as_eigr_, libaswing), Float64); as_eigr_idx = as_eigr
        aqeig = unsafe_wrap(Array, as_eigr_idx, (NEIGX+1)); as_eigr_idx += sizeof(Float64)*(NEIGX+1)
        AQEIG = OffsetArrays.OffsetArray(aqeig, 0:NEIGX)
        EPHASE = unsafe_wrap(Array, as_eigr_idx, 1); as_eigr_idx += sizeof(Float64)
        freq = unsafe_wrap(Array, as_eigr_idx, (NFREQX+1)); as_eigr_idx += sizeof(Float64)*(NFREQX+1)
        FREQ = OffsetArrays.OffsetArray(freq, 0:NFREQX)
        fdat = unsafe_wrap(Array, as_eigr_idx, (NFREQX+1, NDFFX))
        FDAT = OffsetArrays.OffsetArray(fdat, 0:NFREQX, 1:NDFFX)

        return new(LROM, NVTOT, NEIGEN, NEIDAT, IPEIG,
            IPDAT, IPFRQ, NFREQ, NFDAT, NROM, NROMC, KROM, ZESEED, ZLAMDA,
            ZLADAT, ZQEIG, ZLEIG, ZWB, ZWBC, ZROML, ZROMCL, ZROM, ZROMC, ZQ,
            ZPARAM, ZQSENS, ZQSENSF, ZAN, ZQJOIN, ZAFORCE, ZRFORCE, ZGFORCE,
            ZTFORCE, ZAMOMNT, ZRMOMNT, ZGMOMNT, ZTMOMNT, ZFLIFT, ZFDRAG, ZFUNS,
            ZDATS, ZQSFUNS, ZQSDATS, AQEIG, EPHASE, FREQ, FDAT)
    end
end


text:
@test_approx_eq_eps(pagerank(g5)[3], 0.318, 0.001)
@test_throws ErrorException pagerank(g5, 2)
@test_throws ErrorException pagerank(g5, 0.85, 2)


text:
using ProgressMeter

function main()
    out_dir="result"
    PARAMDAT="PARAMDAT.json"
    xmax,ymax,nodes,vecAx,vecAy = read_allgrid()
    nt,dt,init_rho,init_u,init_v,init_p,specific_heat_ratio,bdcon,every_outnum,out_file_front,out_ext,restart_file,restartnum,time_integ,init_small,norm_ok,nt_lusgs = input_para(PARAMDAT)

    Qbase,cellxmax,cellymax,restartnum = set_initQbase(xmax,ymax,restart_file,init_rho,init_u,init_v,init_p,specific_heat_ratio,out_file_front,out_ext,out_dir,restartnum)

    # init Delta_Qcon_hat
    Delta_Qcon_hat = zeros(cellxmax,cellymax,4)

    # main loop
    println("threads num ")
    println(Threads.nthreads())

    prog = Progress(nt,1)
    for k in 1:nt
        next!(prog)

        volume = set_volume(nodes,cellxmax,cellymax)

        Qbase = set_boundary(Qbase,cellxmax,cellymax,vecAx,vecAy,bdcon)
        Qcon = base_to_conservative(Qbase,cellxmax,cellymax,specific_heat_ratio)

        Qcon_hat = setup_Qcon_hat(Qcon,cellxmax,cellymax,volume)

        cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus = cal_jacobi(Qbase,Qcon,cellxmax,cellymax,specific_heat_ratio,vecAx,vecAy)
        cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus = setup_cell_flux_hat(Qcon,cellxmax,cellymax,cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus)

        E_hat,F_hat = FVS(cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus,cellxmax,cellymax)

        RHS = setup_RHS(cellxmax,cellymax,E_hat,F_hat,Qcon_hat)

        evalnum = k+restartnum
        if time_integ == "1"
            # exlicit scheme
            Qcon_hat = time_integration_explicit(dt,Qcon_hat,RHS,cellxmax,cellymax)
        elseif time_integ == "2"
            # implicit scheme
            A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, A_beta_sig, B_beta_sig = one_wave(Qbase,Qcon,cellxmax,cellymax,vecAx,vecAy,specific_heat_ratio,volume)
            Lx,Ly,Ux,Uy,D = set_LDU(dt,Qcon_hat,A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, A_beta_sig, B_beta_sig,RHS,cellxmax,cellymax,volume)

            Delta_Qcon_hat_temp = zeros(cellxmax,cellymax,4)

            for lusgs_t in 1:nt_lusgs
                nt_temp = 100
                norm2 = 0.0

                for i in 1:nt_temp
                    #println(Delta_Qcon_hat[2,2,:])
                    Delta_Qcon_hat = lusgs(dt,RHS,cellxmax,cellymax,volume,Lx,Ly,Ux,Uy,D,Delta_Qcon_hat)
                    #println(Delta_Qcon_hat[2,2,:])

                    res = set_res(Delta_Qcon_hat,Delta_Qcon_hat_temp,cellxmax,cellymax)
                    norm2 = check_converge(res,RHS,cellxmax,cellymax,init_small)

                    Delta_Qcon_hat_temp = deepcopy(Delta_Qcon_hat)

                    #throw(UndefVarError(:x))

                    #println(norm2)
                    #println(Delta_Qcon_hat[2,2,:])
                    #println(dt*RHS[2,2,:])
                    println("\n ---------------------------------- \n")
                end

                if norm2[1] < norm_ok && norm2[4] < norm_ok
                    println("\n ---------------------------------- \n")
                    println("nt : "*string(round(evalnum))*"     lusgs_t : "*string(lusgs_t))
                    println("density res:"*string(norm2[1]) * "  energy res:"*string(norm2[4]))
                    println("\n ---------------------------------- \n")
                    break
                end

            end
            #Delta_Qcon_hat = lusgs(dt,RHS,cellxmax,cellymax,volume,Lx,Ly,Ux,Uy,D,Delta_Qcon_hat)
            #Delta_Qcon_hat = lusgs(dt,RHS,cellxmax,cellymax,volume,Lx,Ly,Ux,Uy,D,Delta_Qcon_hat)

            # lusgs = update
            for i in 2:cellxmax-1
                for j in 2:cellymax-1
                    for l in 1:4
                        Qcon_hat[i,j,l] = Qcon_hat[i,j,l] + Delta_Qcon_hat[i,j,l]
                    end
                end
            end

            #println(dt*RHS[2,2,:])

            # 内部反復
            #=
            nt_inner = 10
            Qcon_hat_n = copy(Qcon_hat)
            Qcon_hat_m1 = copy(Qcon_hat)
            Delta_Qcon_hat_m = copy(Delta_Qcon_hat)

            for inner_t in 1:nt_inner
                Delta_Qcon_hat_temp = zeros(cellxmax,cellymax,4)

                for lusgs_t in 1:nt_lusgs
                    Qcon_m1 = Qhat_to_Q(Qcon_hat_m1,cellxmax,cellymax,volume)
                    Qbase_m1 = conservative_to_base(Qcon_m1,cellxmax,cellymax,specific_heat_ratio)

                    cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus = cal_jacobi(Qbase_m1,Qcon_m1,cellxmax,cellymax,specific_heat_ratio,vecAx,vecAy)
                    cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus = setup_cell_flux_hat(Qcon_m1,cellxmax,cellymax,cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus)
                    E_hat,F_hat = FVS(cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus,cellxmax,cellymax)

                    # lusgsと違うので注意
                    RHS = set_RHS_inner(cellxmax,cellymax,dt,Qcon_hat_n,Qcon_hat_m1,E_hat,F_hat)
                    Lx,Ly,Ux,Uy,Dinv = set_LDU_inner(dt,cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus,cellxmax,cellymax,volume)

                    Delta_Qcon_hat_m = inner_iteration(dt,RHS,cellxmax,cellymax,volume,Lx,Ly,Ux,Uy,Dinv,Delta_Qcon_hat_m)

                    res = set_res(Delta_Qcon_hat,Delta_Qcon_hat_temp,cellxmax,cellymax)
                    norm2 = check_converge(res,RHS,cellxmax,cellymax,init_small)
                    Delta_Qcon_hat_temp = copy(Delta_Qcon_hat_m)

                    if norm2[1] < norm_ok && norm2[4] < norm_ok
                        println("\n ---------------------------------- \n")
                        println("nt : "*string(round(inner_t)))
                        println("density res:"*string(norm2[1]) * "  energy res:"*string(norm2[4]))
                        println("\n ---------------------------------- \n")
                        break
                    end
                end

                # inner update
                for i in 2:cellxmax-1
                    for j in 2:cellymax-1
                        for l in 1:4
                            Qcon_hat_m1[i,j,l] = Qcon_hat_m1[i,j,l] + Delta_Qcon_hat_m[i,j,l]
                        end
                    end
                end
            end

            Qcon_hat = copy(Qcon_hat_m1)
            =#

            #println(Qcon_hat[:,:,1])
        else
            println(time_integ)
            println("pressure check time_integ")
            throw(UndefVarError(:x))
        end

        Qcon = Qhat_to_Q(Qcon_hat,cellxmax,cellymax,volume)
        Qbase = conservative_to_base(Qcon,cellxmax,cellymax,specific_heat_ratio)

        if round(evalnum) % every_outnum == 0
            println("nt_______________________________"*string(round(evalnum)))
            output_result(evalnum,Qbase,cellxmax,cellymax,specific_heat_ratio,out_file_front,out_ext,out_dir)
        end

        if Qcon[2,2,4] < 0
            println("pressure are minus !!")
            break
        end
        if Qcon[2,2,1] == NaN
            println(" NaN !!")
            break
        end

    end
end

# -- main --
main()


text:
module TestClockResolution
    import Benchmarks
    using Base.Test
    using Compat

    res = Benchmarks.estimate_clock_resolution(1)

    @test isa(res, UInt)
    @test 1 <= res <= 10_000

    res_2 = Benchmarks.estimate_clock_resolution()

    @test res_2 <= res
end


text:
"""
    ThresholdCompression(threshold=1) <: CompressionStrategy

[`CompressionStrategy`](@ref) that compresses a vector by threshold projection. Every entry
in the vector with a value below the threshold is either set to zero, or increased to the
threshold. The probabilty of setting it to zero is equal to `abs(value) / threshold`.
"""
struct ThresholdCompression{T} <: CompressionStrategy
    threshold::T
end
ThresholdCompression() = ThresholdCompression(1)

function compress!(t::ThresholdCompression, v)
    w = localpart(v)
    for (add, val) in pairs(w)
        prob = abs(val) / t.threshold
        if prob < 1 # projection is only necessary if abs(val) < s.threshold
            val = ifelse(prob > cRand(), t.threshold * sign(val), zero(val))
            w[add] = val
        end
    end
    return v
end

"""
    DoubleOrNothing(; threshold=1, prob=0.5) <: CompressionStrategy

[`CompressionStrategy`](@ref) that compresses a vector by a double or nothing approach. For
each value in the vector below the `threshold`, a coin is flipped. If heads, the value is
doubled, if tails, the value is removed from the vector. The fairness of the coin is
controlled with `prob` - lower values indicate lower probabilties to keep.
"""
struct DoubleOrNothing{T} <: CompressionStrategy
    threshold::T
    prob::Float64
end
DoubleOrNothing(; threshold=1, prob=0.5) = DoubleOrNothing(threshold, prob)

function compress!(d::DoubleOrNothing, v)
    w = localpart(v)
    for (add, val) in pairs(w)
        if abs(val) < d.threshold
            val = ifelse(cRand() > d.prob, zero(valtype(w)), val / d.prob)
            w[add] = val
        end
    end
    return v
end

"""
    DoubleOrNothingWithTarget(; target, threshold=1.0) <: CompressionStrategy

Like [`DoubleOrNothing`](@ref), but attempts to keep the number of elements in the vector
after compression as close to `target` as possible. Works like `DoubleOrNothing`, but the
`prob` parameter is chosen dynamically to try to reach the `target` length.
"""
struct DoubleOrNothingWithTarget{T} <: CompressionStrategy
    target::Int
    threshold::T
end
function DoubleOrNothingWithTarget(; target, threshold=1)
    return DoubleOrNothingWithTarget(target, threshold)
end

function compress!(d::DoubleOrNothingWithTarget, v)
    w = localpart(v)
    prob = d.target / length(v)
    if prob < 1
        compress!(DoubleOrNothing(d.threshold, prob), v)
    end
    return v
end

"""
    DoubleOrNothingWithThreshold(; threshold_hi=1.0, threshold_lo=1e-3, prob=0.5) <: CompressionStrategy

A combination of [`ThresholdProjection`](@ref) and [`DoubleOrNothing`](@ref). If a value is
below `threshold_lo`, [`ThresholdProjection`](@ref) is applied, and if it is below
`threshold_hi`, [`DoubleOrNothing`](@ref) is. The `prob` parameter has the same function as
in [`DoubleOrNothing`](@ref).
"""
struct DoubleOrNothingWithThreshold{T} <: CompressionStrategy
    threshold_hi::T
    threshold_lo::T
    prob::Float64
end
function DoubleOrNothingWithThreshold(; threshold_hi=1.0, threshold_lo=1e-3, prob=0.5)
    return DoubleOrNothingWithThreshold(promote(threshold_hi, threshold_lo)..., prob)
end

function compress!(d::DoubleOrNothingWithThreshold, v)
    w = localpart(v)
    for (add, val) in pairs(w)
        if abs(val) < d.threshold_lo
            # Threshold projection
            prob = abs(val) / d.threshold_lo
            w[add] = ifelse(prob > cRand(), d.threshold_lo * sign(val), zero(val))
        elseif abs(val) < d.threshold_hi
            # Double or nothing projection
            val = ifelse(cRand() > d.prob, zero(valtype(w)), val / d.prob)
            w[add] = val
        end
    end
    return v
end


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

export boundingsphere

export WelzlMTF
export WelzlPivot
export Ritter

abstract type BoundingSphereAlg end

"""
    WelzlMTF()

Welzl algorithm with move to front heuristic.
See Algorithm I in https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf.
In almost all situations it is better to use [`WelzlPivot`](@ref) instead.

## Pros
* Fast for small examples

## Cons
* Prone to numerical stability issues
"""
struct WelzlMTF <: BoundingSphereAlg end

"""
    WelzlPivot(;max_iterations=1000)

Welzl algorithm with pivoting. See Algorithm II in https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf.

## Pros
* Fast

## Cons
* In very rare cases can be numerically instable
"""
struct WelzlPivot <: BoundingSphereAlg
    max_iterations::Int
end

function WelzlPivot(;max_iterations=1000)
    WelzlPivot(max_iterations)
end

"""
    center, radius = boundingsphere(pts [, algorithm=WelzlPivot()])

Compute the smallest sphere that contains each point in `pts`.

# Arguments

* pts: A list of points. Points should be vectors with floating point entries.
* algorithm: An optional algorithm to do the computation. See names(BoundingSphere) to get
"""
function boundingsphere end

function boundingsphere!(pts, alg::WelzlMTF=WelzlMTF())
    bdry = create_boundary_device(pts, alg)
    ball, support_count = welzl!(pts, bdry, alg)
    r = radius(ball)
    c = center(ball)
    c, r
end

function boundingsphere!(pts, alg::BoundingSphereAlg)
    bdry = create_boundary_device(pts, alg)
    ball = welzl!(pts, bdry, alg)
    r = radius(ball)
    c = center(ball)
    c, r
end

@noinline function boundingsphere(pts, alg::BoundingSphereAlg=WelzlMTF())
    boundingsphere!(copy(pts), alg)
end


text:
module SimpleLife
export life_step, life_step!, random_life

"""
`_get(A,i,j,wrap=fasle)` returns `A[i,j]` or `0` if subscripts are out of bounds.
With `wrap==true` then the indices are taken modulo the number of rows/columns
as appropriate.
"""
function _get(A::Matrix{T}, i::Int, j::Int, wrap::Bool=false)::T where T<:Number
    r,c = size(A)
    if (1<=i<=r) && (1<=j<=c)
        return @inbounds A[i,j]
    end

    if !wrap
        return zero(T)
    end

    i = mod(i,r)
    j = mod(j,c)
    if i==0
        i=r
    end
    if j==0
        j=c
    end
    return @inbounds A[i,j]
end

"""
`_neighbor_count(A,i,j,wrap=false)` counts how many neighbors of entry `i,j`
of matrix `A` are nonzero.
"""
function _neighbor_count(A::Matrix{Int},i::Int,j::Int,wrap::Bool=false)::Int
    return sum(_get(A,p,q,wrap) for p=i-1:i+1 for q=j-1:j+1 if (p,q) != (i,j))
end

"""
`_new_entry(A,i,j)` computes the new entry in cell `i,j`
if the current state is given by matrix `A`.
"""
function _new_entry(A::Matrix{Int},i::Int,j::Int,wrap::Bool=false)
    nc = _neighbor_count(A,i,j,wrap)
    if A[i,j]==0
        if nc==3
            return 1
        else
            return 0
        end
    end
    # so A[i,j] is 1
    if nc==2 || nc==3
        return 1
    else
        return 0
    end
end

"""
`life_step(A,wrap=false)` gives the next generation after state `A`
in Conway's *Game of Life*. With `wrap==true` treats the board toroidally.
"""
function life_step(A::Matrix{Int},wrap::Bool=false)::Matrix{Int}
    r,c = size(A)
    B = zeros(Int,r,c)
    for i=1:r
        for j=1:c
            @inbounds B[i,j] = _new_entry(A,i,j,wrap)
        end
    end
    return B
end

"""
`life_step!(A,wrap=false)` overwrites `A` with `life_step(A,wrap)`.
This function returns `true` if there has been a change to `A`;
otherwise, returns `false` if there was no change.
"""
function life_step!(A::Matrix{Int}, wrap::Bool=false)
    r,c = size(A)
    B = life_step(A,wrap)
    change = A!=B
    for i=1:r
        for j=1:c
            @inbounds A[i,j] = B[i,j]
        end
    end
    return change
end

"""
`random_life(r,c,p=0.5)` creates a random `r`-by-`c` instance of a life board
with density `p`.

`random_life(n,p=0.5)` is equivalent to `random_life(n,n,p)`.
"""
function random_life(r::Int, c::Int, p::AbstractFloat=0.5)
    M = rand(r,c) .< p
    return Int.(M)
end

random_life(n::Int,p::AbstractFloat=0.5) = random_life(n,n,p)

"""
`_life_check(A)` determines if the integer matrix `A` only contains
zeros and ones.
"""
function _life_check(A::Matrix{Int})
    return 0 <= minimum(A) <= maximum(A) <= 1
end

include("visualize.jl")

end  # module SimpleLife


text:
export DataSetofSets,
    DataFull,
    DataRoyaltyProduce,
    DataDrillOnly,
    split_thetas,
    merge_thetas

# -------------------------------------------------
# Data structure
# -------------------------------------------------

function check_theta(d::AbstractDataStructure, theta)
    _nparm(d) == length(theta) || throw(DimensionMismatch("_nparm(d) = $(_nparm(d)) != length(theta) = $(length(theta))"))
    all(isfinite.(theta)) || throw(error("theta = $theta not finite"))
    return true
end

"check coef restrictions fit in parameter vector"
function check_coef_restr_in_parm_vec(cf, d::AbstractDataSet)
    n = _nparm(d)
    issubset(Set(cf), OneTo(n)) || throw(error("coefs $cf not in 1:$n"))
    allunique(cf) || throw(error("coefs $cf not unique"))
    return true
end

struct DataSetofSets{
    D<:Union{AbstractDataDrill,EmptyDataSet},  R<:Union{DataRoyalty,EmptyDataSet}, P<:Union{DataProduce,EmptyDataSet},
    } <: AbstractDataSetofSets
    data::Tuple{D,R,P}
    coef_link_pdxn::Vector{Int}
    coef_link_drill::Vector{Int}
    function DataSetofSets(d::D, r::R, p::P, cfp::V, cfd::V) where {D,R,P,V<:Vector{Int}}
        drp = (d, r, p)
        lengths = length.(drp)
        issubset(lengths, (0, maximum(lengths),) ) || throw(DimensionMismatch("datasets must same or 0 lengths. lengths = $lengths"))
        length(cfp) == length(cfd) || throw(DimensionMismatch())
        check_coef_restr_in_parm_vec(cfp, p)
        check_coef_restr_in_parm_vec(cfd, d)
        any(lengths .> 0 ) || throw(error("one dataset must have nonzero length"))
        return new{D,R,P}(drp, cfp, cfd)
    end
end

# some versions
const DataFull = DataSetofSets{<:AbstractDataDrill, <:DataRoyalty, <:DataProduce}
const DataRoyaltyProduce = DataSetofSets{EmptyDataSet, <:DataRoyalty, <:DataProduce}
const DataDrillOnly = DataSetofSets{<:AbstractDataDrill, EmptyDataSet, EmptyDataSet}

# If no coef_links provided
DataSetofSets(d,r,p) = DataSetofSets(d,r,p,zeros(Int,0),zeros(Int,0))

function DataSetofSets(d,r,p, cfp::Vector, cfd::Vector)
    DataSetofSets(d,r,p, map(f -> f(p), cfp), map(f -> f(d), cfd))
end

function DataSetofSets(d,r,p, ct::Tuple{Vector,Vector})
    DataSetofSets(d,r,p,ct...)
end

function DataSetofSets(d,r,p,cfl::Vector{<:Tuple})
    return DataSetofSets(d,r,p, first.(cfl), last.(cfl))
end

function SimulationDraws(data::DataSetofSets, M)
    n = num_i(data)
    k = _nparm(drill(data))
    return SimulationDraws(M,n,k)
end

function DataSetofSets(dsos::DataSetofSets, ddrill::DataDrill)
    return DataSetofSets(ddrill, royalty(dsos), produce(dsos), coef_link_pdxn(dsos), coef_link_drill(dsos))
end

# -------------------------------------------------
# methods
# -------------------------------------------------

data(d::DataSetofSets) = d.data
coef_link_drill(d::DataSetofSets) = d.coef_link_drill
coef_link_pdxn(d::DataSetofSets) = d.coef_link_pdxn
coef_links(d::DataSetofSets) = zip(coef_link_pdxn(d), coef_link_drill(d))
length(d::DataSetofSets) = 3
num_i(d::DataSetofSets) = maximum(length.(data(d)))

iterate(d::DataSetofSets, state...) = iterate(data(d), state...)
Broadcast.broadcastable(d::DataSetofSets) = data(d)

drill(  d::DataSetofSets) = d.data[1] # drill
royalty(d::DataSetofSets) = d.data[2] # royalty
produce(d::DataSetofSets) = d.data[3] # produce

_nparm(d::DataFull)      = sum(_nparm.(d)) - length(coef_links(d)) - 1
_nparm(d::DataSetofSets) = sum(_nparm.(d)) - length(coef_links(d))

# @deprecate DataDrill(  d::DataSetofSets) drill(d)
# @deprecate DataRoyalty(d::DataSetofSets) royalty(d)
# @deprecate DataProduce(d::DataSetofSets) produce(d)

# -------------------------------------------------
# Indexing into parameter vector
# -------------------------------------------------

# if dataset is empty
idx_drill(  ::EmptyDataSet) = 1:0
idx_royalty(::EmptyDataSet) = 1:0
idx_produce(::EmptyDataSet) = 1:0

# standard indexing, no overlap
idx_drill(  data::DataSetofSets) = idx_drill(drill(data))
idx_royalty(data::DataSetofSets) = last(idx_drill(data)) .+ idx_royalty(royalty(data))
idx_produce(data::DataSetofSets) = last(idx_royalty(data)) .+ idx_produce(produce(data))

# theta_ρ(data::DataSetofSets{EmptyDataSet,<:DataRoyalty}, theta) = theta[1]
theta_ρ(data::DataSetofSets{<:AbstractDataDrill}, theta) = theta[_nparm(drill(data))]
theta_ρ(data::DataSetofSets{<:EmptyDataSet}, theta) = theta[1]

function merge_thetas(thetas::NTuple{3,AbstractVector}, data::DataFull)
    length.(thetas) == _nparm.(data) ||
        throw(DimensionMismatch("len thetas = $(length.(thetas)) but nparm = $(_nparm.(data))"))
    thet_d, thet_r, thet_p = thetas

    data_p = produce(data)
    idx_p = OneTo(_nparm(data_p))
    idx_p_drops = coef_link_pdxn(data)

    thet_p_short = [thet_p[i] for i in idx_p if i ∉ idx_p_drops]

    return vcat(thet_d, thet_r[2:end], thet_p_short)
end

merge_thetas(thetas, data) = vcat(thetas...)

@deprecate theta_linked(thetas::NTuple{3,AbstractVector}, data::DataFull) merge_thetas(thetas,data)

# full datasets
function idx_royalty(data::DataFull)
    kd = _nparm(drill(data))
    return (kd-1) .+ idx_royalty(royalty(data))
end

function idx_produce(data::DataFull)
    dd, dr, dp = drill(data), royalty(data), produce(data)
    kd = _nparm(dd)
    kr = _nparm(dr)
    idx = collect((kd+kr-1) .+ idx_produce(dp))
    for (idxp, idxd) in coef_links(data)
        idx[idxp] = idxd
        idx[idxp+1:end] .-= 1
    end
    return idx
end

function theta_indexes(data::DataSetofSets)
    d = idx_drill(  data)
    r = idx_royalty(data)
    p = idx_produce(data)
    return d, r, p
end

# to ensure that we get a copy... NOT a view
theta_produce(d::DataFull, theta) = theta[idx_produce(d)]

function split_thetas(data::DataSetofSets, theta::AbstractVector)
    d = theta_drill(  data, theta)
    r = theta_royalty(data, theta)
    p = theta_produce(data, theta)
    return d, r, p
end

@deprecate thetas(data::DataSetofSets, theta::AbstractVector) split_thetas(data,theta)
@deprecate split_theta(data, theta) split_thetas(data, theta)

function coefnames(data::DataSetofSets)
    nms = coefnames.(data)
    return merge_thetas(nms, data)
end


text:
using BilevelJuMP
using Test, Gurobi, MathOptInterface, JuMP, Dualization
using MathOptFormat

const MOI  = MathOptInterface
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges
const MOIT = MathOptInterface.Test

# original problem
# min -4x -3y
# s.t.
#      2x + y <= 4
#       x +2y <= 4
#       x     >= 0
#           y >= 0
#
# sol x = y = 4/3
# obj = 28/3

# bilevel modifications
# min -4x -3y
# s.t.
# y = argmin_y y
#      2x + y <= 4
#       x +2y <= 4
#       x     >= 0
#           y >= 0
#
# sol: x = 2, y = 0
# obj_upper = -8
# obj_lower =  0

@testset "Simple LP" begin
    optimizer = Gurobi.Optimizer()
    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOI.empty!(bridged)
    @test MOI.is_empty(bridged)

    # add 10 variables - only diagonal is relevant
    X = MOI.add_variables(bridged, 2)

    c1 = MOIU.normalize_and_add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(2.0, X[1]),
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    c2 = MOIU.normalize_and_add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1]),
            MOI.ScalarAffineTerm(2.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    b1 = MOIU.normalize_and_add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1])
        ], 0.0), MOI.GreaterThan(0.0))

    b2 = MOIU.normalize_and_add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.GreaterThan(0.0))

    MOI.set(bridged,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
        )
    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(bridged)

    obj = MOI.get(bridged, MOI.ObjectiveValue())

    @test obj ≈ -9.33333 atol = 1e-2

    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)

    @test Xr ≈ [1.3333, 1.3333] atol = 1e-2

end

@testset "Simple BLP" begin
    optimizer = Gurobi.Optimizer()
    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOI.empty!(bridged)
    @test MOI.is_empty(bridged)

    X = MOI.add_variables(bridged, 2)

    MOI.set(bridged, MOI.VariableName(), X[1], "x_u")
    MOI.set(bridged, MOI.VariableName(), X[2], "y_u")

    MOI.set(bridged,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
        )
    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    optimizer2 = Gurobi.Optimizer()
    bridged2 = MOIB.full_bridge_optimizer(optimizer2, Float64)
    MOI.empty!(bridged2)
    @test MOI.is_empty(bridged2)

    X2 = MOI.add_variables(bridged2, 2)
    MOI.set(bridged2, MOI.VariableName(), X2[1], "x_l")
    MOI.set(bridged2, MOI.VariableName(), X2[2], "y_l")

    c1 = MOIU.normalize_and_add_constraint(bridged2,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(2.0, X2[1]),
            MOI.ScalarAffineTerm(1.0, X2[2])
        ], 0.0), MOI.LessThan(4.0))
    MOI.set(bridged2, MOI.ConstraintName(), c1, "l_c1")

    c2 = MOIU.normalize_and_add_constraint(bridged2,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X2[1]),
            MOI.ScalarAffineTerm(2.0, X2[2])
        ], 0.0), MOI.LessThan(4.0))
    MOI.set(bridged2, MOI.ConstraintName(), c2, "l_c2")

    b1 = MOIU.normalize_and_add_constraint(bridged2,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X2[1])
        ], 0.0), MOI.GreaterThan(0.0))
    MOI.set(bridged2, MOI.ConstraintName(), b1, "l_b1")

    b2 = MOIU.normalize_and_add_constraint(bridged2,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X2[2])
        ], 0.0), MOI.GreaterThan(0.0))
    MOI.set(bridged2, MOI.ConstraintName(), b2, "l_b12")

    MOI.set(bridged2,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0], [X2[2]]), 0.0)
        )
    MOI.set(bridged2, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    links = Dict(X[1] => X2[1], X[2] => X2[2])

    parametric = [X2[1]]

    blp, _, _ = BilevelJuMP.build_bilivel(bridged, bridged2, links, parametric)

    # MOI.optimize!(bridged)

    # obj = MOI.get(bridged, MOI.ObjectiveValue())

    # @test obj ≈ -9.33333 atol = 1e-2

    # Xr = MOI.get(bridged, MOI.VariablePrimal(), X)

    # @test Xr ≈ [1.3333, 1.3333] atol = 1e-2

    # lp_model = MathOptFormat.MOF.Model()
    # MOI.copy_to(lp_model, blp)
    # MOI.write_to_file(lp_model, "my_model.LP")
    # @show pwd()

    optimizer = Gurobi.Optimizer()
    bridged3 = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOI.empty!(bridged3)
    @test MOI.is_empty(bridged3)
    MOI.copy_to(bridged3, blp)

    MOI.optimize!(bridged3)

end

@testset "Simple BLP2" begin
    optimizer = Gurobi.Optimizer()
    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOI.empty!(bridged)
    @test MOI.is_empty(bridged)

    x = MOI.add_variable(bridged)
    y = MOI.add_variable(bridged)

    MOI.set(bridged, MOI.VariableName(), x, "x_u")
    MOI.set(bridged, MOI.VariableName(), y, "y_u")

    MOI.set(bridged,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 3.0], [x, y]), 0.0)
        )
    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    optimizer2 = Gurobi.Optimizer()
    bridged2 = MOIB.full_bridge_optimizer(optimizer2, Float64)
    MOI.empty!(bridged2)
    @test MOI.is_empty(bridged2)

    x2 = MOI.add_variable(bridged2)
    y2 = MOI.add_variable(bridged2)

    MOI.set(bridged2, MOI.VariableName(), x2, "x_l")
    MOI.set(bridged2, MOI.VariableName(), y2, "y_l")

    c1 = MOIU.normalize_and_add_constraint(bridged2,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, x2),
            MOI.ScalarAffineTerm(1.0, y2)
        ], 0.0), MOI.LessThan(8.0))
    MOI.set(bridged2, MOI.ConstraintName(), c1, "l_c1")

    c2 = MOIU.normalize_and_add_constraint(bridged2,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, x2),
            MOI.ScalarAffineTerm(4.0, y2)
        ], 0.0), MOI.GreaterThan(8.0))
    MOI.set(bridged2, MOI.ConstraintName(), c2, "l_c2")

    c3 = MOIU.normalize_and_add_constraint(bridged2,
    MOI.ScalarAffineFunction([
        MOI.ScalarAffineTerm(1.0, x2),
        MOI.ScalarAffineTerm(2.0, y2)
    ], 0.0), MOI.LessThan(13.0))
    MOI.set(bridged2, MOI.ConstraintName(), c3, "l_c3")

    b1 = MOIU.normalize_and_add_constraint(bridged2,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, x2)
        ], 0.0), MOI.GreaterThan(1.0))
    MOI.set(bridged2, MOI.ConstraintName(), b1, "l_b1")

    b2 = MOIU.normalize_and_add_constraint(bridged2,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, x2)
        ], 0.0), MOI.LessThan(6.0))
    MOI.set(bridged2, MOI.ConstraintName(), b2, "l_b12")

    MOI.set(bridged2,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-1.0], [y2]), 0.0)
        )
    MOI.set(bridged2, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    links = Dict(x => x2, y => y2)

    parametric = [x2]

    blp, _, _ = BilevelJuMP.build_bilivel(bridged, bridged2, links, parametric)

    # MOI.optimize!(bridged)

    # obj = MOI.get(bridged, MOI.ObjectiveValue())

    # @test obj ≈ -9.33333 atol = 1e-2

    # Xr = MOI.get(bridged, MOI.VariablePrimal(), X)

    # @test Xr ≈ [1.3333, 1.3333] atol = 1e-2

    # lp_model = MathOptFormat.MOF.Model()
    # MOI.copy_to(lp_model, blp)
    # MOI.write_to_file(lp_model, "my_model.LP")
    # @show pwd()

    optimizer = Gurobi.Optimizer()
    bridged3 = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOI.empty!(bridged3)
    @test MOI.is_empty(bridged3)
    MOI.copy_to(bridged3, blp)

    MOI.optimize!(bridged3)

end

@testset "Simple BLP JuMP" begin

    model = BilevelModel()

    @variable(UpperToLower(model), x)
    @variable(LowerToUpper(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @objective(Lower(model), Min, y)

    @constraint(Lower(model), 2x+y <= 4)
    @constraint(Lower(model), x+2y <= 4)
    @constraint(Lower(model), x >= 0)
    @constraint(Lower(model), y >= 0)

    optimize!(model, Gurobi.Optimizer())

    primal_status(model)

    termination_status(model)

    @test objective_value(model) == -8

    @test value(x) ==  2
    @test value(y) ==  0

end

@testset "Simple BLP2 JuMP" begin

    model = BilevelModel()

    @variable(UpperToLower(model), x)
    @variable(LowerToUpper(model), y)

    @objective(Upper(model), Min, x + 3y)

    @objective(Lower(model), Min, -y)

    @constraint(Lower(model), x+y <= 8)
    @constraint(Lower(model), x+4y >= 8)
    @constraint(Lower(model), x+2y <= 13)
    @constraint(Lower(model), x >= 1)
    @constraint(Lower(model), x <= 6)
    @constraint(Lower(model), y >= 0)

    optimize!(model, Gurobi.Optimizer())

    primal_status(model)

    termination_status(model)

    @test objective_value(model) == 12

    @test value(x) == 6
    @test value(y) == 2

end

# alternative sintax
# @constraint(model, lower, [y] \in ArgMin([x]))


text:
export add_walsh, add_residue

function add_walsh(mgr::Ptr{Manager}, f_pp::Ptr{Ptr{Node}}, g_pp::Ptr{Ptr{Node}}, n::Int)
    result = ccall((:Cudd_addWalsh, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Ptr{Node}}, Ptr{Ptr{Node}}, Cuint), mgr, f_pp, g_pp, n)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_residue(mgr::Ptr{Manager}, n::Int, m::Int, options::Int, top::Int)
    result = ccall((:Cudd_addResidue, libcudd),
        Ptr{Node}, (Ptr{Manager}, Cuint, Cuint, Cuint, Cuint), mgr, n, m, options, top)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end


text:
#=
Visual illustration of the central limit theorem in 3d

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

References
----------
Based off the original python  file clt3d.py
=#
using Plots
pyplot()
using Distributions
using KernelDensity

beta_dist = Beta(2.0, 2.0)

function gen_x_draws(k)
    bdraws = rand(beta_dist, 3, k)

    # == Transform rows, so each represents a different distribution == #
    bdraws[1, :] -= 0.5
    bdraws[2, :] += 0.6
    bdraws[3, :] -= 1.1

    # == Set X[i] = bdraws[j, i], where j is a random draw from {1, 2, 3} == #
    js = rand(1:3, k)
    X = Array{Float64}(k)
    for i=1:k
        X[i]=  bdraws[js[i], i]
    end

    # == Rescale, so that the random variable is zero mean == #
    m, sigma = mean(X), std(X)
    return (X .- m) ./ sigma
end

nmax = 5
reps = 100000
ns = 1:nmax

# == Form a matrix Z such that each column is reps independent draws of X == #
Z = Array{Float64}(reps, nmax)
for i=ns
    Z[:, i] = gen_x_draws(reps)
end

# == Take cumulative sum across columns
S = cumsum(Z, 2)

# == Multiply j-th column by sqrt j == #
Y = S .* (1. ./ sqrt(ns))'

# == Plot == #
a, b = -3, 3
gs = 100
xs = linspace(a, b, gs)

x_vec = []
y_vec = []
z_vec = []
colors = []
for n=ns
    kde_est = kde(Y[:, n])
    _xs, ys = kde_est.x, kde_est.density
    push!(x_vec, collect(_xs))
    push!(y_vec, ys)
    push!(z_vec, collect(n*ones( length(_xs))))
    push!(colors, RGBA(0, 0, 0, 1-(n-1)/nmax))
end

plot(x_vec, z_vec, y_vec, color=reshape(colors,1,length(colors)), legend=:none)
plot!(xlims=(a,b), xticks=[-3; 0; 3])
plot!(ylims=(1, nmax), yticks=ns, ylabel="n")
plot!(zlims=(0, 0.4), zticks=[0.2; 0.4])


text:
#
# N-body models
#

"""
A `ModelingToolkit.ODESystem` for the Newtonian N-body Problem.

The order of the states follows:
`[x₁, y₁, z₁, ..., xₙ, yₙ, zₙ, ẋ₁, ẏ₁, ż₁, ..., ẋₙ, ẏₙ, żₙ]`.

The order of the parameters follows:
`[G, m₁, m₂, ..., mₙ]`.

!!! warning
    Be careful about specifying `stm=true` for systems with
    `N ≥ 3`! If state transition matrix dynamics are enabled,
    you can calculate the total number of system states with
    `N*6 + (N*6)^2`. Note that this increases exponentially as
    `N` grows! For `N == 6`, unless you're using parallelization,
    your computer may run for several hours.

# Extended Help
The N-body problem is a model which describes how `N` bodies
will move with respect to a common origin. This problem
typically involves many bodies which act due to one force:
electromagentism, gravity, etc. This model applies most
closely to many celestial bodies moving due to gravity.
That's about right for a model in a package called
`AstrodynamicalModels`!

### Usage

```julia
# One body for ALL the planets in our solar system 😎
model = NBP(9)
```
"""
@memoize function NBP(N::Int; stm=false, structural_simplify=true, name=:NBP)

    N > 0 || throw(ArgumentError("`N` must be a number greater than zero!"))
    @parameters t G m[1:N]
    @variables x[1:N](t) y[1:N](t) z[1:N](t) ẋ[1:N](t) ẏ[1:N](t) ż[1:N](t)
    δ = Differential(t)

    r = [[x[i], y[i], z[i]] for i in 1:N]
    v = [[ẋ[i], ẏ[i], ż[i]] for i in 1:N]

    poseqs = reduce(vcat, [
        δ.(r[i]) .~ v[i]
        for i ∈ 1:N
    ])

    veleqs = reduce(vcat, [
        δ.(v[i]) .~ sum((
            ((m[i] * m[j]) / norm(r[j] .- r[i])^3 * (r[j] .- r[i])) .* (j == N ? G / m[i] : 1)
            for j ∈ 1:N if j ≢ i
        ))  for i ∈ 1:N
    ])

    eqs = vcat(poseqs, veleqs)

    if stm
        if N ≥ 3
            @warn """
            You requested state transition matrix dynamics for $N bodies.
            This will result in $(N*6 + (N*6)^2) states! That may take
            a long time! Consider setting `stm=false`, and using
            `ModelingToolkit.calculate_jacobian` instead.
            """
        end

        @variables Φ[1:length(eqs),1:length(eqs)](t)
        Φ = Symbolics.scalarize(Φ)
        A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r...,v...))

        LHS = map(δ, Φ)
        RHS = map(simplify, A * Φ)

        eqs = vcat(eqs, [LHS[i] ~ RHS[i] for i in 1:length(LHS)])
    end

    if string(name) == "NBP" && stm
        modelname = Symbol("NBPWithSTM")
    else
        modelname = name
    end

    sys = ODESystem(
        eqs, t, stm  ? vcat(r...,v...,Φ...) : vcat(r...,v...), vcat(G, m...);
        name = modelname
    )
    return structural_simplify ? ModelingToolkit.structural_simplify(sys) : sys
end

"""
Returns an `ODEFunction` for NBP dynamics.
Results are cached with `Memoize.jl`.
The order of states and parameters in the
`ODEFunction` arguments are equivalent to the
order of states and parameters for the system
produced with `NBP(N)`. As a general rule,
the order of the states follows: `[x₁, y₁, z₁,
..., xₙ, yₙ, zₙ, ẋ₁, ẏ₁, ż₁, ..., ẋₙ, ẏₙ, żₙ]`.

!!! note
    Unlike `R2BP` and `CR3BP`, `jac` is set to
    `false` by default. The number of states for `NBP`
    systems can be very large for relatively small
    numbers of bodies (`N`). Enabling `jac=true`
    by default would cause unnecessarily long
    waiting times for this function to return for
    `N ≥ 3` or so. If `N=2` and `stm=true`,
    setting `jac=true` could still result in several
    minutes of calculations, depending on the computer
    you're using.

!!! warning
    Be careful about specifying `stm=true` for systems with
    `N ≥ 3`! If state transition matrix dynamics are enabled,
    you can calculate the total number of system states with
    `N*6 + (N*6)^2`. Note that this increases exponentially as
    `N` grows! For `N == 6`, unless you're using parallelization,
    your computer may run for several hours.

# Extended Help

### Usage

The `stm`, `structural_simplify`, and `name` keyword arguments
are passed to `NBP`. All other keyword arguments are passed
directly to `SciMLBase.ODEFunction`.

```julia
f = NBPFunction(3; stm=false, structural_simplify=true, name=:NBP, jac=false, sparse=false)
let u = randn(3*6), p = randn(1 + 3), t = 0
    f(u, p, t)
end
```
"""
@memoize function NBPFunction(N::Int; stm=false, structural_simplify=true, name=:R2BP, kwargs...)
    defaults = (; jac=false)
    options  = merge(defaults, kwargs)
    if N ≥ 2 && stm && options.jac
        @warn """
        With state transition matrix dynamics appended to the state
        vector, and Jacobian compuations reqeusted, your system
        will require $((((N*6) + (N*6)^2) + ((N*6) + (N*6)^2)^2)) scalar
        calculations! Consider setting `jac=false`, `stm=false`, or both.
        """
    end
    return ODEFunction(
        NBP(N; stm=stm, structural_simplify=structural_simplify, name=name);
        options...
    )
end


text:
"""
    computeEigenvalues(
        m::AIMSteppingMethod,
        p::QuadraticEigenvalueProblem{N,T},
        c::AIMCache{N,Polynomial{T}};
        plr_polish::Bool = true,
        plr_epsilon::Real = convert(T, 1.0e-10)
        ) where {N <: Unsigned, T <: Number}

Compute the eigenvalues for the problem `p` with corresponding cache `c`.

# Input
- `m::AIMSteppingMethod`: The stepping method to use.
- `p::QuadraticEigenvalueProblem`: The previously defined problem data.
- `c::AIMCache`: The cache constructed from p.
- `plr_polish::Bool`: Tell PolynomialRoots to divide the original polynomial by each root found and polish the results using the full polynomial.
- `plr_epsilon::Real`: The stopping criterion described in Skowron & Gould paper. This is not the precision with which the roots will be calculated.

# Output
An object of type `Array{T,1}` containing the computed eigenvalues.
"""
function computeEigenvalues(
    m::AIMSteppingMethod,
    p::QuadraticEigenvalueProblem{N,T},
    c::AIMCache{N,Polynomial{T}};
    plr_polish::Bool=true,
    plr_epsilon::Real=1.0e-10
    ) where {N <: Unsigned,T <: Number}

    # Compute the AIM "quantization condition"
    δ = computeDelta!(m, p, c)

    # Solve the quantization condition to obtain the eigenvalues
    eigenvalues = PolynomialRoots.roots(coeffs(δ), polish=plr_polish, epsilon=plr_epsilon)

    if isempty(eigenvalues)
        println("Warning: The computed mode array is empty. This means that no roots of the polynomial equation in ω were found.")
    end

    return eigenvalues
end

"""
    computeEigenvalues(
        m::AIMSteppingMethod,
        p::NumericAIMProblem{N,T},
        c::AIMCache{N,T},
        guess::T;
        nls_xtol::Real = convert(T, 1.0e-10),
        nls_ftol::Real = convert(T, 1.0e-10),
        nls_iterations::Int = 1000
        ) where {N <: Unsigned, T <: Complex}

Compute a single eigenvalue for the problem `p` with corresponding cache `c`.

# Input
- `m::AIMSteppingMethod`: The stepping method to use.
- `p::NumericAIMProblem`: The previously defined problem data.
- `c::AIMCache`: The cache constructed from p.
- `nls_xtol::Real`: Norm difference in x between two successive iterates under which convergence is declared.
- `nls_ftol::Real`: Infinite norm of residuals under which convergence is declared.
- `nls_iterations::Int`: Maximum number of iterations performed by NLsolve.

# Output
An object of type `SolverResults` returned by `nlsolve`. See [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) for further details.
"""
function computeEigenvalues(
    m::AIMSteppingMethod,
    p::NumericAIMProblem{N,T},
    c::AIMCache{N,T},
    guess::T;
    nls_xtol::Real=1.0e-10,
    nls_ftol::Real=1.0e-10,
    nls_iterations::Int=1000
    ) where {N <: Unsigned,T <: Complex}

    # This function is passed to NLsolve to find the roots of δ
    function f!(F, x)
        y = computeDelta!(m, p, c, Complex(x[1], x[2]))
        F[1] = real(y)
        F[2] = imag(y)
    end

    # compute the roots of δ using NLsolve
    return nlsolve(f!, [real(guess), imag(guess)], ftol=nls_ftol, xtol=nls_xtol, iterations=nls_iterations)
end

"""
    computeEigenvalues(
        m::AIMSteppingMethod,
        p::NumericAIMProblem{N,T},
        c::AIMCache{N,T},
        guess::T;
        roots_atol::Real = 1.0e-10,
        roots_rtol::Real = 1.0e-10,
        roots_xatol::Real = 1.0e-10,
        roots_xrtol::Real = 1.0e-10,
        roots_maxevals::Int = 100
        ) where {N <: Unsigned, T <: Real}

Compute a single eigenvalue for the problem `p` with corresponding cache `c`.
For details on convergence settings see [Roots.jl](https://juliahub.com/docs/Roots/o0Xsi/1.0.7/reference/#Convergence).

# Input
- `m::AIMSteppingMethod`: The stepping method to use.
- `p::NumericAIMProblem{N,T}`: The previously defined problem data.
- `c::AIMCache{N,T}`: The cache constructed from p.
- `grid::Tuple{T,T}`: A tuple consisting of (start point, end point).
- `roots_atol::Real`: Absolute tolerance.
- `roots_rtol::Real`: Relative tolerance.
- `roots_xatol::Real`: Floating point comparison absolute tolerance.
- `roots_xrtol::Real`: Floating point comparison relative tolerance.
- `roots_maxevals::Int`: Number of algorithm iterations performed.

# Output
An object of type T containing the found eigenvalue.
"""
function computeEigenvalues(
    m::AIMSteppingMethod,
    p::NumericAIMProblem{N,T},
    c::AIMCache{N,T},
    guess::T;
    roots_atol::Real=1.0e-10,
    roots_rtol::Real=1.0e-10,
    roots_xatol::Real=1.0e-10,
    roots_xrtol::Real=1.0e-10,
    roots_maxevals::Int=100
    ) where {N <: Unsigned,T <: Real}

    try
        find_zero(
            x -> computeDelta!(m, p, c, x),
            guess,
            atol=roots_atol,
            rtol=roots_rtol,
            xatol=roots_xatol,
            xrtol=roots_xrtol,
            maxevals=roots_maxevals
            )
    catch
        println("find_zeros was unable to converge to an eigenvalue")
    end
end


text:
using ExtendedKronigPennyMatrix
using Test
using Unitful

@testset "ExtendedKronigPennyMatrix.jl" begin
   include("Alternates.jl")
   include("basic.jl")
   # include("KronigPennyTest.jl")
end


text:
import NamedNamedTuples: @NamedNamedTuple
using Test

import Setfield: @set

@testset "NamedNamedTuples.jl" begin
    # Write your tests here.
    @NamedNamedTuple MyType

    mt = MyType((x=3, y="hey"))
    @test mt.x == 3
    @test mt.y == "hey"
    @test_throws ErrorException mt.z

    mt2 = @set mt.x = 7
    @test mt2.x == 7
    @test mt2.y == "hey"

    # allow field type to change
    mt3 = @set mt.x = "seven"
    @test mt3.x == "seven"
    @test mt3.y == "hey"

    @test_throws MethodError mt.x = 3
    @test_throws ErrorException mt.x = (x=1, y="one")
    @test_throws ErrorException mt.y = "hey"
end


text:
# tests for ApproxManifoldProducts.jl

using ApproxManifoldProducts
using Test

include("basics.jl")

include("ex_1D.jl")

include("ex_2D_rot.jl")

include("testManiProductBigSmall.jl")

include("testBasicManiProduct.jl")

#


text:
# ---
# title: 363. Max Sum of Rectangle No Larger Than K
# id: problem363
# author: Indigo
# date: 2021-05-29
# difficulty: Hard
# categories: Binary Search, Dynamic Programming, Queue
# link: <https://leetcode.com/problems/max-sum-of-rectangle-no-larger-than-k/description/>
# hidden: true
# ---
#
# Given a non-empty 2D matrix _matrix_ and an integer _k_ , find the max sum of
# a rectangle in the _matrix_ such that its sum is no larger than _k_.
#
# **Example:**
#
#
#
#     Input: matrix = [[1,0,1],[0,-2,3]], k = 2
#     Output: 2
#     Explanation:  Because the sum of rectangle [[0, 1], [-2, 3]] is 2,
#                  and 2 is the max number no larger than k (k = 2).
#
# **Note:**
#
#   1. The rectangle inside the matrix must have an area > 0.
#   2. What if the number of rows is much larger than the number of columns?
#
#
## @lc code=start
using LeetCode

function max_sum_submatrix(matrix::AbstractMatrix{Int}, k::Int)
    m, n = size(matrix)
    m > n && return max_sum_submatrix(matrix', k)
    pref_sum = OffsetArray(fill(0, m + 1, n + 1), -1, -1)

    for i in 1:m, j in 1:n
        pref_sum[i, j] = pref_sum[i - 1, j] + pref_sum[i, j - 1] + matrix[i, j] -
                         pref_sum[i - 1, j - 1]
    end
    res = typemin(Int)
    for tp in 1:m, bot in tp:m
        ss = SortedSet{Int}()
        push!(ss, 0)
        for r in 1:n
            rows_sum = pref_sum[bot, r] - pref_sum[tp - 1, r]
            tk = searchsortedfirst(ss, rows_sum - k)
            tk != pastendsemitoken(ss) && (res = max(res, rows_sum - deref((ss, tk))))
            push!(ss, rows_sum)
        end
    end
    return res
end
## @lc code=end