text:
struct BlastResult
    qseqid::String
    sseqid::String
    pident::Float64
    length::Int
    mismatch::Int
    gapopen::Int
    qstart::Int
    qend::Int
    sstart::Int
    send::Int
    evalue::Float64
    bitscore::Float64
end

function BlastResult(line::AbstractString)
    cols = split(line, "\t")

    qseqid = cols[1]
    sseqid = cols[2]
    pident = parse(Float64, cols[3])
    length = parse(Int, cols[4])
    mismatch = parse(Int, cols[5])
    gapopen = parse(Int, cols[6])
    qstart = parse(Int, cols[7])
    qend = parse(Int, cols[8])
    sstart = parse(Int, cols[9])
    send = parse(Int, cols[10])
    evalue = parse(Float64, cols[11])
    bitscore = parse(Float64, cols[12])

    return BlastResult(qseqid, sseqid, pident, length, mismatch, gapopen, qstart, qend, sstart, send, evalue, bitscore)
end

qseqid(record::BlastResult) = record.qseqid
sseqid(record::BlastResult) = record.sseqid

function blastLCA(;filepath::AbstractString, outpath::AbstractString, sqlite::SQLite.DB, taxonomy::Taxonomy.DB, method::Function, header::Bool=false, ranks=[:superkingdom, :phylum, :class, :order, :family, :genus, :species], rmselfhit=true)
    f = open(filepath,"r")
    o = open(outpath,"w")
    try
        lca_ch = blastLCA(f; sqlite=sqlite, taxonomy=taxonomy, method=method, header=header, ranks=ranks, rmselfhit=rmselfhit)
        for (qseqid, taxon, lineage) in lca_ch
            id = taxid(taxon)
            lineage_txt = sprint(io -> print_lineage(io, lineage))
            write(o, "$(qseqid)\t$(id)\t$(lineage_txt)\n")
        end
    finally
        close(f)
        close(o)
    end
end

function blastLCA(f::IOStream; sqlite::SQLite.DB, taxonomy::Taxonomy.DB, method::Function, header::Bool=false, ranks=[:superkingdom, :phylum, :class, :order, :family, :genus, :species], rmselfhit=true)

    blastresult_ch = Channel{BlastResult}(500)
    lcainput_ch = Channel{Tuple{String,Dict{Taxon,BlastResult}}}(500)
    lineage_ch = Channel{Tuple{String,Taxon,Lineage}}(500)

    t1 = @async parse_blastresult!(blastresult_ch, f, header, rmselfhit)
    t2 = @async put_blastresults!(lcainput_ch, blastresult_ch, sqlite, taxonomy)
    t3 = @async lca_blastresults!(lineage_ch, lcainput_ch, method, ranks)

    bind(blastresult_ch, t1)
    bind(lcainput_ch, t2)
    bind(lineage_ch, t3)

    return lineage_ch
end

function parse_blastresult!(out_channel::Channel{BlastResult}, f::IOStream, header::Bool, rmselfhit::Bool)
    header ? readline(f) : nothing
    for line in eachline(f)
        record = BlastResult(line)
        if rmselfhit && qseqid(record) == sseqid(record)
            continue
        end
        put!(out_channel, record)
    end
    close(out_channel)
end

function put_blastresults!(out_channel::Channel{Tuple{String,Dict{Taxon,BlastResult}}}, in_channel::Channel{BlastResult}, sqlite::SQLite.DB, taxonomy::Taxonomy.DB)

    results = Dict{Taxon,BlastResult}()

    while true
        record = take!(in_channel)
        taxid = get(sqlite, sseqid(record), nothing)

        if taxid === nothing
            @warn "record $(sseqid(record)) has no taxid in $(sqlite.file)"
            taxon = nothing
        else
            taxon = get(taxid, taxonomy, nothing)
        end

        if taxon === nothing
            @warn "There is no taxon correspondinig to $(taxid)!\ncontinue..."
        else
            if ! haskey(results, taxon)
                results[taxon] = record
            else
                if record.bitscore > results[taxon].bitscore
                    results[taxon] = record
                end
            end
        end

        try
            next = fetch(in_channel)
            if qseqid(record) != qseqid(next)
                put!(out_channel, (qseqid(record),results))
                results = Dict{Taxon,BlastResult}()
            end
        catch e
            put!(out_channel, (qseqid(record),results))
            break
        end
    end
    close(out_channel)
end

function lca_blastresults!(out_channel::Channel{Tuple{String,Taxon,Lineage}}, in_channel::Channel{Tuple{String,Dict{Taxon,BlastResult}}}, method::Function, ranks)
    for results in in_channel
        records = results[2]
        taxon = method(records)
        lineage = reformat(Lineage(taxon), ranks)
        put!(out_channel, (results[1], taxon, lineage))
    end
    close(out_channel)
end


text:
# Test of DSMC with Harmonic potential
using Random: MersenneTwister
# using BenchmarkTools
using OptEvapCool
using StatsPlots

# Assumes constant trapping frequencies

function harmonic_plots(conditions, sensor, final_cloud, dir, T₀, depth = 10000)
    species = conditions.species

    # THEORY
    ωx, ωy, ωz = 2π * 150,   2π * 150,   2π * 15

    depth = OptEvapCool.time_parametrize(depth)
    duration = last(sensor.time)
    τ_bg = 100#conditions.τ_bg
    N₀ = size(conditions.positions, 2) * conditions.F

    t_series, N_series, T_series, Γ_series = harmonic_theory(
        species, ωx, ωy, ωz, T₀, N₀, depth, duration, τ_bg
    )

    # Plotting
    temperature_plt, T_final = plot_temperature(sensor)
    plot!(t_series, T_series, label = "Theory")

    speed_hist = plot_speed(final_cloud)
    max_speed = maximum(OptEvapCool.speeds(final_cloud))
    plot!(equilibrium_speeds(species.m, T_final, max_speed)...,
        label = "Theory")

    number_plt = plot_number(sensor)
    plot!(t_series, N_series, label = "Theory")
    energy_plt = plot_energy(sensor)
    # TODO - energy theory
    collrate_plt = plot_collrate(sensor)
    plot!(t_series, Γ_series, label = "Theory")

    savefig(temperature_plt, "$dir/temp.png")
    savefig(energy_plt, "$dir/energy.png")
    savefig(speed_hist, "$dir/speed.png")
    savefig(collrate_plt, "$dir/collrate.png")
    savefig(number_plt, "$dir/number.png")
end

function harmonic_none(T₀, Np, duration)
    results = harmonic_test(T₀, Np, duration, name = "harmonic-no-loss")
    harmonic_plots(results..., T₀)
    return nothing
end

function harmonic_background(T₀, Np, duration, τ_bg = 5)
    results = harmonic_test(T₀, Np, duration,
        name = "harmonic-bg-loss", τ_bg = τ_bg)
    harmonic_plots(results..., T₀)
    return nothing
end

function harmonic_spontaneous(T₀, Np, duration, η = 10)
    depth = η * kB * T₀
    results = harmonic_test(T₀, Np, duration,
        name = "harmonic-spontaneous", depth = depth)

    harmonic_plots(results..., T₀, depth)
    return nothing
end

function harmonic_forced(T₀, Np, duration, τ)
    depth = exponential_ramp(start, stop, τ)
    results = harmonic_test(T₀, Np, duration,
        name = "harmonic-forced", depth = depth, τ_bg = 100)
    harmonic_plots(results..., T₀, depth)
    return nothing
end

function harmonic_test(T₀, Np, duration;
        τ_bg = Inf, K = 0, depth = 10000 #= ≈Inf =#, name = "harmonic")

    F = 10; Nc = 5;

    # TEMPORARY
    F = Np / 10000

    ωx, ωy, ωz = 2π * 150,   2π * 150,   2π * 15
    field = HarmonicField(ωx, ωy, ωz)

    Nt = ceil(Int64, Np / F)

    species = Rb87
    m = Rb87.m

    positions = harmonic_boltzmann_positions(Nt, m, T₀, ωx, ωy, ωz)
    velocities = boltzmann_velocities(Nt, m, T₀)

    # Function to make measurements on the system
    sensor = GlobalSensor()
    measure = measurer(sensor)

    poten = potential(field)
    evap = energy_evap(depth, poten)

    conditions = SimulationConditions(species, F, positions, velocities,
        acceleration(field), poten, K = K, evap = evap, τ_bg = τ_bg)

    max_dt = 0.05 * 2π / max(ωx, ωy, ωz)
    # Run evolution
    final_cloud = evolve(conditions, duration;
        Nc = Nc, max_dt = max_dt, measure = measure)

    # Save files
    ft = filetime()
    dir = "./results/$ft-$name"
    mkpath(dir)

    # Save data to CSV file
    savecsv(sensor, "$dir/sensor-data.csv")

    # Save simulation parameters
    parameter_str = """
    # Harmonic potential test
    Name: $name

    ## Conditions
    Initial temperature: $T₀ K
    Initial number of real particles: $Np
    Duration: $duration
    Background loss time constant: $τ_bg
    Three-body loss rate constant: $K
    Evaporation?: $(evap != no_evap)
    Trapping frequencies:
        $ωx
        $ωy
        $ωz

    ## Numerical parameters
    F = $F
    Nc = $Nc
    Max. timestep: $max_dt
    """

    io = open("$dir/parameters.txt", "w");
    write(io, parameter_str)
    close(io)

    return conditions, sensor, final_cloud, dir
end

ext: jl


text:
#!/usr/bin/env julia

using Luxor, Colors

function get_path(str)
  p = textpath(str)
  o = getpathflat()
  stroke()
  sethue("red")
  x, y = 0, 0
  for e in o
        if e.element_type == Cairo.CAIRO_PATH_MOVE_TO
            (x, y) = e.points
            move(x, y)
        elseif e.element_type == Cairo.CAIRO_PATH_LINE_TO
            (x, y) = e.points
            # straight lines
            line(x, y)
            stroke()
            circle(x, y, 1, :stroke)
        elseif e.element_type == Cairo.CAIRO_PATH_CURVE_TO
            (x1, y1, x2, y2, x3, y3) = e.points
            #bezier control lines
            circle(x1, y1, 1, :stroke)
            circle(x2, y2, 1, :stroke)
            circle(x3, y3, 1, :stroke)
            move(x, y)
            curve(x1, y1, x2, y2, x3, y3)
            stroke()
            (x, y) = (x3, y3) # update current point
        elseif e.element_type == Cairo.CAIRO_PATH_CLOSE_PATH
            closepath()
        else
            error("unknown CairoPathEntry " * repr(e.element_type))
            error("unknown CairoPathEntry " * repr(e.points))
        end
    end
end

function get_path_test(fname)
    Drawing(500, 500, fname)
    origin()
    sethue("blue")
    fontface("Aachen")
    fontsize(80)
    setline(0.5)
    setopacity(.7)
    get_path("Luxor")
    finish()
    println("finished test: output in $(fname)")
end

get_path_test("/tmp/get-path-test.pdf")


text:
function JetProp2DAcoVTIDenQ_DEO2_FDTD(;
        v = Float32[],
        b = Float32[],
        ϵ = Float32[],
        η = Float32[],
        f = 0.85,
        srcfieldfile = joinpath(tempdir(), "field-$(uuid4()).bin"),
        comptype = nothing,
        compscale = 1e-2,
        isinterior = false,
        nz_subcube = 32,
        nx_subcube = 32,
        sz = 0.0,
        sx = 0.0,
        st = 0.0,
        interpmethod = :hicks,
        rz = [0.0],
        rx = [0.0],
        z0 = 0.0,
        x0 = 0.0,
        dz = 10.0,
        dx = 10.0,
        freqQ = 5.0,
        qMin = 0.1,
        qInterior = 100.0,
        padz = 0.0,
        padx = 0.0,
        dtmod,
        dtrec,
        ntrec,
        nbz_cache = 512,
        nbx_cache = 8,
        nbz_inject = 16,
        nbx_inject = 16,
        nsponge = 50,
        wavelet = WaveletCausalRicker(f=5.0),
        freesurface = false,
        imgcondition = "standard",
        nthreads = Sys.CPU_THREADS,
        reportinterval = 500)
    # active and passive earth model properties.  The active set is in the model-space
    active_modelset = Dict{String,Int}()
    passive_modelset = Dict{String,Array{Float32,2}}()
    i = 1
    for (n,x) in (("v",v), ("b",b), ("ϵ",ϵ), ("η",η))
        if length(x) > 0
            passive_modelset[n] = x
        else
            active_modelset[n] = i # so that m[:,:,i] corresponds to property n∈("v","b","ϵ","η")
            i += 1
        end
    end
    @assert length(active_modelset) > 0

    # active and passive wavefields (an active wavefield is serialized to disk)
    active_modelset_keys = keys(active_modelset)
    local active_wavefields, modeltype
    if "v" ∈ active_modelset_keys && "b" ∉ active_modelset_keys && "ϵ" ∈ active_modelset_keys && "η" ∈ active_modelset_keys
        active_wavefields = ["pold","mold","pspace","mspace"]
        modeltype = WaveFD.Prop2DAcoVTIDenQ_DEO2_FDTD_Model_VEA()
    elseif "v" ∈ active_modelset_keys && "b" ∉ active_modelset_keys && "ϵ" ∉ active_modelset_keys && "η" ∉ active_modelset_keys
        active_wavefields = ["pspace","mspace"]
        modeltype = WaveFD.Prop2DAcoVTIDenQ_DEO2_FDTD_Model_V()
    else
        error("unsupported model-space")
    end

    # dtrec must be integer multiple of dtmod
    @assert abs(dtrec - dtmod*round(dtrec/dtmod,RoundNearest)) < eps(Float32)

    # source location (sz,sx,st)
    sz,sx,st = map(val->Float64[val...], (sz,sx,st))

    # domain and range
    dom = JetSpace(Float32, size(first(values(passive_modelset)))..., length(active_modelset))
    rng = JetSpace(Float32, ntrec, length(rz))

    # type conversions
    z0, x0, padz, padx = map(val->Float64(val), (z0, x0, padz, padx)) # used for computing integer grid locations for source injection
    rz, rx = map(val->convert(Array{Float64}, val), (rz, rx))
    dz, dx, dtmod, dtrec, f, freqQ, qMin, qInterior = map(val->Float32(val), (dz, dx, dtmod, dtrec, f, freqQ, qMin, qInterior)) # used in F.D. calculations
    @assert length(rz) == length(rx)
    if isa(wavelet, Array) == true
        wavelet = convert(Array{Float32}, wavelet)
    end

    # ginsu view of earth model... ginsu is aware of the sponge
    nsponge_top = freesurface ? 0 : nsponge
    padz_top = freesurface ? 0.0 : padz

    ginsu = Ginsu((z0,x0), (dz,dx), size(b), (sz,sx), (rz,rx), ((padz_top,padz),(padx,padx)), ((nsponge_top,nsponge),(nsponge,nsponge)), T=Float32)

    # if srcfieldfile is specified make sure its containing folder exists
    if srcfieldfile != ""
        srcfieldpath = split(srcfieldfile, "/")[1:end-1]
        if length(srcfieldpath) > 0
            srcfieldpath = join(srcfieldpath, "/")
            mkpath(srcfieldpath)
        end
    end

    # compression for nonlinear source wavefields
    C = WaveFD.comptype(comptype, Float32)[1]
    compressor = Dict{String,WaveFD.Compressor{Float32,Float32,C,2}}()

    # we need to serialize for the data, but not for the linearization
    _active_wavefields = "pold" ∈ active_wavefields ? active_wavefields : [active_wavefields;"pold"]

    for active_wavefield in _active_wavefields
        compressor[active_wavefield] = WaveFD.Compressor(Float32, Float32, C, size(ginsu,interior=isinterior),
            (nz_subcube,nx_subcube), compscale, ntrec, isinterior)
    end

    # imaging condition
    icdict = Dict(
        lowercase("standard") => WaveFD.ImagingConditionStandard(),
        lowercase("FWI") => WaveFD.ImagingConditionWaveFieldSeparationFWI(),
        lowercase("RTM") => WaveFD.ImagingConditionWaveFieldSeparationRTM())

    if lowercase(imgcondition) ∉ keys(icdict)
        error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
    end

    Jet(
        dom = dom,
        rng = rng,
        f! = JopProp2DAcoVTIDenQ_DEO2_FDTD_f!,
        df! = JopProp2DAcoVTIDenQ_DEO2_FDTD_df!,
        df′! = JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!,
        s = (
            modeltype = modeltype,
            active_modelset = active_modelset,
            passive_modelset = passive_modelset,
            active_wavefields = active_wavefields,
            f = Float32(f),
            srcfieldfile = srcfieldfile,
            srcfieldhost = Ref(""),
            chksum = Ref(zero(UInt32)),
            compressor = compressor,
            isinterior = isinterior,
            sz = sz,
            sx = sx,
            st = st,
            interpmethod = interpmethod,
            rz = rz,
            rx = rx,
            z0 = z0,
            x0 = x0,
            dz = dz,
            dx = dx,
            freqQ = freqQ,
            qMin = qMin,
            qInterior = qInterior,
            ginsu = ginsu,
            nsponge = nsponge,
            dtmod = dtmod,
            dtrec = dtrec,
            ntrec = ntrec,
            nbz_cache = nbz_cache,
            nbx_cache = nbx_cache,
            nbz_inject = nbz_inject,
            nbx_inject = nbx_inject,
            wavelet = wavelet,
            freesurface = freesurface,
            imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()),
            nthreads = nthreads,
            reportinterval = reportinterval,
            stats = Dict{String,Float64}("MCells/s"=>0.0, "%io"=>0.0, "%inject/extract"=>0.0, "%imaging"=>0.0)))
end

@doc """
    JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; kwargs...)
    JopLnProp2DAcoVTIDenQ_DEO2_FDTD(; m₀, kwargs...)

Create a `Jets` nonlinear or linearized operator for 2D pseudo- visco-acoustic, vertical transverse isotropy,
variable density modeling.

# Model Parameters
This propagator operates with four model parameters, as shown in the table below.

| Parameter | Description                                       |
|:---------:|:--------------------------------------------------|
| `v`       | P wave velocity                                   |
| `ϵ`       | Modified Thomsen's weak anisotropy parameter      |
| `η`       | Modified Alkhalifah's weak anisotropy parameter η |
| `b`       | buoyancy (reciprocal density)                     |

## Pseudo-acoustic approximation
With the pseudo-acoustic approximation we employ, the transformation to self-adjoint form requires a
parameter `f` representing the average ratio of shear-wave to p-wave velocity, and a modfication
to Alkhalifah's weak anisotropy parameter `η`. We show formulas for these two parameters below.
The parameter `f` is specified as a scalar for the entire model, with default value 0.85, implying
a shear velocity around 38% of the p wave velocity.
* `η = sqrt[2(ϵ-δ)/(f+2ϵ)]`
* `f = 1 - Vₛ² / Vₚ²`

For more information about the pseudo-acoustic approximation we employ please see
https://library.seg.org/doi/10.1190/segam2016-13878451.1.

## Active and passive model parameters
There are two modes of operation that define different sets of **active** and **passive** parameters,
*velocity only* and *velocity and anisotropy*.

An **active** parameter can be inverted for using the Jacobian linearized operator machinery,
and a **passive** parameter is constant. The two modes are:

| Mode                    | Active Parameters | Passive Parameters |
|:------------------------|:------------------|:-------------------|
| velocity only           | [`v`]             | [`ϵ`, `η`, `b`]    |
| velocity and anisotropy | [`v`, `ϵ`, `η`]   | [`b`]              |

To make a parameter **passive**, you pass the array for that parameter in the constructor for the
operator, implying that it is in state and does not change with the action of the operator.
Parameters that are *not* passed to the constructor are assumed to be **active**, and will be part
of the model that the operator acts on, stored as a 3D array with the following indices:
`[Z coord][X coord][Active Parameter]`

# Examples

## Model and acquisition geometry setup
1. load modules Jets, WaveFD, and JetPackWaveFD
1. set up the model discretization, coordinate size and spacing
1. set up the acquisition geometry, including the time discretization, locations for source and receivers, and the source wavelet
1. create constant models for `v`, `ϵ`, `η`, `b`
```
using Jets, WaveFD, JetPackWaveFD
nz,nx = 200, 100                      # spatial discretization size
dz,dx = 20.0, 20.0                    # spatial discretization sampling
ntrec = 1101                          # number of temporal samples in recorded data
dtrec = 0.0100                        # sample rate for recorded data
dtmod = 0.0025                        # sample rate for modeled data
wavelet = WaveletCausalRicker(f=5.0)  # type of wavelet to use (source signature)
sz = dz                               # source Z location
sx = dx*(nx/2)                        # source X location
rx = dx*[0:0.5:nx-1;];                # array of receiver X locations
rz = 2*dz*ones(length(0:0.5:nx-1));   # array of receiver Z locations
b = ones(Float32, nz, nx);            # constant buoyancy
ϵ = 0.1 * ones(Float32, nz, nx);      # constant epsilon
η = 0.2 * ones(Float32, nz, nx);      # constant eta
```

## Construct and apply the nonlinear operator (`v` is active; `ϵ, η, b` are passive)
1. create the nonlinear operator `F`
1. create the model vector and set parameter `v`
1. perform nonlinear forward modeling with model `m₀` and return the resulting modeled data in `d`
```
F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ϵ = ϵ, η = η, ntrec = ntrec, dtrec = dtrec,
    dtmod = dtmod, dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
d = F*m₀;              # forward nonlinear op
```

## Construct and apply the nonlinear operator (`v, ϵ, η` are active; `b` is passive)
1. create the nonlinear operator `F`
1. create the model vector and set parameters `v, ϵ, η`
1. perform nonlinear forward modeling with model `m₀` and return the resulting modeled data in `d`
```
F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, dtrec = dtrec, dtmod = dtmod,
    dz = dz, dx = dx, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
d = F*m₀;              # forward nonlinear op
```

## Construct and apply the linearized Jacobian operator (method 1)
For this example we assume in the model that `v, ϵ, η` are active parameters, and `b` is passive.
1. create the nonlinear operator `F` directly by construction of a jet.
1. create the model vector and set parameters `v, ϵ, η`
1. create the Jacobian operator `J` by linearization of `F` at point `m₀`
1. create a random model perturbation vector `δm`.
1. perform linearized forward (Born) modeling on the model perturbation vector `δm` and
    return the resulting data perturbation in `δd`.
1. perform linearized adjoint (Born) migration on the data perturbation vector `δd` and
    return the resulting model perturbation in `δm`.

Note that the Jacobian operators `J` and `J'` require the serialized nonlinear forward wavefield,
and this is generated automatically whenever required. If you watch the logging to standard out
from this example, you will first see the finite difference evolution for the nonlinear forward,
followed by the linearized forward, and finally the linearized adjoint.
```
F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(; b = b, ntrec = ntrec, dz = dz, dx = dx,
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
m₀ = zeros(domain(F));
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = jacobian(F, m₀)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op
```

## Construct and apply the linearized Jacobian operator (method 2)
For this example we assume in the model that `v, ϵ, η` are active parameters, and `b` is passive.
1. create the model vector and set parameters `v, ϵ, η`
1. create the Jacobian operator `J` at point `m₀` directly by construction of a jet.
1. create a random model perturbation vector `δm`.
1. perform linearized forward (Born) modeling on the model perturbation vector `δm` and
    return the resulting data perturbation in `δd`.
1. perform linearized adjoint (Born) migration on the data perturbation vector `δd` and
    return the resulting model perturbation in `δm`.
```
m₀ = zeros(Float32, nz, nx, 3);
m₀[:,:,1] .= 1500;
m₀[:,:,2] .= 0.1;
m₀[:,:,3] .= 0.2;
J = JopLnProp2DAcoVTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, ntrec = ntrec, dz = dz, dx = dx,
    dtrec = dtrec, dtmod = dtmod, wavelet = wavelet, sz = sz, sx = sx, rz = rz, rx = rx)
δm = rand(domain(J));
δd = J*δm;             # forward linearized op
δm = J'*δd;            # adjoint linearized op
```

# Required Parameters
* `m₀` the point at which the jet is linearized. Note this argument is required in the constuctor for
    `JopLnProp2DAcoVTIDenQ_DEO2_FDTD` but not `JopNlProp2DAcoVTIDenQ_DEO2_FDTD`. This constuctor is shown
    in the last example above. Please note that you must consider which parameters are active and passive,
    per the discussion on model parameters above and examples.
* `dtmod` The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate
    with the following expression: `dt = 0.75 * 0.38 * max(dx,dz) / maximum(v)`. Note that the usual
    Courant–Friedrichs–Lewy condition (CFL condition) for stability in finite difference modeling is
    modified *heuristically* to include the impact of the visco-acoustic implementation of this operator,
    requiring a 25% smaller smaple rate.
* `dtrec` The number of time samples in the recorded data. Note requirement that `dtrec`
    be an even multiple of `dtmod`.
* `ntrec` The number of time samples in the recorded data. Note that the number of samples in the modeled
    data is determined from `ntrec`, `dtrec`, and `dtmod`.

# Optional Parameters
Defaults for arguments are shown inside square brackets.

* `b [ones(Float32,0,0)]` the buoyancy (reciprocal density) array.
* `f [0.85]` See the discussion above regarding the **Pseudo-acoustic approximation** used here
* `srcfieldfile [joinpath(tempdir(), "field-$(uuid4()).bin")]` the full path to a scratch file used for
    the serializationof the compressed nonlinear source wavefield.
* `comptype [nothing]` the type of compression to use for the serialization of
    the nonlinear source wavefield. The type of compression must be one of:
    * `nothing` - no compression.
    * `Float32` - if wavefield is `Float64`, then do a simple conversion to Float32.
    * `UInt32` - compression using CvxCompress (windowing + 2D wavelet transform +
        thresholding + quantization + run-length-encoding).
* `compscale [1e-2]` determines the thresholding for the compression of the nonlinear source
    wavefield prior to serialization. Larger values mean more aggressive compression. You can
    likely increase compscale to 1.0 before you start to notice differences in output from
    Jacobian operations.
* `nz_subcube, nx_subcube [32]` The Z and X sizes of windows used for compression
    of the nonlinear source wavefield with `CvxCompress`. Note the requirement
    `[8 <= n*_subcube <=256]`.
* `isinterior [false]` boolean flag that indicates how the nonlinear source wavefield is
    saved. For large models, operation will be faster with `isinterior = true`, but
    the linearization correctness test may fail.
    * `true` the entire model including absorbing boundaries is serialized and deserialized
    * `false` the interior part of the model excluding absorbing boundaries is serialized
        and deserialized
* `sz [0.0]` Array of source Z coordinate. Note that if multiple sources are provided,
    they will will be injected simultaneously during finite difference evolution.
* `sx [0.0]` Array of source X coordinate.
* `st [0.0]` Array of source delay times.
* `interpmethod [:hicks]` Type of physical interpolation for sources and receivers. For locations
    that are not on the physical grid coordinates, interpolation is used either to inject or extract
    information. `interpmethod` must be one of:
    * `:hicks` Hicks 2D sinc interpolation (up to 8x8 nonzero points per location)
    * `:linear` bilinear interpolation (up to 2x2 nonzero points per location)
* `rz [[0.0]]` 2D array of receiver Z coordinates
* `rx [[0.0]]` 2D array of receiver Z coordinates
* `z0 [0.0]` Origin of physical coordinates for the Z dimension.
* `x0 [0.0]` Origin of physical coordinates for the X dimension.
* `dz [10.0]` Spacing of physical coordinates in the Z dimension.
* `dx [10.0]` Spacing of physical coordinates in the X dimension.
* `freqQ [5.0]` The center frequency for the Maxwell body approximation to dissipation only attenuation.
    Please see `JetPackWaveFD` package documentation for more information concerning the attenuation model.
* `qMin [0.1]` The minimum value for Qp at the boundary of the model used in our Maxwell body approximation
    to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the
    attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the
    boundaries of the computational domain.
    Please see `JetPackWaveFD` package documentation for more information concerning the attenuation model.
* `qInterior [100.0]` the value for Qp in the interior of the model used in our Maxwell body approximation
    to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is
    a physically meaningful value.
    Please see `JetPackWaveFD` package documentation for more information concerning the attenuation model.
* `padz, padx [0.0], [0.0]` - apply extra padding to the survey determined aperture in `Ginsu`.
    Please see `Ginsu` for more information.
* `nbz_cache, nbx_cache [512], [8]` The size of cache blocks in the Z and X dimensions.
    In general the cache block in the Z (fast) dimension should be ≥ the entire size of that dimension,
    and the cache block size in the slower dimension is generally small in order to allow the entire
    block to fit in cache.
* `nbz_inject, nbx_inject = [16], [16]` The number of blocks in the Z and X dimensions for threading
    the wavefield injection.
* `nsponge [50]` The number of grid cells to use for the absorbing boundary. For high fidelity modeling
    this should be > 60 grid points, but can be significantly smaller for some use cases like low frequency
    full waveform inversion.
* `wavelet [WaveletCausalRicker(f=5.0)]` The source wavelet, can be specified as either a Wavelet type
    or an array.
* `freesurface [false]` Determines if a free surface (`true`) or absorbing (`false`) top boundary condition
    is applied.
* `imgcondition` ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI",
    and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition
    in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for
    RTM). Note the true adjoint only exists for "standard" imaging condition currently.
* `nthreads [Sys.CPU_THREADS]` The number of threads to use for OpenMP parallelization of the modeling.
* `reportinterval [500]` The interval at which information about the propagtion is logged.

See also: `Ginsu`, `WaveletSine`, `WaveletRicker`, `WaveletMinPhaseRicker`, `WaveletDerivRicker`,
    `WaveletCausalRicker`, `WaveletOrmsby`, `WaveletMinPhaseOrmsby`
"""
JopNlProp2DAcoVTIDenQ_DEO2_FDTD(;kwargs...) = JopNl(JetProp2DAcoVTIDenQ_DEO2_FDTD(;kwargs...))

@doc (@doc JopNlProp2DAcoVTIDenQ_DEO2_FDTD)
JopLnProp2DAcoVTIDenQ_DEO2_FDTD(; m₀, kwargs...) = JopLn(JetProp2DAcoVTIDenQ_DEO2_FDTD(;kwargs...), m₀)

export JopNlProp2DAcoVTIDenQ_DEO2_FDTD
export JopLnProp2DAcoVTIDenQ_DEO2_FDTD

function JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::AbstractArray; kwargs...)
    # make propagator
    nz_ginsu,nx_ginsu = size(kwargs[:ginsu])
    z0_ginsu,x0_ginsu = origin(kwargs[:ginsu])
    p = Prop2DAcoVTIDenQ_DEO2_FDTD(
        nz = nz_ginsu,
        nx = nx_ginsu,
        nbz = kwargs[:nbz_cache],
        nbx = kwargs[:nbx_cache],
        dz = kwargs[:dz],
        dx = kwargs[:dx],
        dt = kwargs[:dtmod],
        nsponge = kwargs[:nsponge],
        nthreads = kwargs[:nthreads],
        freesurface = kwargs[:freesurface],
        freqQ = kwargs[:freqQ],
        qMin = kwargs[:qMin],
        qInterior = kwargs[:qInterior])

    # wave-fields
    wavefields = Dict("pcur"=>WaveFD.PCur(p), "pold"=>WaveFD.POld(p), "pspace"=>WaveFD.PSpace(p), "mcur"=>WaveFD.MCur(p), "mold"=>WaveFD.MOld(p), "mspace"=>WaveFD.MSpace(p))

    # earth model
    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))

    # ginsu'd earth model (active-set)
    for prop in keys(kwargs[:active_modelset])
        sub!(model_ginsu[prop], kwargs[:ginsu], @view(m[:,:,kwargs[:active_modelset][prop]]), extend=true)
    end
    for prop in keys(kwargs[:passive_modelset])
        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
    end
    model_ginsu["f"] .= kwargs[:f]

    ginsu_interior_range = interior(kwargs[:ginsu])

    # we need to serialize "pold" for the data but, not (necessarily) for the linearization
    active_wavefields = "pold" ∈ kwargs[:active_wavefields] ? kwargs[:active_wavefields] : [kwargs[:active_wavefields]; "pold"]

    it0, ntmod_wav = WaveFD.default_ntmod(kwargs[:dtrec], kwargs[:dtmod], kwargs[:st], kwargs[:ntrec])

    # source wavelet for injection, one for each source location
    wavelet_realization = realizewavelet(kwargs[:wavelet], kwargs[:sz], kwargs[:sx], kwargs[:st], kwargs[:dtmod], ntmod_wav)

    # Get source and receiver interpolation coefficients
    local iz_sou, ix_sou, c_sou
    if kwargs[:interpmethod] == :hicks
        iz_sou, ix_sou, c_sou = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:sz], kwargs[:sx])
    else
        iz_sou, ix_sou, c_sou = WaveFD.linearcoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:sz], kwargs[:sx])
    end
    blks_sou = WaveFD.source_blocking(nz_ginsu, nx_ginsu, kwargs[:nbz_inject], kwargs[:nbx_inject], iz_sou, ix_sou, c_sou)

    c_sou_scaled = Array{Array{Float32,2},1}(undef, length(c_sou))
    for i = 1:length(c_sou_scaled)
        c_sou_scaled[i] = similar(c_sou[i])
        for ix = 1:size(c_sou_scaled[i], 2), iz = 1:size(c_sou_scaled[i], 1)
            jz = iz_sou[i][iz]
            jx = ix_sou[i][ix]
            c_sou_scaled[i][iz,ix] = c_sou[i][iz,ix] * kwargs[:dtmod]^2 * model_ginsu["v"][jz,jx]^2 / model_ginsu["b"][jz,jx]
        end
    end
    blks_sou_scaled = WaveFD.source_blocking(nz_ginsu, nx_ginsu, kwargs[:nbz_inject], kwargs[:nbx_inject], iz_sou, ix_sou, c_sou_scaled)

    local iz_rec, ix_rec, c_rec
    if length(d) > 0
        if kwargs[:interpmethod] == :hicks
            iz_rec, ix_rec, c_rec = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
        else
            iz_rec, ix_rec, c_rec = WaveFD.linearcoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
        end
    end

    # disk-file for source-field serialization
    iofield = Dict{String,IOStream}()
    if kwargs[:srcfieldfile] != ""
        for active_wavefield in active_wavefields
            filename = "$(kwargs[:srcfieldfile])-$(active_wavefield)"
            try
                if isfile(filename) == true
                    rm(filename)
                end
                iofield[active_wavefield] = open(filename, "w")
            catch
                @info "Unable to open $(filename) on proc $(myid()) - $(gethostname())"
                rethrow()
            end
            open(kwargs[:compressor][active_wavefield])
        end
    end

    itskip = round(Int, kwargs[:dtrec]/kwargs[:dtmod])
    time1 = time()
    cumtime_io, cumtime_ex, cumtime_pr = 0.0, 0.0, 0.0
    kwargs[:reportinterval] == 0 || @info "nonlinear forward on $(gethostname()), srcfieldfile=$(kwargs[:srcfieldfile])"

    set_zero_subnormals(true)
    for it = 1:ntmod_wav
        if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod_wav)
            JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1, cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
        end

        # propagate and wavefield swap
        cumtime_pr += @elapsed WaveFD.propagateforward!(p)
        wavefields["pcur"],wavefields["pold"] = wavefields["pold"],wavefields["pcur"]
        wavefields["mcur"],wavefields["mold"] = wavefields["mold"],wavefields["mcur"]

        # inject source wavelet
        cumtime_ex += @elapsed begin
            WaveFD.injectdata!(wavefields["pcur"], blks_sou_scaled, wavelet_realization, it)
            WaveFD.injectdata!(wavefields["pspace"], blks_sou, wavelet_realization, it)
        end

        if it >= it0 && rem(it-it0,itskip) == 0
            if length(d) > 0
                # extract receiver data
                cumtime_ex += @elapsed WaveFD.extractdata!(d, wavefields["pold"], div(it-it0,itskip)+1, iz_rec, ix_rec, c_rec)
            end

            # scale spatial derivatives by v^2/b to make them temporal derivatives
            cumtime_io += @elapsed WaveFD.scale_spatial_derivatives!(p)

            if kwargs[:srcfieldfile] != ""
                cumtime_io += @elapsed for active_wavefield in active_wavefields
                    if kwargs[:isinterior]
                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], div(it-it0,itskip)+1,
                            wavefields[active_wavefield], ginsu_interior_range)
                    else
                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], div(it-it0,itskip)+1,
                            wavefields[active_wavefield])
                    end
                end
            end
        end
    end
    set_zero_subnormals(false)
    JopProp2DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod_wav, time()-time1, cumtime_io, cumtime_ex)

    if kwargs[:srcfieldfile] != ""
        for active_wavefield in active_wavefields
            close(iofield[active_wavefield])
            close(kwargs[:compressor][active_wavefield])
        end
    end

    free(p)

    nothing
end

function JopProp2DAcoVTIDenQ_DEO2_FDTD_f!(d::AbstractArray, m::AbstractArray{Float32}; kwargs...)
    d .= 0
    isvalid, _chksum = isvalid_srcfieldfile(m, kwargs[:srcfieldhost][], kwargs[:srcfieldfile]*"-pold", kwargs[:chksum][])
    if !isvalid
        JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d, m; kwargs...)
        kwargs[:chksum][] = _chksum
        kwargs[:srcfieldhost][] = gethostname()
    else
        field = Array{Float32}(undef,size(kwargs[:ginsu], interior=kwargs[:isinterior]))
        local iz, ix, c
        if kwargs[:interpmethod] == :hicks
            iz, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx],
                origin(kwargs[:ginsu], interior=kwargs[:isinterior])...,
                size(kwargs[:ginsu], interior=kwargs[:isinterior])..., kwargs[:rz], kwargs[:rx])
        else
            iz, ix, c = WaveFD.linearcoeffs(kwargs[:dz], kwargs[:dx],
                origin(kwargs[:ginsu], interior=kwargs[:isinterior])...,
                size(kwargs[:ginsu], interior=kwargs[:isinterior])..., kwargs[:rz], kwargs[:rx])
        end

        iofield = open("$(kwargs[:srcfieldfile])-pold")
        open(kwargs[:compressor]["pold"])
        for it = 1:kwargs[:ntrec]
            WaveFD.compressedread!(iofield, kwargs[:compressor]["pold"], it, field)
            WaveFD.extractdata!(d, field, it, iz, ix, c)
        end
        close(iofield)
        close(kwargs[:compressor]["pold"])
    end
    d
end

function JopProp2DAcoVTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArray; kwargs...)
    nz_ginsu,nx_ginsu = size(kwargs[:ginsu])
    z0_ginsu,x0_ginsu = origin(kwargs[:ginsu])
    p = Prop2DAcoVTIDenQ_DEO2_FDTD(
        nz = nz_ginsu,
        nx = nx_ginsu,
        nbz = kwargs[:nbz_cache],
        nbx = kwargs[:nbx_cache],
        dz = kwargs[:dz],
        dx = kwargs[:dx],
        dt = kwargs[:dtmod],
        nsponge = kwargs[:nsponge],
        nthreads = kwargs[:nthreads],
        freesurface = kwargs[:freesurface],
        freqQ = kwargs[:freqQ],
        qMin = kwargs[:qMin],
        qInterior = kwargs[:qInterior])

    # wave-fields
    pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)

    # earth model
    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))

    # ginsu'd earth model (active-set)
    for prop in keys(kwargs[:active_modelset])
        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
    end
    for prop in keys(kwargs[:passive_modelset])
        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
    end
    earth_ginsu["f"] .= kwargs[:f]

    δm_ginsu = Dict{String,Array{Float32,2}}()
    for prop in keys(kwargs[:active_modelset])
        δm_ginsu[prop] = sub(kwargs[:ginsu], @view(δm[:,:,kwargs[:active_modelset][prop]]), extend=false)
    end

    ginsu_interior_range = interior(kwargs[:ginsu])

    # pre-compute receiver interpolation coefficients
    local iz,ix,c
    if kwargs[:interpmethod] == :hicks
        iz, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
    else
        iz, ix, c = WaveFD.linearcoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
    end

    # if necessary, re-run the nonlinear forward
    isvalid, _chksum = isvalid_srcfieldfile(kwargs[:mₒ], kwargs[:srcfieldhost][], kwargs[:srcfieldfile]*"-pold", kwargs[:chksum][])
    if !isvalid
        JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(Array{Float32}(undef,0,0), kwargs[:mₒ]; kwargs...)
        kwargs[:srcfieldhost][] = gethostname()
        kwargs[:chksum][] = _chksum
    end

    # local disk file for source wavefield deserialization
    wavefields = Dict{String,Array{Float32,2}}()
    iofields = Dict{String,IOStream}()
    for active_wavefield in kwargs[:active_wavefields]
        iofields[active_wavefield] = open(kwargs[:srcfieldfile]*"-"*active_wavefield)
        open(kwargs[:compressor][active_wavefield])
        wavefields[active_wavefield] = zeros(Float32, nz_ginsu, nx_ginsu)
    end

    ntmod = WaveFD.default_ntmod(kwargs[:dtrec], kwargs[:dtmod], kwargs[:ntrec])

    δdinterp = zeros(Float32, ntmod, size(δd, 2))

    itskip = round(Int, kwargs[:dtrec]/kwargs[:dtmod])
    kwargs[:reportinterval] == 0 || @info "linear forward on $(gethostname()), srcfieldfile=$(kwargs[:srcfieldfile])"
    time1 = time()
    cumtime_io, cumtime_ex, cumtime_im, cumtime_pr = 0.0, 0.0, 0.0, 0.0

    set_zero_subnormals(true)
    for it = 1:ntmod
        if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod)
            JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1,
                cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "forward")
        end

        # propagate and swap wavefields
        cumtime_pr += @elapsed WaveFD.propagateforward!(p)
        pcur,pold = pold,pcur

        if rem(it-1,itskip) == 0
            # read source field from disk
            cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                if kwargs[:isinterior]
                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
                        wavefields[active_wavefield], ginsu_interior_range)
                else
                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
                        wavefields[active_wavefield])
                end
            end

            # born injection
            cumtime_im += @elapsed WaveFD.forwardBornInjection!(p, kwargs[:modeltype], δm_ginsu, wavefields)
        end

        # extract data at receivers
        cumtime_ex += @elapsed WaveFD.extractdata!(δdinterp, pold, it, iz, ix, c)
    end
    set_zero_subnormals(false)
    JopProp2DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1,
        cumtime_io, cumtime_ex, cumtime_im)

    δd .= 0
    WaveFD.interpforward!(WaveFD.interpfilters(kwargs[:dtmod], kwargs[:dtrec], 0, WaveFD.LangC(), kwargs[:nthreads]), δd, δdinterp)

    for active_wavefield in kwargs[:active_wavefields]
        close(iofields[active_wavefield])
        close(kwargs[:compressor][active_wavefield])
    end

    free(p)

    δd
end

function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractArray; kwargs...)
    nz_ginsu,nx_ginsu = size(kwargs[:ginsu])
    z0_ginsu,x0_ginsu = origin(kwargs[:ginsu])
    p = Prop2DAcoVTIDenQ_DEO2_FDTD(
        nz = nz_ginsu,
        nx = nx_ginsu,
        nbz = kwargs[:nbz_cache],
        nbx = kwargs[:nbx_cache],
        dz = kwargs[:dz],
        dx = kwargs[:dx],
        dt = kwargs[:dtmod],
        nsponge = kwargs[:nsponge],
        nthreads = kwargs[:nthreads],
        freesurface = kwargs[:freesurface],
        freqQ = kwargs[:freqQ],
        qMin = kwargs[:qMin],
        qInterior = kwargs[:qInterior])

    # wave-fields
    pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)

    # ginsu'd earth model
    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))

    # active model-set
    for prop in keys(kwargs[:active_modelset])
        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
    end
    # passive model-set
    for prop in keys(kwargs[:passive_modelset])
        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
    end
    earth_ginsu["f"] .= kwargs[:f]

    δm_ginsu = Dict{String,Array{Float32,2}}()
    for prop in keys(kwargs[:active_modelset])
        δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
    end

    ginsu_interior_range = interior(kwargs[:ginsu])

    # Get receiver interpolation coefficients
    local iz, ix, c
    if kwargs[:interpmethod] == :hicks
        iz, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
    else
        iz, ix, c = WaveFD.linearcoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
    end
    for i = 1:length(c)
        for jx = 1:size(c[i], 2), jz = 1:size(c[i], 1)
            kz = iz[i][jz]
            kx = ix[i][jx]
            c[i][jz,jx] *= kwargs[:dtmod]^2 * earth_ginsu["v"][kz,kx]^2 / earth_ginsu["b"][kz,kx]
        end
    end
    blks = WaveFD.source_blocking(nz_ginsu, nx_ginsu, kwargs[:nbz_inject], kwargs[:nbx_inject], iz, ix, c)

    ntmod = WaveFD.default_ntmod(kwargs[:dtrec], kwargs[:dtmod], kwargs[:ntrec])

    # if necessary, re-run the nonlinear forward
    isvalid, _chksum = isvalid_srcfieldfile(kwargs[:mₒ], kwargs[:srcfieldhost][], kwargs[:srcfieldfile]*"-pold", kwargs[:chksum][])
    if !isvalid
        JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(Array{Float32}(undef,0,0), kwargs[:mₒ]; kwargs...)
        kwargs[:srcfieldhost][] = gethostname()
        kwargs[:chksum][] = _chksum
    end

    # local disk file for source wavefield deserialization
    wavefields = Dict{String,Array{Float32,2}}()
    iofields = Dict{String,IOStream}()
    for active_wavefield in kwargs[:active_wavefields]
        iofields[active_wavefield] = open(kwargs[:srcfieldfile]*"-"*active_wavefield)
        open(kwargs[:compressor][active_wavefield])
        wavefields[active_wavefield] = zeros(Float32, nz_ginsu, nx_ginsu)
    end

    # adjoint interpolation of receiver wavefield
    δdinterp = zeros(Float32, ntmod, size(δd, 2))
    WaveFD.interpadjoint!(WaveFD.interpfilters(kwargs[:dtmod], kwargs[:dtrec], 0, WaveFD.LangC(), kwargs[:nthreads]), δdinterp, δd)

    itskip = round(Int, kwargs[:dtrec]/kwargs[:dtmod])
    kwargs[:reportinterval] == 0 || @info "linear adjoint on $(gethostname()), srcfieldfile=$(kwargs[:srcfieldfile])"
    time1 = time()
    cumtime_io, cumtime_ex, cumtime_im, cumtime_pr = 0.0, 0.0, 0.0, 0.0

    set_zero_subnormals(true)
    for it = ntmod:-1:1
        if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod)
            JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
        end

        # propagate and wavefield swap
        cumtime_pr += @elapsed WaveFD.propagateadjoint!(p)
        pcur,pold = pold,pcur

        # inject receiver data
        cumtime_ex += @elapsed WaveFD.injectdata!(pcur, blks, δdinterp, it)

        if rem(it-1,itskip) == 0
            # read source field from disk
            cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                if kwargs[:isinterior]
                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
                        wavefields[active_wavefield], ginsu_interior_range)
                else
                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
                        wavefields[active_wavefield])
                end
            end

            # born accumulation
            cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
        end
    end
    set_zero_subnormals(false)
    JopProp2DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im)

    # undo ginsu
    for prop in keys(kwargs[:active_modelset])
        super!(@view(δm[:,:,kwargs[:active_modelset][prop]]), kwargs[:ginsu], δm_ginsu[prop], accumulate=false)
    end

    for active_wavefield in kwargs[:active_wavefields]
        close(iofields[active_wavefield])
        close(kwargs[:compressor][active_wavefield])
    end

    free(p)

    δm
end

modelindex(F::Jop{T}, key) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoVTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key]

function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoVTIDenQ_DEO2_FDTD_f!)}}
    s = state(A)
    isvalid, _chksum = isvalid_srcfieldfile(jet(A).mₒ, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[])
    if !isvalid
        JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(Array{Float32}(undef,0,0), m; s...)
        s.chksum[] = _chksum
        s.srcfieldhost[] = gethostname()
    end
    open(s.compressor["pold"])
    I = srcillum!(γ, s.srcfieldfile*"-pold", s.compressor["pold"], s.isinterior, s.ginsu, s.ntrec, s.nthreads)
    close(s.compressor["pold"])
    I
end

function JopProp2DAcoVTIDenQ_DEO2_FDTD_stats(stats, ginsu, it, cumtime_total, cumtime_io, cumtime_ex, cumtime_im=0.0)
    stats["MCells/s"] = megacells_per_second(size(ginsu)..., it-1, cumtime_total)
    stats["%io"] = cumtime_total > 0 ? cumtime_io/cumtime_total*100.0 : 0.0
    stats["%inject/extract"] = cumtime_total > 0 ? cumtime_ex/cumtime_total*100.0 : 0.0
    stats["%imaging"] = cumtime_total > 0 ? cumtime_im/cumtime_total*100.0 : 0.0
end

Jets.perfstat(J::T) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoVTIDenQ_DEO2_FDTD_f!)}} = state(J).stats

@inline function JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_ln(ginsu, it, ntmod, cumtime_total, cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, d::AbstractArray{T}, mode) where {T}
    jt = occursin("adjoint", mode) ? ntmod-it+1 : it
    kt = fmt("5d", jt)
    nt = fmt("5d", ntmod)
    mcells = fmt("7.2f", megacells_per_second(size(ginsu)..., jt-1, cumtime_total))
    IO = fmt("5.2f", cumtime_total > 0 ? cumtime_io/cumtime_total*100.0 : 0.0)
    EX = fmt("5.2f", cumtime_total > 0 ? cumtime_ex/cumtime_total*100.0 : 0.0)
    IM = fmt("5.2f", cumtime_total > 0 ? cumtime_im/cumtime_total*100.0 : 0.0)
    PR = fmt("5.2f", cumtime_total > 0 ? cumtime_pr/cumtime_total*100.0 : 0.0)
    rmsd = fmt("10.4e", length(d) > 0 ? sqrt(norm(d)^2 / length(d)) : zero(T))
    rmsp = fmt("10.4e", sqrt(norm(pcur)^2 / length(pcur)))

    @info "Prop2DAcoVTIDenQ_DEO2_FDTD, $mode, time step $kt of $nt $mcells MCells/s (IO=$IO, EX=$EX, IM=$IM, PR=$PR) -- rms d,p; $rmsd $rmsp"
end

@inline function JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_nl(ginsu, it, ntmod, cumtime_total, cumtime_io, cumtime_ex, cumtime_pr, pcur, d::AbstractArray{T}) where {T}
    kt = fmt("5d", it)
    nt = fmt("5d", ntmod)
    mcells = fmt("7.2f", megacells_per_second(size(ginsu)..., it-1, cumtime_total))
    IO = fmt("5.2f", cumtime_total > 0 ? cumtime_io/cumtime_total*100.0 : 0.0)
    EX = fmt("5.2f", cumtime_total > 0 ? cumtime_ex/cumtime_total*100.0 : 0.0)
    PR = fmt("5.2f", cumtime_total > 0 ? cumtime_pr/cumtime_total*100.0 : 0.0)
    rmsd = fmt("10.4e", length(d) > 0 ? sqrt(norm(d)^2 / length(d)) : zero(T))
    rmsp = fmt("10.4e", sqrt(norm(pcur)^2 / length(pcur)))

    @info "Prop2DAcoVTIDenQ_DEO2_FDTD, nonlinear forward, time step $kt of $nt $mcells MCells/s (IO=$IO, EX=$EX, PR=$PR) -- rms d,p; $rmsd $rmsp"
end

function Base.close(jet::Jet{D,R,typeof(JopProp2DAcoVTIDenQ_DEO2_FDTD_f!)}) where {D,R}
    rm("$(state(jet).srcfieldfile)-pold", force=true)
    rm("$(state(jet).srcfieldfile)-mold", force=true)
    rm("$(state(jet).srcfieldfile)-pspace", force=true)
    rm("$(state(jet).srcfieldfile)-mspace", force=true)
end

max_stars_repo_name: ChevronETC/JetPackWaveFD

max_issues_repo_name: ChevronETC/JetPackWaveFD

max_forks_repo_name: ChevronETC/JetPackWaveFD


text:
using Shearlab, FFTW, Plots
function reorderCoeff(A::AbstractArray{T,2}) where T
    if size(A,1)==1
        includingZeros = reshape(sort(abs.(A), dims =
                                      length(size(A)),rev=true), (length(A,)))
    else
        includingZeros = sort(abs.(A), dims = length(size(A)),rev=true)
    end
    firstZero = findfirst(includingZeros.==0)
    if typeof(firstZero)<:Nothing
        return includingZeros
    else
        return includingZeros[1:firstZero]
    end
end
reorderCoeff(A::AbstractArray{T,1}) where T = sort(abs.(A), rev=true)

N = 32
xlocs = range(-1, 1, length=32)
rot(θ) = [cos(θ) -sin(θ); sin(θ) cos(θ)]
x=1; y=0; θ=π/3; σx = .4;σy = .1
f(xloc,yloc,σx,σy,θ) = [exp.(-1/2*[x,y]'*rot(θ)'*([σx^2 0; 0 σy^2])^(-1)*rot(θ)*[x,y]) for x in xloc, y in yloc ]

gaussianExample = f(xlocs,xlocs,.6,.3,π/6)
heatmap(gaussianExample)

shears = Shearlab.getshearletsystem2D(n, n, 4,typeBecomes=Float32)
size(shears.shearlets)
sheared = Shearlab.sheardec2D(Float32.(gaussianExample), shears)[:,:,1:(end-1)]

mutilated = copy(gaussianExample);
nZeros =floor(Int, 12/32*n); mutilated[1:nZeros,:] = zeros(nZeros,size(mutilated,2))
heatmap(mutilated)

mutsheared = Shearlab.sheardec2D(Float32.(mutilated), shears)[:,:,1:(end-1)]

reorderedShear = reorderCoeff(reshape(sheared,:))
reorderedMutShear = reorderCoeff(reshape(mutsheared,:))
pyplot()
plot(reorderedShear)
plot(log.(1:length(reorderedShear)), max.(-25, log.([reorderedShear reorderedMutShear])), labels=["original", "mutilated"])
savefig("mutilated_Gaussian_Decay_Shearlets_32_4Scales.pdf")


text:
module MixedIntegerExperiments

using JuMP
using Gurobi
using CoinOptServices
using MultilinearOpt
using Polyhedra, CDDLib
using Plots
using DrakeVisualizer
using ColorTypes
using AxisArrayVariables
using AxisArrays
using StaticArrays
using Parameters

include("util.jl")
include("region_descriptions.jl")
include("disjoint_union.jl")
include("mccormick.jl")
include("plotting.jl")
include("shouldbeelsewhere.jl")
include("trajopt.jl")
# include("boxrobotsinterop.jl")

export
    RegionDescription,
    EnvironmentRegion,
    MIQPTrajOptParams,
    ContactPointDescription,
    BoxRobotWithRotation2D,
    ContactPointState,
    BoxRobotWithRotation2DState,
    ContactPointInput,
    axis_aligned_bounding_box,
    contact_region,
    isfree,
    axis_aligned_free_box_region,
    plot_polyhedron!,
    plot_environment,
    plot_allowable_forces,
    plot_kinematic_regions,
    miqp_trajopt,
    rms_constraint_violation

end # module


text:
"""
    M = time_to_mean(t, T, Mo=0., epoch=0.)

Computes the mean anomaly of an orbit with period `T` from a time `t`, the orbit's mean anomaly
at epoch `Mo`, and the epoch.
"""
function time_to_mean(t, T, Mo=0., epoch=0.)
    return Mo + 2π*(t - epoch)/T
end

function time_to_mean(t, orb::Orbit)
    M = time_to_mean(t, orb.period, orb.Mo, orb.epoch)
    if orb.e < 1
        return wrap_angle(M)
    else
        return M
    end
end

"""
    t = mean_to_time(M, T, Mo=0, epoch=0, tmin=nothing)

Computes the time at which an orbit with period `T`, mean anomaly
at epoch `Mo`, and epoch `epoch` has the specified mean anomaly `M`.

Because this time is not uniquely defined for elliptical orbits, `tmin` can be specified
to ensure that the returned time falls between `tmin` and `tmin + T`.
"""
function mean_to_time(M, T, Mo=0, epoch=0, tmin=nothing)
    t = epoch + (M-Mo)*T/(2π)
    if tmin === nothing
        return t
    else
        return t + T*ceil((tmin-t)/T)
    end
end

mean_to_time(M, orb::Orbit, tmin=nothing) = mean_to_time(M, orb.period, orb.Mo, orb.epoch, orb.e<1 ? tmin : nothing)

"""
    θ = mean_to_true(M, e)

Computes the true anomaly of an orbit with eccentricity `e` at the specified mean anomaly `M`.
"""
function mean_to_true(M, e)
    if e == 1       # Parabolic case
        throw(ArgumentError("Parabolic case (e=1) not implemented"))
    elseif e < 1    # Elliptical case
        E = kepler(M, e)
        return 2*atan(√(1+e)*sin(E/2),  √(1-e)*cos(E/2))
    else            # Hyperbolic case
        H = kepler(M, e)
        return 2*atan(√(e+1)*sinh(H/2), √(e-1)*cosh(H/2))
    end
end

"""
    θ = time_to_true(t, e, T, Mo=0., epoch=0.)

Computes the mean anomaly of an orbit with eccentricity `e` and period `T` from a time `t`, the orbit's mean anomaly
at epoch `Mo`, and the epoch.
"""
function time_to_true(t, e, T, Mo=0, epoch=0)
    M = time_to_mean(t, T, Mo, epoch)
    return mean_to_true(M, e)
end

time_to_true(t, orb::Orbit) = time_to_true(t, orb.e, orb.period, orb.Mo, orb.epoch)

"""
    M = true_to_mean(θ, e)

Computes the mean anomaly of an orbit with eccentricity `e` at the specified true anomaly `θ`.
"""
function true_to_mean(θ, e)
    if e == 1       # Parabolic case
        throw(ArgumentError("Parabolic case (e=1) not implemented"))
    elseif e < 1    # Elliptical case
        E = 2*atan(sin(θ/2)*√(1-e), cos(θ/2)*√(1+e))
        return wrap_angle(E - e*sin(E))
    else            # Hyperbolic case
        H = 2*atanh(tan(θ/2)*√((e-1)/(e+1)))
        # H = copysign(acosh((e+cos(θ))/(1+e*cos(θ))), θ)
        return e*sinh(H) - H
    end
end

true_to_mean(θ, orb::Orbit) = true_to_mean(θ, orb.e)

"""
    t = true_to_time(θ, e, T, Mo=0, epoch=0, tmin=nothing)

Computes the time at which an orbit with eccentricity `e`, period `T`, mean anomaly
at epoch `Mo`, and epoch `epoch` has the specified true anomaly `θ`.

Because this time is not uniquely defined for elliptical orbits, `tmin` can be specified
to ensure that the returned time falls between `tmin` and `tmin + T`.
"""
function true_to_time(θ, e, T, Mo=0, epoch=0, tmin=nothing)
    M = true_to_mean(θ, e)
    return mean_to_time(M, T, Mo, epoch, e<1 ? tmin : nothing)
end

# true_to_time(θ, e, a, μ, Mo=0, epoch=0, tmin=nothing) = true_to_time(θ, e, period(a, μ), Mo, epoch, tmin)
true_to_time(θ, orb::Orbit, tmin=nothing) = true_to_time(θ, orb.e, orb.period, orb.Mo, orb.epoch, tmin)


text:
#=
MIT License

Optimal Testing and Containment Strategies for Universities in Mexico amid COVID-19

Copyright © 2021 Test and Contain. Luis Benavides-Vázquez, Héctor Alonso Guzmán-Gutiérrez, Jakob Jonnerby, Philip Lazos, Edwin Lock, Francisco J. Marmolejo-Cossío, Ninad Rajgopal,and José Roberto Tello-Ayala. https://www.testandcontain.com/

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#
include("frontier.jl")

function count(pop, config)
    k = length(pop.p)
    budget, step, guaranteed = config.budget, config.stepsize, config.guaranteed
    gsize, bsizes = config.groupsizes, config.bucketsizes
    solutions = Dict()
    s = [0 for i in 1:2k]  # preallocate memory for solution
    @showprogress for groups in GroupIterator(gsize, k)
        for tests in TestIterator(budget, step, guaranteed, pop.n .÷ groups, k)
            s[1:k] .= tests  # test portion
            s[k+1:2k] = [tests[i]===0 ? 0 : groups[i] for i in 1:k]  # group portion
            solutions[copy(s)] = true
        end
    end
    return length(solutions)
end

"""Count number of solutions, loading parameters from `filename` file."""
function compute(filename)
    pop, config = load_from_json(filename)
    return count(pop, config)
end

function run()
    # Process commandline arguments
    if length(ARGS) === 0
        println("Please specify the input data file.")
        println("\nUsage: run `julia path/to/data.json` in the terminal.")
        return
    end
    input = ARGS[1]

    # Compute solutions
    println("Number of solutions: ", compute(input))
    return
end

run()

max_stars_repo_name: hguzmang/eaamo

max_issues_repo_name: hguzmang/eaamo

max_forks_repo_name: hguzmang/eaamo


text:
using Documenter, JPkg

makedocs(;
    modules=[JPkg],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/mattg3004/JPkg.jl/blob/{commit}{path}#L{line}",
    sitename="JPkg.jl",
    authors="Matthew",
    assets=String[],
)

deploydocs(;
    repo="github.com/mattg3004/JPkg.jl",
)


text:
"""
    LogSpaced(min::Float64, max::Float64, n::Int64)

This function evaluates ``n`` points, logarithmically spaced between
    ``min`` and ``max``.
"""
function LogSpaced(min::Float64, max::Float64, n::Int64)
    logmin = log10(min)
    logmax = log10(max)
    logarray = Array(LinRange(logmin, logmax, n))
    return exp10.(logarray)
end

"""
    ΔLogSpaced(min::Float64, max::Float64, n::Int64)

This function evaluates ``n`` points, logarithmically spaced between
    ``min`` and ``max``.
"""
function ΔLogSpaced(min::Float64, max::Float64, n::Int64)
    λmin = log10(min)
    λmax = log10(max)
    Δℓ = zeros(n)
    for i in 1:n
        Δℓ[i] = 10^(λmin+i*(λmax-λmin)/n)-10^(λmin+(i-1)*(λmax-λmin)/n)
    end
    return Δℓ
end

"""
    BinSearch(x::Float64, Array::Vector{Float64})

Given a value ``z`` and an Array, determines the couple of array elements where
``z`` lies and returns the index corresponding to the first value.
"""
function BinSearch(x::Float64, Array::Vector{Float64})
    idx = 1
    if x <= first(Array)
        idx = 1
    elseif x >= last(Array)
        idx = length(Array)-1
    else
        while !(x >= Array[idx] && x <= Array[idx+1])
            idx += 1
        end
    end
    return idx
end

"""
    CustomRegression(x::Vector{Float64}, y::Vector{Float64})

Given two arrays ``x`` and ``y``, performs the linear regression and returns the
coefficients ``c`` and ``m`` of the fitted line.
"""
function CustomRegression(x::Vector{Float64}, y::Vector{Float64})
    x_mean = mean(x)
    y_mean = mean(y)
    m = 0
    den = 0
    for i in 1:length(x)
        m += (x[i]-x_mean)*(y[i]-y_mean)
        den  += (x[i]-x_mean)^2
    end
    m /= den
    c = y_mean - m* x_mean
    return c, m
end

"""
    SimpsonWeightArray(n::Int64)

This function evaluates an array with the Simpson weight, for Simpson
Integration
"""
function SimpsonWeightArray(n::Int64)
    number_intervals = floor((n-1)/2)
    weight_array = zeros(n)
    if n == number_intervals*2+1
        for i in 1:number_intervals
            weight_array[Int((i-1)*2+1)] += 1/3
            weight_array[Int((i-1)*2+2)] += 4/3
            weight_array[Int((i-1)*2+3)] += 1/3
        end
    else
        weight_array[1] += 0.5
        weight_array[2] += 0.5
        for i in 1:number_intervals
            weight_array[Int((i-1)*2+1)+1] += 1/3
            weight_array[Int((i-1)*2+2)+1] += 4/3
            weight_array[Int((i-1)*2+3)+1] += 1/3
        end
        weight_array[length(weight_array)]   += 0.5
        weight_array[length(weight_array)-1] += 0.5
        for i in 1:number_intervals
            weight_array[Int((i-1)*2+1)] += 1/3
            weight_array[Int((i-1)*2+2)] += 4/3
            weight_array[Int((i-1)*2+3)] += 1/3
        end
        weight_array ./= 2
    end
    return weight_array
end

"""
    SimpsonWeightMatrix(n::Int64)

This function evaluates an array with the Simpson weight, for Simpson
Integration of the Lensing Efficiency
"""
function SimpsonWeightMatrix(n::Int64)
    number_intervals = floor((n-1)/2)
    weight_matrix = zeros(n, n)
    for i in 1:n
        weight_matrix[i,i:n] = SimpsonWeightArray(n-i+1)
    end
    return weight_matrix
end

"""
    Difference(InputArray::Vector{Float64})

This function evaluates the n-th discrete difference of a given 1-D array.
"""
function Difference(InputArray::Vector{Float64})
    OutputArray = zeros(length(InputArray)-1)
    for i in 1:length(InputArray)-1
        OutputArray[i] = InputArray[i+1] - InputArray[i]
    end
    return OutputArray
end

"""
    UnevenTrapzWeightArray(InputArray::Vector{Float64})

This function evaluates the array of weights for the integration with uneven
spaced points and the Trapezoidal rule.
"""
function UnevenTrapzWeightArray(InputArray::Vector{Float64})
    WeightArray = zeros(length(InputArray))
    DifferenceArray = Difference(InputArray)
    for idx in 2:length(WeightArray)-1
        WeightArray[idx] = (DifferenceArray[idx]+DifferenceArray[idx-1])/2
    end
    WeightArray[1] = DifferenceArray[1]/2
    WeightArray[length(InputArray)] = DifferenceArray[length(InputArray)-1]/2
    return WeightArray
end

function UnevenTrapzWeightMatrix(InputMatrix::AbstractArray{Float64, 2})
    WeightMatrix = zeros(size(InputMatrix))
    for lidx in 1:length(WeightMatrix[:,1])
        WeightMatrix[lidx,:] = UnevenTrapzWeightArray(InputMatrix[lidx,:])
    end
    return WeightMatrix
end

function uᵢⱼ(n::Int, i::Int, j::Int)
    u = zeros(floor(Int,0.5*n*(n+1)))
    u[floor(Int64, (j-1)*n+i-0.5*j*(j-1))] = 1
    return u
end

function eᵢ(n::Int, i::Int)
    e = zeros(n)
    e[i] = 1
    return e
end

function Eᵢⱼ(n::Int, i::Int, j::Int)
    E = zeros((n,n))
    E[i,j] = 1
    return E
end

function Tᵢⱼ(n::Int, i::Int, j::Int)
    if i != j
        T = Eᵢⱼ(n,i,j) + Eᵢⱼ(n,j,i)
    else
        T = Eᵢⱼ(n,i,i)
    end
    return T
end

function EliminationMatrix(n::Int)
    L = zeros((floor(Int64,0.5*n*(n+1)), n*n))
    for i in 1:n
        for j in 1:i
            L += kron(uᵢⱼ(n,i,j), transpose(vec(Eᵢⱼ(n,i,j))))
        end
    end
    return L
end

"""
    DuplicationMatrix(n::Int)

Duplication matrix ``\\boldsymbol{D}_n``.
"""
function DuplicationMatrix(n::Int)
    D = zeros((n*n, floor(Int64,0.5*n*(n+1))))
    for i in 1:n
        for j in 1:i
            D += transpose(kron(uᵢⱼ(n,i,j), transpose(vec(Tᵢⱼ(n,i,j)))))
        end
    end
    return D
end


text:
floattype(T::Type{<:Integer}) = float(T)
floattype(T::Type{<:AbstractFloat}) = T
floattype(_) = Float64

Base.@propagate_inbounds @inline Base.getindex(v::AbstractVector, ::typeof(!), i) = v[i]
Base.@propagate_inbounds @inline Base.getindex(v::AbstractVector, ::typeof(!), i::AbstractRange) = view(v, i)
Base.@propagate_inbounds @inline Base.getindex(v::AbstractMatrix, ::typeof(!), i) = view(v,:,i)
Base.@propagate_inbounds @inline Base.getindex(v::AbstractArray, ::typeof(!), i) = view(v, Base.front(ntuple(i -> Colon(), ndims(v)))..., i)

Base.@propagate_inbounds @inline Base.setindex!(v::AbstractVector, val, ::typeof(!), i) = v[i] = val
Base.@propagate_inbounds @inline Base.setindex!(v::AbstractVector, val, ::typeof(!), i::AbstractRange) = v[i] .= val
Base.@propagate_inbounds @inline Base.setindex!(v::AbstractMatrix, val, ::typeof(!), i) = v[:,i] .= val
Base.@propagate_inbounds @inline Base.setindex!(v::AbstractArray, val, ::typeof(!), i) = view(v, Base.front(ntuple(i -> Colon(), ndims(v)))..., i) .= val

"Return the length of the last axis"
lastlength(x) = size(x, ndims(x))
actuallastlength(x) = lastlength(x)

"""
    getwindow(T, m, i)

Extract a window of length `m` staring at `i`. The returned object will be an `UnsafeView` into `T`.
"""
Base.@propagate_inbounds @inline getwindow(T,m,i) = T[!, (0:m-1) .+ i]


text:
using Documenter
using VLConstraintBasedModelGenerationUtilities

makedocs(
    sitename = "Model",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [VLConstraintBasedModelGenerationUtilities],
    pages = [
        "Home" => "index.md",
    ],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/varnerlab/VLConstraintBasedModelGenerationUtilities.jl.git",
    devurl = "stable",
    devbranch = "main",
)


text:
################################################################################
## Linear Regression
##   y ~ N(b0 + b1 * x, s2)
##   b0, b1 ~ N(0, 1000)
##   s2 ~ invgamma(0.001, 0.001)
################################################################################

using Mamba

## Data
data = Dict(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)

## Log-transformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector
logfgrad = function(x::DenseVector)
  b0 = x[1]
  b1 = x[2]
  logs2 = x[3]
  r = data[:y] - b0 - b1 * data[:x]
  logf = (-0.5 * length(data[:y]) - 0.001) * logs2 -
           (0.5 * dot(r, r) + 0.001) / exp(logs2) -
           0.5 * b0^2 / 1000 - 0.5 * b1^2 / 1000
  grad = [
    sum(r) / exp(logs2) - b0 / 1000,
    sum(data[:x] .* r) / exp(logs2) - b1 / 1000,
    -0.5 * length(data[:y]) - 0.001 + (0.5 * dot(r, r) + 0.001) / exp(logs2)
  ]
  logf, grad
end

## MCMC Simulation with Metropolis-Adjusted Langevin Algorithm
## Without (1) and with (2) a user-specified proposal covariance matrix
n = 5000
sim1 = Chains(n, 3, names = ["b0", "b1", "s2"])
sim2 = Chains(n, 3, names = ["b0", "b1", "s2"])
epsilon = 0.1
Sigma = eye(3)
theta1 = MALAVariate([0.0, 0.0, 0.0], epsilon, logfgrad)
theta2 = MALAVariate([0.0, 0.0, 0.0], epsilon, Sigma, logfgrad)
for i in 1:n
  sample!(theta1)
  sample!(theta2)
  sim1[i, :, 1] = [theta1[1:2]; exp(theta1[3])]
  sim2[i, :, 1] = [theta2[1:2]; exp(theta2[3])]
end
describe(sim1)
describe(sim2)


text:
# ops.jl - overloaded operators for constucting symbolic expressions.
#
# A couple of examples:
#
#    :x ⊕ :y      ==>  :(x + y)
#    2 ⊗ :(x ⊕ y) ==>  :(2 * (x + y))

# import Base: +, -, *, /, .+, .-, .*, ./

# ⊕(ex::Symbolic, v::Numeric) = :($ex + $v)
# ⊕(v::Numeric, ex::Symbolic) = :($v + $ex)
# ⊕(ex1::Symbolic, ex2::Symbolic) = :($ex1 + $ex2)
# ⊕(x, y) = x + y

# (-)(ex::Symbolic, v::Numeric) = :($ex - $v)
# (-)(v::Numeric, ex::Symbolic) = :($v - $ex)
# (-)(ex1::Symbolic, ex2::Symbolic) = :($ex1 - $ex2)

# ⊗(ex::Symbolic, v::Numeric) = :($ex * $v)
# ⊗(v::Numeric, ex::Symbolic) = :($v * $ex)
# ⊗(ex1::Symbolic, ex2::Symbolic) = :($ex1 * $ex2)
# ⊗(x, y) = x * y

# (/)(ex::Symbolic, v::Numeric) = :($ex / $v)
# (/)(v::Numeric, ex::Symbolic) = :($v / $ex)
# (/)(ex1::Symbolic, ex2::Symbolic) = :($ex1 / $ex2)

# elementwise operations

⊕(ex::Symbolic, v::Numeric) = :($ex .+ $v)
⊕(v::Numeric, ex::Symbolic) = :($v .+ $ex)
⊕(ex1::Symbolic, ex2::Symbolic) = :($ex1 .+ $ex2)
⊕(x, y) = x .+ y

# (.-)(ex::Symbolic, v::Numeric) = :($ex .- $v)
# (.-)(v::Numeric, ex::Symbolic) = :($v .- $ex)
# (.-)(ex1::Symbolic, ex2::Symbolic) = :($ex1 .- $ex2)

⊗(ex::Symbolic, v::Numeric) = :($ex .* $v)
⊗(v::Numeric, ex::Symbolic) = :($v .* $ex)
⊗(ex1::Symbolic, ex2::Symbolic) = :($ex1 .* $ex2)
⊗(x, y) = x .* y

# (./)(ex::Symbolic, v::Numeric) = :($ex ./ $v)
# (./)(v::Numeric, ex::Symbolic) = :($v ./ $ex)
# (./)(ex1::Symbolic, ex2::Symbolic) = :($ex1 ./ $ex2)

## Note: this is also correct, but generates more ugly expressions
##
## for op in [:(+), :(-), :(*), :(/)]
##     for T in [Number, Array]
##         @eval begin
##             Base.$(op)(ex::Symbolic, v::$T) = :(($$op)($ex,$v))
##             Base.$(op)(v::$T, ex::Symbolic) = :(($$op)($v,$ex))
##             Base.$(op)(ex1::Symbolic, ex2::Symbolic) = :(($$op($ex1,$ex2)))
##         end
##     end
## end


text:
const HUMAN1_PUBLICATION_PROJ_NAME = "Human1_Publication_Data_Scripts"
const HUMAN1_PROBLICATION_ZENODO_LINK = "https://zenodo.org/record/3583004/files/Human1_Publication_Data_Scripts.zip?download=1"

# Some usful iders
const BIOMASS_IDER = "biomass_human"
const COST_IDER = "prot_pool_exchange"
const ATPM_IDER = "HMR_3964";
const EXCH_SUBSYS_HINT = "Exchange/demand"

# Human1
const CELL_NAMES = ["HS_578T", "RPMI_8226", "HT29", "MALME_3M", "SR",
               "UO_31", "MDMAMB_231", "HOP62", "NCI_H226", "HOP92", "O_786"];
const CONSTRAINTS_NAMES = ["media", "glucose", "L-lactate", "threonine"];
const CONSTRAINTS_LEVELS = [0, 1, 2, 3]

# Common constants
const MAX_ABS_BOUND = 100
const MAX_CONC = 99999 # Inf conc means that the metabolite will be never limiting the growth
const ZEROTH = 1e-8

const MODEL_TYPE_KEY = :model
const MODEL_TYPE_SUFFIX = "__Model"
const EC_MODEL_TYPE_KEY = :ecModel
const EC_MODEL_TYPE_SUFFIX = "__ecModel"


text:
import Optimize
using Krylov
using OptimizationProblems

# Nonconvex problems var ≥ 10 (50 pbs)
Problems = [broydn7d, brybnd, chainwoo, chnrosnb_mod, cosine, cragglvy, curly,
dixmaane, dixmaanf, dixmaang, dixmaanh, dixmaani, dixmaanj, dixmaank, dixmaanl,
dixmaanm, dixmaann, dixmaano, dixmaanp, edensch, eg2, errinros_mod, extrosnb,
fletcbv2, fletcbv3_mod, fletchcr, fminsrf2, freuroth, genhumps, genrose,
genrose_nash, indef_mod, liarwhd, morebv, ncb20, ncb20b, noncvxu2, noncvxun,
nondia, penalty2, powellsg, schmvett, scosine, sinquad, sparsine, sparsqur,
srosenbr, tointgss, tquartic, woods]

# Convex problems var ≥ 10 (17 pbs)
# Problems = [arglina, arglinb, arglinc, arwhead, bdqrtic, clplatea, clplateb, clplatec,
# dixon3dq, dqdrtic, dqrtic, engval1, nondquar, power, quartc, tridia, vardim]

n = 100

function OptPbs()
    solvers = Dict{Symbol, Function}(:CR => LCR,:CG => LCG)
    bmark_args = Dict{Symbol, Any}(:skipif => model -> !unconstrained(model))

    Optimize.bmark_solvers(solvers,
                (MathProgNLPModel(p(n), name=string(p)) for p in Problems);
                bmark_args...)

end

stats = OptPbs()
tf = font(10) # titlefont
f = font(10)

pb_type = "nonconvex problems"
p = performance_profile(hcat([p[:, 1] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("Performance profile: #f in the solution of ", pb_type),
                            titlefont = tf, legendfont = f, guidefont = f,
                            legend=:bottomright, palette=:blues) # Profile for #f
Plots.xlabel!("Within this factor of the best (log₂ scale)")
Plots.savefig(p, string("profil_f_", pb_type, ".pdf"))

p = performance_profile(hcat([p[:, 2] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("Performance profile: #g in the solution of ", pb_type),
                            titlefont = tf, legendfont = f, guidefont = f,
                            legend=:bottomright, palette=:blues) # Profile for #g
Plots.xlabel!("Within this factor of the best (log₂ scale)")
Plots.savefig(p, string("profil_g_", pb_type, ".pdf"))

p = performance_profile(hcat([p[:, 3] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("Performance profile: #Hv in the solution of ", pb_type),
                            titlefont = tf, legendfont = f, guidefont = f,
                            legend=:bottomright, palette=:blues) # Profile for #Hv
Plots.xlabel!("Within this factor of the best (log₂ scale)")
Plots.savefig(p, string("profil_Hv_", pb_type, ".pdf"))

p = performance_profile(hcat([sum(p, 2) for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("Performance profile: #f + #g  + #Hv in the solution of ", pb_type),
                            titlefont = tf, legendfont = f, guidefont = f,
                            legend=:bottomright, palette=:blues) # Profile for #f + #g + #Hv
Plots.xlabel!("Within this factor of the best (log₂ scale)")
Plots.savefig(p, string("profil_f+g+Hv_", pb_type, ".pdf"))


text:
const MOIU = MOI.Utilities

const AFF = MOI.VectorAffineFunction{Float64}

mutable struct Solution
    ret_val::Union{Nothing,Int}
    primal::Vector{Float64}
    dual_eq::Vector{Float64}
    dual_ineq::Vector{Float64}
    slack::Vector{Float64}
    objective_value::Float64
    dual_objective_value::Float64
    solve_time::Float64
    iter::Int
end
# The values used by ECOS are from -7 to 10 so -10 should be safe
function Solution()
    return Solution(
        nothing,
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        NaN,
        NaN,
        NaN,
        0,
    )
end

MOIU.@product_of_sets(Zeros, MOI.Zeros)
MOIU.@product_of_sets(
    Cones,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    PermutedExponentialCone
)

MOIU.@struct_of_constraints_by_set_types(
    ZerosOrNot,
    MOI.Zeros,
    Union{MOI.Nonnegatives,MOI.SecondOrderCone,PermutedExponentialCone},
)

const OptimizerCache = MOI.Utilities.GenericModel{
    Cdouble,
    MOIU.ObjectiveContainer{Cdouble},
    MOIU.VariablesContainer{Cdouble},
    ZerosOrNot{Cdouble}{
        MOIU.MatrixOfConstraints{
            Cdouble,
            MOIU.MutableSparseMatrixCSC{
                Cdouble,
                Clong,
                MOI.Utilities.ZeroBasedIndexing,
            },
            Vector{Cdouble},
            Zeros{Cdouble},
        },
        MOIU.MatrixOfConstraints{
            Cdouble,
            MOIU.MutableSparseMatrixCSC{
                Cdouble,
                Clong,
                MOI.Utilities.ZeroBasedIndexing,
            },
            Vector{Cdouble},
            Cones{Cdouble},
        },
    },
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    zeros::Union{Nothing,Zeros{Cdouble}}
    cones::Union{Nothing,Cones{Cdouble}}
    sol::Solution
    silent::Bool
    options::Dict{String,Any}
    function Optimizer(; kwargs...)
        optimizer = new(nothing, nothing, Solution(), false, Dict{String,Any}())
        if length(kwargs) > 0
            @warn(
                "Passing keyword attributes is deprecated. Use " *
                "`set_optimizer_attribute` instead.",
            )
        end
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawOptimizerAttribute(String(key)), value)
        end
        return optimizer
    end
end

function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges)
    return [PermutedExponentialBridge{Cdouble}]
end

function _rows(optimizer::Optimizer, ci::MOI.ConstraintIndex{AFF,MOI.Zeros})
    return MOIU.rows(optimizer.zeros, ci)
end
function _rows(optimizer::Optimizer, ci::MOI.ConstraintIndex{AFF})
    return MOIU.rows(optimizer.cones, ci)
end

MOI.get(::Optimizer, ::MOI.SolverName) = "ECOS"

function MOI.supports(::Optimizer, param::MOI.RawOptimizerAttribute)
    return hasfield(Csettings, Symbol(param.name))
end

function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
    optimizer.options[param.name] = value
    return
end

function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
    # TODO: This gives a poor error message if the name of the parameter is
    # invalid.
    return optimizer.options[param.name]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    return optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

function MOI.is_empty(optimizer::Optimizer)
    return optimizer.zeros === nothing && optimizer.cones === nothing
end

function MOI.empty!(optimizer::Optimizer)
    optimizer.zeros = nothing
    optimizer.cones = nothing
    optimizer.sol = Solution()
    return
end

function MOI.supports(
    ::Optimizer,
    ::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    },
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{AFF},
    ::Type{
        <:Union{
            MOI.Zeros,
            MOI.Nonnegatives,
            MOI.SecondOrderCone,
            PermutedExponentialCone,
        },
    },
)
    return true
end

function _optimize!(dest::Optimizer, src::OptimizerCache)
    MOI.empty!(dest)
    Ab = MOI.Utilities.constraints(src.constraints, AFF, MOI.Zeros)
    A = Ab.coefficients
    Gh = MOI.Utilities.constraints(src.constraints, AFF, MOI.Nonnegatives)
    G = Gh.coefficients
    q =
        MOI.dimension.(
            MOI.get.(
                Gh,
                MOI.ConstraintSet(),
                MOI.get(
                    Gh,
                    MOI.ListOfConstraintIndices{AFF,MOI.SecondOrderCone}(),
                ),
            ),
        )
    @assert A.n == G.n
    max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    obj =
        MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    objective_constant = MOI.constant(obj)
    c0 = zeros(A.n)
    for term in obj.terms
        c0[term.variable.value] += term.coefficient
    end

    dest.zeros = deepcopy(Ab.sets) # TODO copy(Ab.sets)
    dest.cones = deepcopy(Gh.sets) # TODO copy(Gh.sets)

    options = Dict(Symbol(k) => v for (k, v) in dest.options)
    if dest.silent
        options[:verbose] = false
    end

    inner = _setup(
        A.n,
        G.m,
        A.m,
        Gh.sets.num_rows[1],
        length(q),
        q,
        MOI.get(Gh, MOI.NumberOfConstraints{AFF,PermutedExponentialCone}()),
        ECOSMatrix(-G.nzval, G.colptr, G.rowval),
        ECOSMatrix(-A.nzval, A.colptr, A.rowval),
        max_sense ? -c0 : c0,
        Gh.constants,
        Ab.constants,
    )
    settings(inner, options)
    ret_val = ECOS.solve(inner)
    ecos_prob = unsafe_load(inner)
    stat = unsafe_load(ecos_prob.info)
    solve_time = stat.tsetup + stat.tsolve
    iter = stat.iter
    primal = unsafe_wrap(Array, ecos_prob.x, ecos_prob.n)[:]
    dual_eq = unsafe_wrap(Array, ecos_prob.y, ecos_prob.p)[:]
    dual_ineq = unsafe_wrap(Array, ecos_prob.z, ecos_prob.m)[:]
    slack = unsafe_wrap(Array, ecos_prob.s, ecos_prob.m)[:]
    ECOS.cleanup(inner, 0)
    objective_value = (max_sense ? -1 : 1) * stat.pcost
    dual_objective_value = (max_sense ? -1 : 1) * stat.dcost
    dest.sol = Solution(
        ret_val,
        primal,
        dual_eq,
        dual_ineq,
        slack,
        objective_value,
        dual_objective_value,
        solve_time,
        iter,
    )
    if !MOIU.is_ray(MOI.get(dest, MOI.PrimalStatus()))
        dest.sol.objective_value += objective_constant
    end
    if !MOIU.is_ray(MOI.get(dest, MOI.DualStatus()))
        dest.sol.dual_objective_value += objective_constant
    end
    return
end

function MOI.optimize!(dest::Optimizer, src::OptimizerCache)
    _optimize!(dest, src)
    return MOIU.identity_index_map(src), false
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    cache = OptimizerCache()
    index_map = MOI.copy_to(cache, src)
    _optimize!(dest, cache)
    return index_map, false
end

MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec) = optimizer.sol.solve_time
function MOI.get(optimizer::Optimizer, ::MOI.BarrierIterations)
    return Int64(optimizer.sol.iter)
end
function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    # Strings from https://github.com/ifa-ethz/ecos/blob/master/include/ecos.h
    flag = optimizer.sol.ret_val
    if flag === nothing
        return "Optimize not called"
    elseif flag == ECOS_OPTIMAL
        return "Problem solved to optimality"
    elseif flag == ECOS_OPTIMAL + ECOS_INACC_OFFSET
        return "Problem solved to inaccurate optimality"
    elseif flag == ECOS_PINF
        return "Found certificate of primal infeasibility"
    elseif flag == ECOS_PINF + ECOS_INACC_OFFSET
        return "Found inaccurate certificate of primal infeasibility"
    elseif flag == ECOS_DINF
        return "Found certificate of dual infeasibility"
    elseif flag == ECOS_DINF + ECOS_INACC_OFFSET
        return "Found inaccurate certificate of dual infeasibility"
    elseif flag == ECOS_MAXIT
        return "Maximum number of iterations reached"
    elseif flag == ECOS_NUMERICS
        return "Search direction unreliable"
    elseif flag == ECOS_OUTCONE
        return "s or z got outside the cone, numerics?"
    elseif flag == ECOS_SIGINT
        return "solver interrupted by a signal/ctrl-c"
    else
        @assert flag == ECOS_FATAL
        return "Unknown problem in solver"
    end
end

# Implements getter for result value and statuses
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    flag = optimizer.sol.ret_val
    if flag === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    elseif flag == ECOS.ECOS_OPTIMAL
        return MOI.OPTIMAL
    elseif flag == ECOS.ECOS_PINF
        return MOI.INFEASIBLE
    elseif flag == ECOS.ECOS_DINF
        return MOI.DUAL_INFEASIBLE
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.ITERATION_LIMIT
    elseif flag == ECOS.ECOS_NUMERICS || flag == ECOS.ECOS_OUTCONE
        return MOI.NUMERICAL_ERROR
    elseif flag == ECOS.ECOS_SIGINT
        return MOI.INTERRUPTED
    elseif flag == ECOS.ECOS_FATAL
        return MOI.OTHER_ERROR
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.ALMOST_OPTIMAL
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.ALMOST_INFEASIBLE
    elseif flag == ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
        return MOI.ALMOST_DUAL_INFEASIBLE
    else
        error("Unrecognized ECOS solve status flag: $flag.")
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.objective_value
end
function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.dual_objective_value
end

function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    flag = optimizer.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif flag == ECOS.ECOS_PINF
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS.ECOS_DINF
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.UNKNOWN_RESULT_STATUS
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.NEARLY_FEASIBLE_POINT
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    else
        return MOI.OTHER_RESULT_STATUS
    end
end
function MOI.get(
    optimizer::Optimizer,
    attr::MOI.VariablePrimal,
    vi::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.primal[vi.value]
end
function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{AFF,MOI.Zeros},
)
    MOI.check_result_index_bounds(optimizer, attr)
    return zeros(length(_rows(optimizer, ci)))
end
function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex,
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.slack[_rows(optimizer, ci)]
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    flag = optimizer.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif flag == ECOS.ECOS_PINF
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS.ECOS_DINF
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.UNKNOWN_RESULT_STATUS
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.NEARLY_FEASIBLE_POINT
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
        return MOI.INFEASIBLE_POINT
    else
        return MOI.OTHER_RESULT_STATUS
    end
end
function _dual(optimizer, ci::MOI.ConstraintIndex{AFF,MOI.Zeros})
    return optimizer.sol.dual_eq
end
_dual(optimizer, ci::MOI.ConstraintIndex) = optimizer.sol.dual_ineq
function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{AFF},
)
    MOI.check_result_index_bounds(optimizer, attr)
    return _dual(optimizer, ci)[_rows(optimizer, ci)]
end

MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1


text:
using BinDeps

@BinDeps.setup

polylib = library_dependency("libpolylib", aliases=["libpolylib64-8"])