text:
"""
    ReducedSpaceEvaluator{T, VI, VT, MT, Jacx, Jacu, JacCons, Hess} <: AbstractNLPEvaluator

Reduced-space evaluator projecting the optimization problem
into the powerflow manifold defined by the nonlinear equation ``g(x, u) = 0``.
The state ``x`` is defined implicitly, as a function of the control
``u``. Hence, the powerflow equation is implicitly satisfied
when we are using this evaluator.

Once a new point `u` is passed to the evaluator,
the user needs to call the method `update!` to find the corresponding
state ``x(u)`` satisfying the balance equation ``g(x(u), u) = 0``.

Taking as input a [`PolarForm`](@ref) structure, the reduced evaluator
builds the bounds corresponding to the control `u`,
The reduced evaluator could be instantiated on the host memory, or on a specific device
(currently, only CUDA is supported).

## Examples

```julia-repl
julia> datafile = "case9.m"  # specify a path to a MATPOWER instance
julia> nlp = ReducedSpaceEvaluator(datafile)
A ReducedSpaceEvaluator object
    * device: KernelAbstractions.CPU()
    * #vars: 5
    * #cons: 10
    * constraints:
        - voltage_magnitude_constraints
        - active_power_constraints
        - reactive_power_constraints
    * linear solver: ExaPF.LinearSolvers.DirectSolver()
```

If a GPU is available, we could instantiate `nlp` as

```julia-repl
julia> nlp_gpu = ReducedSpaceEvaluator(datafile; device=CUDADevice())
A ReducedSpaceEvaluator object
    * device: KernelAbstractions.CUDADevice()
    * #vars: 5
    * #cons: 10
    * constraints:
        - voltage_magnitude_constraints
        - active_power_constraints
        - reactive_power_constraints
    * linear solver: ExaPF.LinearSolvers.DirectSolver()

```

## Note
Mathematically, we set apart the state ``x`` from the control ``u``,
and use a third variable ``y`` --- the by-product --- to denote the remaining
values of the network.
In the implementation of `ReducedSpaceEvaluator`,
we only deal with a control `u` and an attribute `buffer`,
storing all the physical values needed to describe the network.
The attribute `buffer` stores the values of the control `u`, the state `x`
and the by-product `y`. Each time we are calling the method `update!`,
the values of the control are copied into the buffer.

"""
mutable struct ReducedSpaceEvaluator{T, VI, VT, MT, Jacx, Jacu, JacCons, Hess} <: AbstractNLPEvaluator
    model::PolarForm{T, VI, VT, MT}
    λ::VT

    u_min::VT
    u_max::VT

    constraints::Vector{Function}
    g_min::VT
    g_max::VT

    # Cache
    buffer::PolarNetworkState{VI, VT}
    # AutoDiff
    state_jacobian::FullSpaceJacobian{Jacx, Jacu}
    obj_stack::AutoDiff.TapeMemory{typeof(active_power_generation), AdjointStackObjective{VT}, Nothing}
    cons_stacks::Vector{AutoDiff.TapeMemory} # / constraints
    constraint_jacobians::JacCons
    hessians::Hess

    # Options
    linear_solver::LinearSolvers.AbstractLinearSolver
    powerflow_solver::AbstractNonLinearSolver
    has_jacobian::Bool
    update_jacobian::Bool
    has_hessian::Bool
end

function ReducedSpaceEvaluator(
    model::PolarForm{T, VI, VT, MT};
    constraints=Function[voltage_magnitude_constraints, active_power_constraints, reactive_power_constraints],
    linear_solver=direct_linear_solver(model),
    powerflow_solver=NewtonRaphson(tol=1e-12),
    want_jacobian=true,
    want_hessian=true,
) where {T, VI, VT, MT}
    # First, build up a network buffer
    buffer = get(model, PhysicalState())
    # Populate buffer with default values of the network, as stored
    # inside model
    init_buffer!(model, buffer)

    u_min, u_max = bounds(model, Control())
    λ = similar(buffer.dx)

    g_min = VT()
    g_max = VT()
    for cons in constraints
        cb, cu = bounds(model, cons)
        append!(g_min, cb)
        append!(g_max, cu)
    end

    obj_ad = pullback_objective(model)
    state_ad = FullSpaceJacobian(model, power_balance)
    cons_ad = AutoDiff.TapeMemory[]
    for cons in constraints
        push!(cons_ad, AutoDiff.TapeMemory(model, cons, VT))
    end

    # Jacobians
    cons_jac = nothing
    if want_jacobian
        cons_jac = ConstraintsJacobianStorage(model, constraints)
    end

    # Hessians
    hess_ad = nothing
    if want_hessian
        hess_ad = HessianStorage(model, constraints)
    end

    return ReducedSpaceEvaluator(
        model, λ, u_min, u_max,
        constraints, g_min, g_max,
        buffer,
        state_ad, obj_ad, cons_ad, cons_jac, hess_ad,
        linear_solver, powerflow_solver, want_jacobian, false, want_hessian,
    )
end
function ReducedSpaceEvaluator(datafile::String; device=KA.CPU(), options...)
    return ReducedSpaceEvaluator(PolarForm(datafile, device); options...)
end

array_type(nlp::ReducedSpaceEvaluator) = array_type(nlp.model)

n_variables(nlp::ReducedSpaceEvaluator) = length(nlp.u_min)
n_constraints(nlp::ReducedSpaceEvaluator) = length(nlp.g_min)

constraints_type(::ReducedSpaceEvaluator) = :inequality
has_hessian(::ReducedSpaceEvaluator) = true

# Getters
get(nlp::ReducedSpaceEvaluator, ::Constraints) = nlp.constraints
function get(nlp::ReducedSpaceEvaluator, ::State)
    x = similar(nlp.λ) ; fill!(x, 0)
    get!(nlp.model, State(), x, nlp.buffer)
    return x
end
get(nlp::ReducedSpaceEvaluator, ::PhysicalState) = nlp.buffer
get(nlp::ReducedSpaceEvaluator, ::AutoDiffBackend) = nlp.autodiff
# Physics
get(nlp::ReducedSpaceEvaluator, ::PS.VoltageMagnitude) = nlp.buffer.vmag
get(nlp::ReducedSpaceEvaluator, ::PS.VoltageAngle) = nlp.buffer.vang
get(nlp::ReducedSpaceEvaluator, ::PS.ActivePower) = nlp.buffer.pg
get(nlp::ReducedSpaceEvaluator, ::PS.ReactivePower) = nlp.buffer.qg
function get(nlp::ReducedSpaceEvaluator, attr::PS.AbstractNetworkAttribute)
    return get(nlp.model, attr)
end

# Setters
function setvalues!(nlp::ReducedSpaceEvaluator, attr::PS.AbstractNetworkValues, values)
    setvalues!(nlp.model, attr, values)
end
function setvalues!(nlp::ReducedSpaceEvaluator, attr::PS.ActiveLoad, values)
    setvalues!(nlp.model, attr, values)
    setvalues!(nlp.buffer, attr, values)
end
function setvalues!(nlp::ReducedSpaceEvaluator, attr::PS.ReactiveLoad, values)
    setvalues!(nlp.model, attr, values)
    setvalues!(nlp.buffer, attr, values)
end

# Transfer network values inside buffer
function transfer!(
    nlp::ReducedSpaceEvaluator, vm, va, pg, qg,
)
    setvalues!(nlp.buffer, PS.VoltageMagnitude(), vm)
    setvalues!(nlp.buffer, PS.VoltageAngle(), va)
    setvalues!(nlp.buffer, PS.ActivePower(), pg)
    setvalues!(nlp.buffer, PS.ReactivePower(), qg)
end

# Initial position
function initial(nlp::ReducedSpaceEvaluator)
    u = similar(nlp.u_min) ; fill!(u, 0.0)
    return get!(nlp.model, Control(), u, nlp.buffer)
end

# Bounds
bounds(nlp::ReducedSpaceEvaluator, ::Variables) = (nlp.u_min, nlp.u_max)
bounds(nlp::ReducedSpaceEvaluator, ::Constraints) = (nlp.g_min, nlp.g_max)

## Callbacks
function update!(nlp::ReducedSpaceEvaluator, u)
    jac_x = nlp.state_jacobian.x
    # Transfer control u into the network cache
    transfer!(nlp.model, nlp.buffer, u)
    # Get corresponding point on the manifold
    conv = powerflow(
        nlp.model,
        jac_x,
        nlp.buffer,
        nlp.powerflow_solver;
        linear_solver=nlp.linear_solver
    )

    if !conv.has_converged
        error("Newton-Raphson algorithm failed to converge ($(conv.norm_residuals))")
        return conv
    end

    # Refresh values of active and reactive powers at generators
    update!(nlp.model, PS.Generators(), PS.ActivePower(), nlp.buffer)
    # Evaluate Jacobian of power flow equation on current u
    AutoDiff.jacobian!(nlp.model, nlp.state_jacobian.u, nlp.buffer)
    # Specify that constraint's Jacobian is not up to date
    nlp.update_jacobian = nlp.has_jacobian
    return conv
end

# TODO: determine if we should include λ' * g(x, u), even if ≈ 0
function objective(nlp::ReducedSpaceEvaluator, u)
    # Take as input the current cache, updated previously in `update!`.
    return cost_production(nlp.model, nlp.buffer)
end

function constraint!(nlp::ReducedSpaceEvaluator, g, u)
    ϕ = nlp.buffer
    mf = 1::Int
    mt = 0::Int
    for cons in nlp.constraints
        m_ = size_constraint(nlp.model, cons)::Int
        mt += m_
        cons_ = @view(g[mf:mt])
        cons(nlp.model, cons_, ϕ)
        mf += m_
    end
end

function _forward_solve!(nlp::ReducedSpaceEvaluator, y, x)
    if isa(y, CUDA.CuArray)
        ∇gₓ = nlp.state_jacobian.x.J
        LinearSolvers.ldiv!(nlp.linear_solver, y, ∇gₓ, x)
    else
        LinearSolvers.ldiv!(nlp.linear_solver, y, x)
    end
end

function _backward_solve!(nlp::ReducedSpaceEvaluator, y::VT, x::VT) where {VT <: AbstractArray}
    ∇gₓ = nlp.state_jacobian.x.J
    if isa(nlp.linear_solver, LinearSolvers.AbstractIterativeLinearSolver)
        # Iterative solver case
        ∇gT = LinearSolvers.get_transpose(nlp.linear_solver, ∇gₓ)
        # Switch preconditioner to transpose mode
        LinearSolvers.update!(nlp.linear_solver, ∇gT)
        # Compute adjoint and store value inside λₖ
        LinearSolvers.ldiv!(nlp.linear_solver, y, ∇gT, x)
    elseif isa(y, CUDA.CuArray)
        ∇gT = LinearSolvers.get_transpose(nlp.linear_solver, ∇gₓ)
        LinearSolvers.ldiv!(nlp.linear_solver, y, ∇gT, x)
    else
        LinearSolvers.rdiv!(nlp.linear_solver, y, x)
    end
end

###
# First-order code
####
#
# compute inplace reduced gradient (g = ∇fᵤ + (∇gᵤ')*λ)
# equivalent to: g = ∇fᵤ - (∇gᵤ')*λ_neg
# (take λₖ_neg to avoid computing an intermediate array)
function reduced_gradient!(
    nlp::ReducedSpaceEvaluator, grad, ∂fₓ, ∂fᵤ, λ, u,
)
    ∇gᵤ = nlp.state_jacobian.u.J
    ∇gₓ = nlp.state_jacobian.x.J

    # λ = ∇gₓ' \ ∂fₓ
    _backward_solve!(nlp, λ, ∂fₓ)

    grad .= ∂fᵤ
    mul!(grad, transpose(∇gᵤ), λ, -1.0, 1.0)
    return
end
function reduced_gradient!(nlp::ReducedSpaceEvaluator, grad, ∂fₓ, ∂fᵤ, u)
    reduced_gradient!(nlp::ReducedSpaceEvaluator, grad, ∂fₓ, ∂fᵤ, nlp.λ, u)
end

# Compute only full gradient wrt x and u
function full_gradient!(nlp::ReducedSpaceEvaluator, gx, gu, u)
    buffer = nlp.buffer
    ∂obj = nlp.obj_stack
    # Evaluate adjoint of cost function and update inplace AdjointStackObjective
    gradient_objective!(nlp.model, ∂obj, buffer)
    copyto!(gx, ∂obj.stack.∇fₓ)
    copyto!(gu, ∂obj.stack.∇fᵤ)
end

function gradient!(nlp::ReducedSpaceEvaluator, g, u)
    buffer = nlp.buffer
    # Evaluate adjoint of cost function and update inplace AdjointStackObjective
    gradient_objective!(nlp.model, nlp.obj_stack, buffer)
    ∇fₓ, ∇fᵤ = nlp.obj_stack.stack.∇fₓ, nlp.obj_stack.stack.∇fᵤ

    reduced_gradient!(nlp, g, ∇fₓ, ∇fᵤ, u)
    return
end

function jacobian_structure(nlp::ReducedSpaceEvaluator)
    m, n = n_constraints(nlp), n_variables(nlp)
    nnzj = m * n
    rows = zeros(Int, nnzj)
    cols = zeros(Int, nnzj)
    jacobian_structure!(nlp, rows, cols)
    return rows, cols
end

function jacobian_structure!(nlp::ReducedSpaceEvaluator, rows, cols)
    m, n = n_constraints(nlp), n_variables(nlp)
    idx = 1
    for i in 1:n # number of variables
        for c in 1:m #number of constraints
            rows[idx] = c ; cols[idx] = i
            idx += 1
        end
    end
end

function _update_full_jacobian_constraints!(nlp)
    if nlp.update_jacobian
        update_full_jacobian!(nlp.model, nlp.constraint_jacobians, nlp.buffer)
        nlp.update_jacobian = false
    end
end

# Works only on the CPU!
function jacobian!(nlp::ReducedSpaceEvaluator, J, u)
    m, n = n_constraints(nlp), n_variables(nlp)
    jac = reshape(J, m, n)

    _update_full_jacobian_constraints!(nlp)

    ∇cons = nlp.constraint_jacobians
    Jx = ∇cons.Jx
    Ju = ∇cons.Ju
    m, nₓ = size(Jx)
    m, nᵤ = size(Ju)
    ∇gᵤ = nlp.state_jacobian.u.J
    ∇gₓ = nlp.state_jacobian.x.J
    # Compute state sensitivities all in once
    μ = zeros(nₓ, nᵤ)
    LinearSolvers.ldiv!(nlp.linear_solver, μ, ∇gᵤ)
    # Compute reduced Jacobian
    copy!(jac, Ju)
    mul!(jac, Jx, μ, -1.0, 1.0)
    return
end

function jprod!(nlp::ReducedSpaceEvaluator, jv, u, v)
    nᵤ = length(u)
    m  = n_constraints(nlp)
    @assert nᵤ == length(v)

    _update_full_jacobian_constraints!(nlp)

    ∇cons = nlp.constraint_jacobians

    Jx = ∇cons.Jx
    Ju = ∇cons.Ju

    ∇gᵤ = nlp.state_jacobian.u.J
    rhs = nlp.buffer.dx
    z = nlp.buffer.balance
    # init RHS
    mul!(rhs, ∇gᵤ, v)
    # Compute z
    _forward_solve!(nlp, z, rhs)

    # jv .= Ju * v .- Jx * z
    mul!(jv, Ju, v)
    mul!(jv, Jx, z, -1.0, 1.0)
    return
end

function full_jtprod!(nlp::ReducedSpaceEvaluator, jvx, jvu, u, v)
    fr_ = 0::Int
    for (cons, stack) in zip(nlp.constraints, nlp.cons_stacks)
        n = size_constraint(nlp.model, cons)::Int
        mask = fr_+1:fr_+n
        vv = @view v[mask]
        # Compute jtprod of current constraint
        jacobian_transpose_product!(nlp.model, stack, nlp.buffer, vv)
        jvx .+= stack.stack.∂x
        jvu .+= stack.stack.∂u
        fr_ += n
    end
end

function jtprod!(nlp::ReducedSpaceEvaluator, jv, u, v)
    ∂obj = nlp.obj_stack
    μ = nlp.buffer.balance
    jvx = ∂obj.stack.jvₓ ; fill!(jvx, 0)
    jvu = ∂obj.stack.jvᵤ ; fill!(jvu, 0)
    full_jtprod!(nlp, jvx, jvu, u, v)
    reduced_gradient!(nlp, jv, jvx, jvu, μ, u)
end

function ojtprod!(nlp::ReducedSpaceEvaluator, jv, u, σ, v)
    ∂obj = nlp.obj_stack
    jvx = ∂obj.stack.jvₓ ; fill!(jvx, 0)
    jvu = ∂obj.stack.jvᵤ ; fill!(jvu, 0)
    # compute gradient of objective
    full_gradient!(nlp, jvx, jvu, u)
    jvx .*= σ
    jvu .*= σ
    # compute transpose Jacobian vector product of constraints
    full_jtprod!(nlp, jvx, jvu, u, v)
    # Evaluate gradient in reduced space
    reduced_gradient!(nlp, jv, jvx, jvu, u)
    return
end

###
# Second-order code
####
# z = -(∇gₓ  \ (∇gᵤ * w))
function _second_order_adjoint_z!(
    nlp::ReducedSpaceEvaluator, z, w,
)
    ∇gᵤ = nlp.state_jacobian.u.J
    rhs = nlp.buffer.dx
    mul!(rhs, ∇gᵤ, w, -1.0, 0.0)
    _forward_solve!(nlp, z, rhs)
end

# ψ = -(∇gₓ' \ (∇²fₓₓ .+ ∇²gₓₓ))
function _second_order_adjoint_ψ!(
    nlp::ReducedSpaceEvaluator, ψ, ∂fₓ,
)
    _backward_solve!(nlp, ψ, ∂fₓ)
end

function _reduced_hessian_prod!(
    nlp::ReducedSpaceEvaluator, hessvec, ∂fₓ, ∂fᵤ, tgt,
)
    nx = get(nlp.model, NumberOfState())
    nu = get(nlp.model, NumberOfControl())
    H = nlp.hessians
    ∇gᵤ = nlp.state_jacobian.u.J
    ψ = H.ψ

    hv = H.tmp_hv

    ## POWER BALANCE HESSIAN
    AutoDiff.adj_hessian_prod!(nlp.model, H.state, hv, nlp.buffer, nlp.λ, tgt)
    ∂fₓ .-= @view hv[1:nx]
    ∂fᵤ .-= @view hv[nx+1:nx+nu]

    # Second order adjoint
    _second_order_adjoint_ψ!(nlp, ψ, ∂fₓ)

    hessvec .+= ∂fᵤ
    mul!(hessvec, transpose(∇gᵤ), ψ, -1.0, 1.0)
    return
end

function hessprod!(nlp::ReducedSpaceEvaluator, hessvec, u, w)
    @assert nlp.hessians != nothing

    nx = get(nlp.model, NumberOfState())
    nu = get(nlp.model, NumberOfControl())
    buffer = nlp.buffer
    H = nlp.hessians

    fill!(hessvec, 0.0)

    # Two vector products
    tgt = H.tmp_tgt
    hv = H.tmp_hv
    z = H.z

    _second_order_adjoint_z!(nlp, z, w)

    # Init tangent
    copyto!(tgt, 1, z, 1, nx)
    copyto!(tgt, nx+1, w, 1, nu)

    ∂f = similar(buffer.pg)
    ∂²f = similar(buffer.pg)

    # Adjoint and Hessian of cost function
    adjoint_cost!(nlp.model, ∂f, buffer.pg)
    hessian_cost!(nlp.model, ∂²f)

    ## OBJECTIVE HESSIAN
    hessian_prod_objective!(nlp.model, H.obj, nlp.obj_stack, hv, ∂²f, ∂f, buffer, tgt)
    ∇²fx = hv[1:nx]
    ∇²fu = hv[nx+1:nx+nu]

    _reduced_hessian_prod!(nlp, hessvec, ∇²fx, ∇²fu, tgt)

    return
end

function hessian_lagrangian_penalty_prod!(
    nlp::ReducedSpaceEvaluator, hessvec, u, y, σ, w, D,
)
    @assert nlp.hessians != nothing

    nx = get(nlp.model, NumberOfState())
    nu = get(nlp.model, NumberOfControl())
    buffer = nlp.buffer
    H = nlp.hessians

    fill!(hessvec, 0.0)

    z = H.z
    _second_order_adjoint_z!(nlp, z, w)

    # Two vector products
    tgt = H.tmp_tgt
    hv = H.tmp_hv

    # Init tangent
    tgt[1:nx] .= z
    tgt[1+nx:nx+nu] .= w

    ∂f = similar(buffer.pg)
    ∂²f = similar(buffer.pg)
    # Adjoint and Hessian of cost function
    adjoint_cost!(nlp.model, ∂f, buffer.pg)
    hessian_cost!(nlp.model, ∂²f)

    ## OBJECTIVE HESSIAN
    hessian_prod_objective!(nlp.model, H.obj, nlp.obj_stack, hv, ∂²f, ∂f, buffer, tgt)
    ∇²Lx = σ .* @view hv[1:nx]
    ∇²Lu = σ .* @view hv[nx+1:nx+nu]

    # CONSTRAINT HESSIAN
    shift = 0
    hvx = @view hv[1:nx]
    hvu = @view hv[nx+1:nx+nu]
    for (cons, Hc) in zip(nlp.constraints, H.constraints)
        m = size_constraint(nlp.model, cons)::Int
        mask = shift+1:shift+m
        yc = @view y[mask]
        AutoDiff.adj_hessian_prod!(nlp.model, Hc, hv, buffer, yc, tgt)
        ∇²Lx .+= hvx
        ∇²Lu .+= hvu
        shift += m
    end
    # Add Hessian of quadratic penalty
    diagjac = similar(y)
    if !iszero(D)
        _update_full_jacobian_constraints!(nlp)
        Jx = nlp.constraint_jacobians.Jx
        Ju = nlp.constraint_jacobians.Ju
        # ∇²Lx .+= Jx' * (D * (Jx * z)) .+ Jx' * (D * (Ju * w))
        # ∇²Lu .+= Ju' * (D * (Jx * z)) .+ Ju' * (D * (Ju * w))
        mul!(diagjac, Jx, z)
        mul!(diagjac, Ju, w, 1.0, 1.0)
        diagjac .*= D
        mul!(∇²Lx, Jx', diagjac, 1.0, 1.0)
        mul!(∇²Lu, Ju', diagjac, 1.0, 1.0)
    end

    # Second order adjoint
    _reduced_hessian_prod!(nlp, hessvec, ∇²Lx, ∇²Lu, tgt)

    return
end

# Return lower-triangular matrix
function hessian_structure(nlp::ReducedSpaceEvaluator)
    n = n_variables(nlp)
    rows = Int[r for r in 1:n for c in 1:r]
    cols = Int[c for r in 1:n for c in 1:r]
    return rows, cols
end

# Utils function
function primal_infeasibility!(nlp::ReducedSpaceEvaluator, cons, u)
    constraint!(nlp, cons, u) # Evaluate constraints
    (n_inf, err_inf, n_sup, err_sup) = _check(cons, nlp.g_min, nlp.g_max)
    return max(err_inf, err_sup)
end
function primal_infeasibility(nlp::ReducedSpaceEvaluator, u)
    cons = similar(nlp.g_min) ; fill!(cons, 0)
    return primal_infeasibility!(nlp, cons, u)
end

# Printing
function sanity_check(nlp::ReducedSpaceEvaluator, u, cons)
    println("Check violation of constraints")
    print("Control  \t")
    (n_inf, err_inf, n_sup, err_sup) = _check(u, nlp.u_min, nlp.u_max)
    @printf("UB: %.4e (%d)    LB: %.4e (%d)\n",
            err_sup, n_sup, err_inf, n_inf)
    print("Constraints\t")
    (n_inf, err_inf, n_sup, err_sup) = _check(cons, nlp.g_min, nlp.g_max)
    @printf("UB: %.4e (%d)    LB: %.4e (%d)\n",
            err_sup, n_sup, err_inf, n_inf)
end

function Base.show(io::IO, nlp::ReducedSpaceEvaluator)
    n = n_variables(nlp)
    m = n_constraints(nlp)
    println(io, "A ReducedSpaceEvaluator object")
    println(io, "    * device: ", nlp.model.device)
    println(io, "    * #vars: ", n)
    println(io, "    * #cons: ", m)
    println(io, "    * constraints:")
    for cons in nlp.constraints
        println(io, "        - ", cons)
    end
    print(io, "    * linear solver: ", nlp.linear_solver)
end

function reset!(nlp::ReducedSpaceEvaluator)
    # Reset adjoint
    fill!(nlp.λ, 0)
    # Reset buffer
    init_buffer!(nlp.model, nlp.buffer)
    return
end


text:
# sorting in data.frame is slow
if false
    N = 100_000_000
    K = 100
    using DataFrames
    srand(1)
    @time df = DataFrame(idstr = rand([@sprintf "id%03d" k for k in 1:(N/K)], N)
        , id = rand(1:K, N)
        , val = rand(1:5,N))

    @time sort!(df,cols=[:id, :val])

    using uCSV
    @time uCSV.write("df.csv", df);

    using CSV
    @time CSV.write("df.csv", df);
end

ext: jl


text:
mutable struct Chunk{T, MaxSize, TypeSize}
    csize ::Int64
    nsize ::Int64
    buff ::Vector{UInt8}
end


text:
"""
# Usage

    solve(α, dx, dt, xStart, xEnd, n, κ, FiniteDiffIm())

Use implicit finite difference to discrete equations.

### References

```tex
@article{Murio2008ImplicitFD,
  title={Implicit finite difference approximation for time fractional diffusion equations},
  author={Diego A. Murio},
  journal={Comput. Math. Appl.},
  year={2008},
  volume={56},
  pages={1138-1145}
}
```

Matlab version: https://github.com/awstown/Fractional-Derivative
"""
struct FiniteDiffIm <: FractionalDiffEqAlgorithm end

function solve(α, dx, dt, xStart, xEnd, n, κ, u0t, uendt, u0, ::FiniteDiffIm)
    x = collect(0:dx:xEnd)
    t = collect(0:dt:n)

    U = zeros(Int64(n/dt + 1), Int64((xEnd - xStart)/dx + 1))

    # Boundry conditions
    U[:, 1] .= u0t
    U[:, end] .= uendt

    U[1, :] .= u0.(x)

    mu = (dt^α)/(dx^2)
    r = mu * gamma(2-α)
    A1 = diagm((1+2*r)*ones(length(x)-2))
    A2 = diagm(1 => -r*ones(length(x)-3))
    A3 = diagm(-1 => -r*ones(length(x)-3))
    A = A1 + A2 + A3

    U[2, 2:end-1] = A \ U[1, 2:end-1]

    j = collect(0:length(t))
    b_j = (j.+1).^(1-α) .- j.^(1-α)
    c_j = b_j[1:end-1] - b_j[2:end]

    V = copy(U[2:end, 2:end-1])
    for k = 1:length(V[:, 1])-1
        V[k+1, :] = A \ ImNextStep(V, U[1, 2:end-1], k, b_j, c_j)'
    end
    U[2:end, 2:end-1] = V
    return U
end

ImNextStep( V, U_0, k, b_j, c_j ) = c_j[1:k]' * reverse(V[1:k, :], dims=1) .+ (b_j[k+1].*U_0)'


text:
using Test
using QuantumOptics
using QuantumOptics.polynomials

@testset "polynomials" begin

# Test Horner scheme
c = [0.2, 0.6, 1.7]
x0 = 1.3
@test horner(c, x0) == c[1] + c[2]*x0 + c[3]*x0^2

# Test Hermite polynomials
an = Vector{Vector{Int}}(undef, 8)
an[1] = [1]
an[2] = [0,2]
an[3] = [-2, 0, 4]
an[4] = [0, -12, 0, 8]
an[5] = [12, 0, -48, 0, 16]
an[6] = [0, 120, 0, -160, 0, 32]
an[7] = [-120, 0, 720, 0, -480, 0, 64]
an[8] = [0, -1680, 0, 3360, 0, -1344, 0, 128]
@test hermite.a(7) == an

A = hermite.A(7)
for n=0:7
    @test A[n+1] ≈ an[n+1]/sqrt(2^n*factorial(n))
end

end # testset


text:
module Mpmm
using ...Utils, ...Types
using DifferentialEquations, LinearAlgebra, Distances
#using aeroCoeff
include("aerodynamic_coeff.jl")
export trajectoryMPMM!, QEfinderMPMM!, iniCond

function mpmm(du,u,p,t)
    Iₓ,w_bar, d, R, g₀,ω_bar, m,target,dist = p

    vx= u[4]
    vy=u[5]
    vz=u[6]
    X₁ = u[1]
    X₂ = u[2]
    X₃ = u[3]
    pₛ = u[7]
   u_bar =[u[4], u[5], u[6]]
   v_bar= u_bar - w_bar
    v = norm(v_bar)
    du_bar = [du[4], du[5],du[6]]
    #println(X₂)
    ρ =  Utils.density(X₂)
    #println("ρ", " ", ρ)
    α = atan(vz/vx)
    β = asin(vy/v)
    #αₜ = asin(sqrt((sin(α)*cos(β))^2+(sin(β))^2))
    #αₜ = asin()
    #println("v", " ", v)
    #println("v", " ", v_bar)
    #println("u", " ", u_bar)
    mach = Utils.machNumber(v,X₂)
    CMα₀ = CM_α0_inter(mach)
    CMα₂ = CM_α2_inter(mach)
    CMα = CMα₀ #+ CMα₂*(sin(αₜ))^2
    #println(du_bar)

   αₑ_bar = yawOfRepose(v_bar, du_bar,pₛ, t, Iₓ, ρ,d,CMα,v)
   #println("αₑ", " ", αₑ_bar)
    αₑ = norm(αₑ_bar)
    αₜ = asin(αₑ)
    #println("αₜ", " ", αₜ )
    #println("αₜ2", " ", sin(αₑ))

    CD₀ = CD0_inter(mach)
    CDα² = CDδ2_inter(mach)
    CDα⁴ = 0.0
    CLα = Cl_α0_inter(mach)
    CLα³ = Cl_α2_inter(mach)
    CLα⁵ = 0.0
    Cmag = CN_pα_inter(mach,(rad2deg(αₜ))^2)
    #println("mach", " ", mach)
    #println("αₜ", " ", αₜ)
    Cₛₚᵢₙ = Cl_p_inter(mach)

    #CD = CD₀ + CDα²*αₑ^2+CDα⁴*αₑ^4
    CD = CD₀ + CDα² * (sin(αₜ))^2
    #CL = CLα + CLα³*αₑ^2+CLα⁵*αₑ^4
    CL = CLα + CLα³ * (sin(αₜ))^2
    #println("CL", " ", CL)
    FDx = - pi*ρ*d^2/(8)*(CD)*v*vx #drag
    FDy = - pi*ρ*d^2/(8)*(CD)*v*vy #drag
    FDz = - pi*ρ*d^2/(8)*(CD)*v*vz #drag
    FLx = pi*ρ*d^2/(8)*(CL)*v^2*αₑ_bar[1] #lift
    FLy = pi*ρ*d^2/(8)*(CL)*v^2*αₑ_bar[2] #lift
    FLz = pi*ρ*d^2/(8)*(CL)*v^2*αₑ_bar[3] #lift
    FMx = -pi*ρ*d^3*pₛ*Cmag/8*cross(αₑ_bar,v_bar)[1]
    FMy = -pi*ρ*d^3*pₛ*Cmag/8*cross(αₑ_bar,v_bar)[2]
    FMz = -pi*ρ*d^3*pₛ*Cmag/8*cross(αₑ_bar,v_bar)[3]
    #println("cmag", " ", Cmag)
    FGx = - g₀*X₁/R*m
    FGy = -g₀*(1-2*X₂/R)*m
    FGz = - g₀*X₃/R*m
    #FGz = 0.0
    FCx = -2*cross(ω_bar,u_bar)[1]*m
    FCy = -2*cross(ω_bar,u_bar)[2]*m
    FCz = -2*cross(ω_bar,u_bar)[3]*m
    #println("FDx", " ", FDx)
    #println("FLx", " ", FLx)
    #println("FMx", " ", FMx)
    #println("FGx", " ", FGx)
    #println("FCx", " ", FCx)
    Fx =  FDx + FLx +FGx +FMx+FCx
    Fy =  FDy + FLy +FGy +FMy+FCy
    Fz =  FDz + FLz +FGz +FMz+FCz
    #println("Fx", " ", Fx)
    #println("FDx", " ", FDx)
    #println("FLx", " ", FLx)
    #println("FMx", " ", FMx)
    #println("FMy", " ", FMy)
    #println("FMz", " ", FMz)
    #println("FGx", " ", FGx)
    #println("FCx", " ", FCx)
    du[1] = u[4] #x
    du[2] = u[5] #y
    du[3] = u[6] #z
    du[4] =Fx/m #ax
    du[5] = Fy/m#ay
    du[6] = Fz/m#az
    du[7] = pi*ρ*d^4*pₛ*v*Cₛₚᵢₙ/(8*Iₓ)#pdot
    global yaw = αₑ_bar

end

function yawOfRepose(v_bar::Array{Float64,1}, du_bar::Array{Float64,1},pₛ::Float64, t::Float64, Iₓ::Float64, ρ::Float64,d::Float64,CMα::Float64,v::Float64)
    if t==0
        return [0.0,0.0,0.0]
    else
       αₑ_bar = -8*Iₓ*pₛ*cross(v_bar,du_bar)/(pi*ρ*d^3*(CMα)*v^4)

        return αₑ_bar
    end
end

#u₀ = norm()

#canon = gun(0.0,0.0,0.0)

#sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6)
#sol = solve(prob, AutoTsit5(Rosenbrock23()))
#sol = solve(prob)
function condition(u,t,integrator) # Event when event_f(u,t) == 0
  #u[2]
  #target.x-u[1]
  #dist-u[1]
  #delta = sqrt((u[1]-integrator.p[8][1])^2+(u[2]-integrator.p[8][2])^2+(u[3]-integrator.p[8][3])^2)-integrator.p[9]
  #println(delta)
  #return delta
  #toto = (euclidean([u[1],u[2],u[3]],integrator.p[8])-integrator.p[9])
  #println(toto)
  integrator.p[8][1]-u[1]
  #integrator.p[8][1]>u[1]
end

function timeFuze(u,t,integrator)
    #euclidean([u[1],u[2],u[3]],integrator.p[8])-integrator.p[9]

    #close=(euclidean([u[1],u[2],u[3]],integrator.p[8])-integrator.p[9])<0.0
    close=(euclidean([u[1],u[2],u[3]],integrator.p[8])-integrator.p[9])-1e-6
    #println((euclidean([u[1],u[2],u[3]],integrator.p[8])-integrator.p[9]))
    #println(close)
    return close
end

affect!(integrator) = terminate!(integrator)
affect_tf!(integrator) = terminate!(integrator)

function trajectoryMPMM!(u0::Array{Float64,1}, tspan::Tuple{Float64,Float64}, p, proj::AbstractPenetrator, target::Target2D)
    #global dist = target.position[1]
    #append!(p,dist)
    prob = ODEProblem(mpmm,u0,tspan, p)
    #cb = DiscreteCallback(condition,affect!)
#println(target.x)
    cb = ContinuousCallback(condition,affect!)
    cb_tf = ContinuousCallback(timeFuze,affect!)
    #cb_tf=DiscreteCallback(timeFuze,affect_tf!)
    cbs = CallbackSet(cb,cb_tf)
    sol = solve(prob,Tsit5(),callback=cbs, reltol=1e-8, abstol=1e-8)
#sol = solve(prob,Tsit5(),callback=cb)
#sol = solve(prob,BS5(),callback=cb)
#sol = solve(prob,AutoTsit5(Rosenbrock23()),callback=cb)
 proj.position[1] = sol.u[end][1]
 proj.position[2] = sol.u[end][2]
 proj.position[3] = sol.u[end][3]
 proj.velocity[1] = sol.u[end][4]
 proj.velocity[2] = sol.u[end][5]
 proj.velocity[3] = sol.u[end][6]
 proj.tof = sol.t[end]
 proj.spin = sol.u[end][7]
 return yaw
end

function iniCond(QE::Float64, ΔAZ::Float64, u₀::Float64, gun::Canon, calibre::Float64)
    p₀ = 2*pi*u₀/(gun.tc*calibre)
    u₀_bar = [u₀*cosd(QE)*cosd(ΔAZ), u₀*sind(QE), u₀*cosd(QE)*sind(ΔAZ)]
    X₀_bar = [gun.lw*cosd(QE)*cosd(ΔAZ), gun.X2w + gun.lw *sind(QE), gun.lw*cosd(QE)*sind(ΔAZ)]
    u0 = [X₀_bar[1], X₀_bar[2], X₀_bar[3], u₀_bar[1], u₀_bar[2], u₀_bar[3],p₀]
    return u0
end

function QEfinderMPMM!(drone::Target2D, proj::Projectile2D, u₀::Float64, g₀::Float64, gun::Canon, w_bar::Array{Float64,1},lat::Float64)

    epsilonAz = 1e6
    epsilonQE = 1e6
    precisie = 0.001
    ddoel = euclidean(proj.position, drone.position)
    tdoel = sqrt(drone.position[1]^2+drone.position[3]^2)/u₀
    QE = (drone.position[2] - proj.position[2] + g₀ /2 *tdoel^2)*tdoel/u₀
    AZ = 0.0
    tspan = (0.0,1000.0)
    R = 6.356766*1e6 #m
    Ω = 7.292115*1e-5 #rad/s
    ω_bar = [Ω*cosd(lat)*cosd(AZ), Ω*sind(lat), -Ω*cosd(lat)*sind(AZ)]
    #p = [proj.inertia[1],w_bar, proj.calibre, R, g₀, ω_bar, proj.mass,drone]
    #αₑ_bar = [0.0,0.0,0.0]
    p = [proj.inertia[1],w_bar, proj.calibre, R, g₀, ω_bar, proj.mass,drone.position,0.0]

    while abs(epsilonAz)>precisie || abs(epsilonQE)>precisie
        u0 = iniCond(QE, AZ, u₀, gun, proj.calibre)
        #global proj = projectile(u0[1],u0[2],u0[3],u0[4],u0[5],u0[6],0.0)
        #proj = projectile(u0[1],u0[2],u0[3],u0[4],u0[5],u0[6],0.0)
        proj.position = [u0[1],u0[2],u0[3]]
        proj.velocity = [u0[4],u0[5],u0[6]]
        proj.tof = 0.0
        trajectoryMPMM!(u0, tspan, p, proj, drone)

        #global epsilonQE = proj.y - drone.y
        epsilonQE = proj.position[2] - drone.position[2]
        #println("epsilonQE", " ", epsilonQE)
        #QE = QE0 + (accuracy)/(range_/QE0)
        #global epsilonAz = (sqrt((proj.z-drone.z)^2+(proj.x-drone.x)^2)*sign(atan(proj.z)/proj.x)-atan(drone.z/drone.x))
        epsilonAz = (sqrt((proj.position[3]-drone.position[3])^2+(proj.position[1]-drone.position[1])^2)*sign(atan(proj.position[3])/proj.position[1])-atan(drone.position[3]/drone.position[1]))
        #println("epsilonAz", " ", epsilonAz)
        #global AZ = AZ - epsilonAz/ddoel
        AZ = AZ - epsilonAz/ddoel
        #global QE = QE - epsilonQE/ddoel
        QE = QE - epsilonQE/ddoel

        #trajectory!(u0, tspan, p, proj, drone)
        #calcRange = euclidean([0.0,0.0,0.0], [proj.x,proj.y, proj.z])
        #global accuracy = range_  - calcRange
        #println("AZ", " ", AZ)
        #println("QE", " ", QE)

    end

return QE,AZ
end

end


text:
using Test, HammerScalpel

function solveModelTest()

    # Solve model with standard parameters
    @time noint    = nopolicy(TT=2);                     # no intervention at all
    @time notest   = withpolicy("notest", TT=2, ndims=3);         # no testing
    @time untarget = withpolicy("untargettest", TT=2, ndims=3);   # untargeted testing
    @time target   = withpolicy("targettest", TT=2, ndims=3);     # targeted testing
    @time isolate  = withpolicy("isolate", TT=2, ndims=3);        # isolation

    return("Success")

end

@test solveModelTest()=="Success"


text:
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using TestEnv

Pkg.activate(dirname(@__DIR__))
Pkg.instantiate()

TestEnv.activate()


text:
export FSimGate, fsim_block

"""
    FSimGate{T<:Number} <: PrimitiveBlock{2}

The two parameter `FSim` gate.

References
-------------------------
* Arute, Frank, et al. "Quantum supremacy using a programmable superconducting processor." Nature 574.7779 (2019): 505-510.
"""
mutable struct FSimGate{T<:Number} <: PrimitiveBlock{2}
    theta::T
    phi::T
end
YaoBlocks.nqudits(fs::FSimGate) = 2

function Base.:(==)(fs1::FSimGate, fs2::FSimGate)
    return fs1.theta == fs2.theta && fs1.phi == fs2.phi
end

function YaoAPI.mat(::Type{T}, fs::FSimGate) where T
    θ, ϕ = fs.theta, fs.phi
    T[1 0          0          0;
     0 cos(θ)     -im*sin(θ) 0;
     0 -im*sin(θ) cos(θ)     0;
     0 0          0          exp(-im*ϕ)]
end

YaoAPI.iparams_eltype(::FSimGate{T}) where T = T
YaoAPI.getiparams(fs::FSimGate{T}) where T = (fs.theta, fs.phi)
function YaoAPI.setiparams!(fs::FSimGate{T}, θ, ϕ) where T
    fs.theta = θ
    fs.phi = ϕ
    return fs
end

YaoBlocks.@dumpload_fallback FSimGate FSimGate
YaoBlocks.Optimise.to_basictypes(fs::FSimGate) = fsim_block(fs.theta, fs.phi)

"""
    fsim_block(θ::Real, ϕ::Real)

The circuit representation of FSim gate.
"""
function fsim_block(θ::Real, ϕ::Real)
    if θ ≈ π/2
        return cphase(2,2,1,-ϕ)*SWAP*rot(kron(Z,Z), -π/2)*put(2,1=>phase(-π/4))
    else
        return cphase(2,2,1,-ϕ)*rot(SWAP,2*θ)*rot(kron(Z,Z), -θ)*put(2,1=>phase(θ/2))
    end
end


text:
# Similar to JuMP/test/solvers.jl

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

mos = false # try_import(:MathOptInterfaceMosek)
csd = try_import(:CSDP)
scs = false # try_import(:SCS)

# See https://github.com/JuliaOpt/JuMP.jl/pull/1396
function with_opt(mod::Module, args...; kwargs...)
    JuMP.OptimizerFactory(mod.Optimizer, args, kwargs)
end

# Semidefinite factories
sdp_factories = JuMP.OptimizerFactory[]
# Need at least Mosek 8 for sosdemo3 to pass; see:
# https://github.com/JuliaOpt/Mosek.jl/issues/92
# and at least Mosek 8.0.0.41 for sosdemo6 to pass; see:
# https://github.com/JuliaOpt/Mosek.jl/issues/98
mos && push!(sdp_factories, with_optimizer(MathOptInterfaceMosek.MosekOptimizer, QUIET=true))
# Currently, need to create a file param.csdp to avoid printing, see https://github.com/JuliaOpt/CSDP.jl/issues/15
csd && push!(sdp_factories, with_opt(CSDP, printlevel=0))
iscsdp(factory) = occursin("CSDP", string(factory.constructor))
# Need 54000 iterations for sosdemo3 to pass on Linux 64 bits
# With 55000, sosdemo3 passes for every platform except Windows 64 bits on AppVeyor
scs && push!(sdp_factories, SCS.SCSOptimizer(eps=1e-6, max_iters=60000, verbose=0))
isscs(factory) = occursin("SCS", string(factory.constructor))


text:
module GalaxyWarp

    using ColorSchemes
    using JLD2
    #using NPZ
    using CairoMakie
    using Colors
    using Images: imrotate
    using Statistics: mean,median, std
    using Distributions: Normal,Uniform,TriangularDist
    using BlackBoxOptim
    using Optim
    using FITSIO
    using Cosmology
    using Unitful, UnitfulAstro, UnitfulAngles
    import PhysicalConstants.CODATA2018: c_0
    using LinearAlgebra: norm,cross,dot
    #using NumericIO
    using QuadGK
    using Parameters
    using ForwardDiff
    using ProgressLogging
    #using BenchmarkTools
    using ImageFiltering

    include("clouds.jl")
    include("modelling.jl")
    include("tools.jl")
    include("simulate.jl")
    include("vizualization.jl")

    #export load_cube
end


text:
# Sample script for plotting the phase space density near a given spatial location.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot

file = "bulk.0000001.vlsv"

meta = load(file)

coordinates = [0.0, 0.0, 0.0]

plot_vdf(meta, coordinates; verbose=true)

# Show the spatial distribution of cells with saved VDF
cellswithVDF = getcellwithvdf(meta)
cellID = meta["CellID"]
cellID = reshape(cellID, meta.ncells[1], meta.ncells[2])
f_saved_index_ = [findfirst(==(cid), cellID) for cid in cellswithVDF]

x = LinRange(meta.coordmin[1], meta.coordmax[1], meta.ncells[1]) ./ Vlasiator.Re
y = LinRange(meta.coordmin[2], meta.coordmax[2], meta.ncells[2]) ./ Vlasiator.Re

xcell = zeros(size(cellswithVDF))
ycell = similar(xcell)

for i in eachindex(f_saved_index_)
   xcell[i] = x[f_saved_index_[i][1]]
   ycell[i] = y[f_saved_index_[i][2]]
end

figure()
pcolormesh(meta, "vg_pressure", colorscale=Linear)
ax = plt.gca()
ax.scatter(xcell, ycell, marker="+", color="w")


text:
function bisection(f::Function, xl::Number, xu::Number, ea::Number, maxIter::Int64)
  # Estimate roots of functions using the bisection method
  xr = NaN;
  fxl = f(xl);

  for i = 1:maxIter
    xrlast = xr;
    xr = 0.5(xl+xu);
    fxr = f(xr);
    relativeError = abs((xr - xrlast) / xr) * 100;

    if fxr * fxl == 0
      break;
    elseif fxr *fxl > 0
      xl = xr;
      fxl = fxr;
    else
      xu = xr;
    end

    if relativeError < ea
      break;
    end
  end

  return xr;
end

max_stars_repo_path: Julia/Root Estimation/bisection.jl

max_stars_repo_name: alexjohnj/numerical-methods

max_issues_repo_path: Julia/Root Estimation/bisection.jl

max_issues_repo_name: alexjohnj/numerical-methods

max_forks_repo_path: Julia/Root Estimation/bisection.jl

max_forks_repo_name: alexjohnj/numerical-methods


text:
using AstroBase
using Test

@testset "States" begin
    r = [
    6068279.27,
    -1692843.94,
    -2516619.18,
    ]/1000

    v = [
    -660.415582,
    5495.938726,
    -5303.093233,
    ]/1000

    rv = [r; v]

    eph = VSOP87()

    ep = UTCEpoch("2016-05-30T12:00:00.000")
    tdb = TDBEpoch(ep)

    iss = State(ep, r, v)
    rot = Rotation(icrf, iau_earth, ep)
    pos_mars, vel_mars = state(eph, ep, earth, mars)
    r_mars = r .- pos_mars
    v_mars = v .- vel_mars
    mars_rot = State(tdb, r, v, frame=iau_mars, body=mars)

    @testset "Operators" begin
        @test iss == State(ep, r, v, frame=icrf, body=earth)
        @test iss ≈ State(ep, r .+ eps(), v .+ eps(), frame=icrf, body=earth)
    end
    @testset "Change frame" begin
        @test State(iss, eph, frame=iau_earth) ≈ State(ep, rot(r, v)..., frame=iau_earth)
        @test iss ≈ State(State(iss, eph, frame=iau_earth), eph, frame=icrf)
    end
    @testset "Change time scale" begin
        @test State(iss, eph, scale=TDB) == State(tdb, r, v)
    end
    @testset "Change central body" begin
        @test State(iss, eph, body=mars) ==
            State(ep, r_mars, v_mars, frame=icrf, body=mars)
        @test State(State(mars_rot, eph, body=earth), eph, body=mars) ≈ mars_rot
    end
    @testset "Change scale and body" begin
        @test State(iss, eph, scale=TDB, body=mars) ==
            State(tdb, r_mars, v_mars, frame=icrf, body=mars)
    end
    @testset "Change frame and body" begin
        @test State(iss, eph, frame=iau_earth, body=mars) ≈
            State(ep, rot(r_mars, v_mars)..., frame=iau_earth, body=mars)
        @test State(State(mars_rot, eph, frame=itrf, body=earth), eph, frame=iau_mars,
            body=mars) ≈ mars_rot
    end
    @testset "Change frame and scale" begin
        @test State(iss, eph, frame=iau_earth, scale=TDB) ≈
            State(tdb, rot(r, v)..., frame=iau_earth)
    end
    @testset "Change scale, frame, and body" begin
        @test State(State(mars_rot, eph, frame=itrf, scale=UT1, body=earth), eph,
            frame=iau_mars, scale=TDB, body=mars) ≈ mars_rot
    end
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DiscretizeRelax/utilities/access_functions.jl
# Defines methods to access attributes via DynamicBoundsBase interface.
#############################################################################

DBB.supports(::DiscretizeRelax, ::DBB.IntegratorName) = true
DBB.supports(::DiscretizeRelax, ::DBB.Subgradient{T}) where {T <: AbstractBoundLoc} = true
DBB.supports(::DiscretizeRelax, ::DBB.Bound{T}) where {T <: AbstractBoundLoc} = true
DBB.supports(::DiscretizeRelax, ::DBB.Relaxation{T}) where {T <: AbstractBoundLoc} = true
DBB.supports(::DiscretizeRelax, ::DBB.IsNumeric) = true
DBB.supports(::DiscretizeRelax, ::DBB.IsSolutionSet) = true
DBB.supports(::DiscretizeRelax, ::DBB.TerminationStatus) = true
DBB.supports(::DiscretizeRelax, ::DBB.Value) = true
DBB.supports(::DiscretizeRelax, ::DBB.ParameterBound{Lower}) = true
DBB.supports(::DiscretizeRelax, ::DBB.ParameterBound{Upper}) = true
DBB.supports(::DiscretizeRelax, ::DBB.ParameterValue) = true
DBB.supports(::DiscretizeRelax, ::DBB.SupportSet) = true
DBB.supports(::DiscretizeRelax, ::DBB.ParameterNumber) = true
DBB.supports(::DiscretizeRelax, ::DBB.StateNumber) = true
DBB.supports(::DiscretizeRelax, ::DBB.SupportNumber) = true
DBB.supports(t::DiscretizeRelax, ::DBB.LocalSensitivityOn) = t.calculate_local_sensitivity

function get_val_loc(t::DiscretizeRelax, index::Int64, time::Float64)
    (index <= 0 && time == -Inf) && error("Must set either index or time.")
    if index > 0
        return t.relax_t_dict_indx[index]
    end
    t.relax_t_dict_flt[time]
end

function get_val_loc_local(t::DiscretizeRelax, index::Int64, time::Float64)
    (index <= 0 && time == -Inf) && error("Must set either index or time.")
    if index > 0
        return t.local_t_dict_indx[index]
    end
    t.local_t_dict_flt[time]
end

DBB.get(t::DiscretizeRelax, v::DBB.IntegratorName) = "Discretize & Relax Integrator" # TO DO... FIX ME
DBB.get(t::DiscretizeRelax, v::DBB.IsNumeric) = false # TO DO... FIX ME
DBB.get(t::DiscretizeRelax, v::DBB.IsSolutionSet) = true
DBB.get(t::DiscretizeRelax, v::DBB.TerminationStatus) = t.error_code
DBB.get(t::DiscretizeRelax, v::DBB.SupportSet) = DBB.SupportSet(t.tsupports)
DBB.get(t::DiscretizeRelax, v::DBB.ParameterNumber) = t.np
DBB.get(t::DiscretizeRelax, v::DBB.StateNumber) = t.nx
DBB.get(t::DiscretizeRelax, v::DBB.SupportNumber) = length(t.tsupports)
DBB.get(t::DiscretizeRelax, v::DBB.LocalSensitivityOn) = t.calculate_local_sensitivity

DBB.getall(t::DiscretizeRelax, v::DBB.ParameterBound{Lower}) = t.pL
DBB.getall(t::DiscretizeRelax, v::DBB.ParameterBound{Upper}) = t.pU

function DBB.set!(t::DiscretizeRelax, v::DBB.LocalSensitivityOn, q::Bool)
    t.calculate_local_sensitivity = q
    return
end

function DBB.set!(t::DiscretizeRelax, v::DBB.SupportSet)
    t.time = v.s
    t.tsupports = v.s
    nothing
end

## Setting problem attributes
function DBB.setall!(t::DiscretizeRelax, v::ParameterBound{Lower}, value::Vector{Float64})
    t.new_decision_box = true
    @inbounds for i in 1:t.np
        t.pL[i] = value[i]
    end
    return
end

function DBB.setall!(t::DiscretizeRelax, v::ParameterBound{Upper}, value::Vector{Float64})
    t.new_decision_box = true
    @inbounds for i in 1:t.np
        t.pU[i] = value[i]
    end
    return
end

function DBB.setall!(t::DiscretizeRelax, v::ParameterValue, value::Vector{Float64})
    t.new_decision_pnt = true
    @inbounds for i in 1:t.np
        t.p[i] = value[i]
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax, v::DBB.ParameterValue)
    out .= t.p[1:t.np]
    return
end

function DBB.getall!(out::Array{Float64,2}, t::DiscretizeRelax, v::DBB.Value)
    copyto!(out, t.local_problem_storage.pode_x)
    return
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax, v::DBB.Value)
    val_loc = get_val_loc_local(t, v.index, v.time)
    out .= t.local_problem_storage.pode_x[:, val_loc]
    return
end

## Inplace integrator acccess functions
function DBB.getall!(out::Vector{Matrix{Float64}}, t::DiscretizeRelax{X,T}, ::DBB.Subgradient{Lower}) where {X, T <: AbstractInterval}
    for i in 1:t.np
        fill!(out[i], 0.0)
    end
    return
end
function DBB.getall!(out::Vector{Matrix{Float64}}, t::DiscretizeRelax{X,T}, ::DBB.Subgradient{Lower}) where {X, T <: MC}
    for i = 1:length(t.storage)
        for j = 1:t.np
            for k = 1:t.nx
                out[j][k,i]= t.storage[i][k].cv_grad[j]
            end
        end
    end
    return
end

function DBB.getall!(out::Vector{Matrix{Float64}}, t::DiscretizeRelax{X,T}, ::DBB.Subgradient{Upper}) where {X, T <: AbstractInterval}
    for i in 1:t.np
        fill!(out[i], 0.0)
    end
    return
end

function DBB.getall!(out::Vector{Matrix{Float64}}, t::DiscretizeRelax{X,T}, ::DBB.Subgradient{Upper}) where {X, T <: MC}
    for i = 1:length(t.storage)
        for j = 1:t.np
            for k = 1:t.nx
                out[j][k,i]= t.storage[i][k].cc_grad[j]
            end
        end
    end
    return
end

function DBB.get(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Subgradient{Lower}) where {X, T <: AbstractInterval}
    fill!(out, 0.0)
    return
end

function DBB.get(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Subgradient{Upper}) where {X, T <: AbstractInterval}
    fill!(out, 0.0)
    return
end

function DBB.get(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Subgradient{Lower}) where {X, T <: MC}
    val_loc = get_val_loc(t, v.index, v.time)
    for i = 1:t.np
        for j = 1:t.nx
            out[j,i] = t.relax_cv_grad[j,val_loc][i]
        end
    end
    return
end

function DBB.get(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Subgradient{Upper}) where {X, T <: MC}
    val_loc = get_val_loc(t, v.index, v.time)
    for i = 1:t.np
        for j = 1:t.nx
            out[j,i] = t.relax_cc_grad[j,val_loc][i]
        end
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Bound{Lower}) where {X, T <: AbstractInterval}
    @inbounds for j in eachindex(out)
        out[j] = t.storage[j][1].lo
    end
    return
end

function DBB.getall!(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Bound{Lower}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].lo
        end
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Bound{Upper}) where {X, T <: AbstractInterval}
    @inbounds for j in eachindex(out)
        out[j] = t.storage[j][1].hi
    end
    return
end
function DBB.getall!(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Bound{Upper}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].hi
        end
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Bound{Lower}) where {X, T <: MC}
    @inbounds for j in eachindex(out)
        out[j] = t.storage[j][1].Intv.lo
    end
    return
end

function DBB.getall!(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Bound{Lower}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].Intv.lo
        end
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Bound{Upper}) where {X, T <: MC}
    @inbounds for j in eachindex(out)
        out[j] = t.storage[j][1].Intv.hi
    end
    return
end
function DBB.getall!(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Bound{Upper}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].Intv.hi
        end
    end
    return
end

function DBB.getall!(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Lower}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].lo
        end
    end
    return
end

function DBB.getall!(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Lower}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].cv
        end
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Lower}) where {X, T <: AbstractInterval}
    @inbounds for j in eachindex(out)
        out[j] = t.storage[j][1].lo
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Lower}) where {X, T <: MC}
    @inbounds for j in eachindex(out)
        out[j] = t.storage[j][1].cv
    end
    return
end

function DBB.getall!(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Upper}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].hi
        end
    end
    return
end

function DBB.getall!(out::Matrix{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Upper}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].cc
        end
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Upper}) where {X, T <: AbstractInterval}
    @inbounds for j in eachindex(out)
        out[j] = t.storage[j][1].hi
    end
    return
end

function DBB.getall!(out::Vector{Float64}, t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Upper}) where {X, T <: MC}
    @inbounds for j in eachindex(out)
        out[j] = t.storage[j][1].cc
    end
    return
end

## Out of place integrator acccess functions
function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Subgradient{Lower}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = Matrix{Float64}[]
    for i = 1:t.np
        push!(out, zeros(Float64, dims[2], dims[1]))
    end
    out
end
function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Subgradient{Lower}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = Matrix{Float64}[]
    for i = 1:t.np
        push!(out, zeros(Float64, dims[2], dims[1]))
        for j = 1:dims[1]
            for k = 1:dims[2]
                out[i][k, j] = t.storage[j][k].cv_grad[i]
            end
        end
    end
    out
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Subgradient{Upper}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = Matrix{Float64}[]
    for i = 1:t.np
        push!(out, zeros(Float64, dims[2], dims[1]))
    end
    out
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Subgradient{Upper}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = Matrix{Float64}[]
    for i = 1:t.np
        push!(out, zeros(Float64, dims[2], dims[1]))
        for j = 1:dims[1]
            for k = 1:dims[2]
                out[i][k, j] = t.storage[j][k].cc_grad[i]
            end
        end
    end
    out
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Bound{Lower}) where {X, T <: AbstractInterval}
    val_loc = get_val_loc(t, v.index, v.time)
    @__dot__ out = lo(t.storage[val_loc])
    return
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Bound{Upper}) where {X, T <: AbstractInterval}
    val_loc = get_val_loc(t, v.index, v.time)
    @__dot__ out = hi(t.storage[val_loc])
    return
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Bound{Lower}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = zeros(Float64, dims[2], dims[1])
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].lo
        end
    end
    out
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Bound{Upper}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = zeros(Float64, dims[2], dims[1])
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].hi
        end
    end
    out
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Bound{Lower}) where {X, T <: MC}
    val_loc = get_val_loc_local(t, v.index, v.time)
    @__dot__ out = lo(t.storage[val_loc])
    return
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Bound{Upper}) where {X, T <: MC}
    val_loc = get_val_loc_local(t, v.index, v.time)
    @__dot__ out = hi(t.storage[val_loc])
    return
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Bound{Lower}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = zeros(Float64, dims[2], dims[1])
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].Intv.lo
        end
    end
    out
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Bound{Upper}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = zeros(Float64, dims[2], dims[1])
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].Intv.hi
        end
    end
    out
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Relaxation{Lower}) where {X, T <: MC}
    val_loc = get_val_loc(t, v.index, v.time)
    @__dot__ out = cv(t.storage[val_loc])
    return
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Relaxation{Upper}) where {X, T <: MC}
    val_loc = get_val_loc(t, v.index, v.time)
    @__dot__ out = cc(t.storage[val_loc])
    return
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Relaxation{Lower}) where {X, T <: AbstractInterval}
    val_loc = get_val_loc(t, v.index, v.time)
    @__dot__ out = lo(t.storage[val_loc])
    return
end

function DBB.get(out::Vector{Float64}, t::DiscretizeRelax{X,T}, v::DBB.Relaxation{Upper}) where {X, T <: AbstractInterval}
    val_loc = get_val_loc(t, v.index, v.time)
    @__dot__ out = hi(t.storage[val_loc])
    return
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Lower}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = zeros(Float64, dims[2], dims[1])
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].lo
        end
    end
    out
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Lower}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = zeros(Float64, dims[2], dims[1])
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].cv
        end
    end
    out
end

function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Upper}) where {X, T <: AbstractInterval}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = zeros(Float64, dims[2], dims[1])
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].hi
        end
    end
    out
end
function DBB.getall(t::DiscretizeRelax{X,T}, ::DBB.Relaxation{Upper}) where {X, T <: MC}
    dims = size(t.storage)
    (length(dims) == 1) && (dims = (dims[1],1))
    out = zeros(Float64,  dims[2], dims[1])
    for i = 1:dims[1]
        for j = 1:dims[2]
            out[j, i] = t.storage[i][j].cc
        end
    end
    out
end

function DBB.get(t::DiscretizeRelax{X,T}, v::DBB.LocalIntegrator) where {X, T}
    return t.local_problem_storage
end

DBB.get(t::DiscretizeRelax{X,T}, v::DBB.AttachedProblem) where {X, T} = t.prob


text:
include("go_compile.jl")
include("rr.jl")
include("unittests.jl")


text:
using Test
a = 1
b = 2
c = a + b
@test c == 3


text:
using  jInv.Mesh
using  jInv.Utils
using  jInv.LinearSolvers
using  jInv.InverseSolve
using  jInv.Vis
using  EikonalInv
using  MAT
using  FWI
using  ForwardHelmholtz
using  Multigrid
#############################################################################################################
modelDir = "../BenchmarkModels";

dataDir = pwd();
include("../../ex2DFWI/readModelAndGenerateMeshMref.jl");
include("../../ex2DFWI/prepareFWIDataFiles.jl");
include("Drivers/setupFWI.jl");

plotting        = false;
large           = false
iterativeSolver = true;
## WorkersFWI can be enlarged here, but note to change numCores in the solvers section such that #workers and #cores make sense in your system.
## Each worker will use numCores cores for its computations.

workersFWI 		= [workers()[1]];
if plotting
	using  PyPlot
end
#############################################################################################################
dataDir = pwd();
resultsDir = pwd();
########################################################################################################
matfile   = matread("../3Dseg12812864.mat")
m         = matfile["VELc"]*1e-3
m         = reshape(m,(128,128,64));

if large==false
	pad     	 = 10;
	ABLPad 		 = pad + 8;
	newSize 	 = [145,145,70];
	omega   	 = [0.5,0.75,1.25,1.75]*2*pi;
	jumpSrc 	 = 16;
	### FOR A QUICK RUN USE THIS INSTEAD OF THE THREE LINES ABOVE:
	newSize 	 = [65,65,33];
	omega   	 = [0.5]*2*pi;
	jumpSrc 	 = 8;
	maxBatchSize     = 27;
else
	pad     	 = 20;
	ABLPad 		 = pad + 5;
	newSize 	 = [236,236,120];
	omega   	 = [0.5,0.8,1.2,1.8,2.5]*2*pi;
	jumpSrc 	 = 25;
	maxBatchSize = 25;
end

offset  = ceil(Int64,(newSize[1]*(16.0/13.5)));
println("Offset is: ",offset)
domain = [0.0,13.5,0.0,13.5,0.0,4.2];

(m,Minv,mref,boundsHigh,boundsLow) = readModelAndGenerateMeshMref(m,pad,newSize,domain);

useFilesForFields = true;

# ###################################################################################################################
dataFilenamePrefix = string(dataDir,"/DATA_SEG",tuple((Minv.n+1)...));
resultsFilename = string(resultsDir,"/SEG");
#######################################################################################################################

if plotting
	limits = [1.5,4.5];
	figure(1,figsize = (22,10))
	plotModel(m,true,Minv,pad,limits);
	figure(2,figsize = (22,10))# ,figsize = (22,10)
	plotModel(mref,true,Minv,pad,limits);

	m_t,M_t = cutBoundaryLayer(mref,Minv,pad);
	writedlm(string(resultsFilename,tuple((M_t.n+1)...),"_mref.dat"),convert(Array{Float16},m_t));
	m_t,M_t = cutBoundaryLayer(m,Minv,pad);
	writedlm(string(resultsFilename,tuple((M_t.n+1)...),"_mtrue.dat"),convert(Array{Float16},m_t));
	m_t = 0;
	M_t = 0;
end

######################## ITERATIVE SOLVER #############################################
if iterativeSolver==true
	levels      = 3;
	numCores 	= 24;
	BLAS.set_num_threads(numCores);
	maxIter     = 50;
	relativeTol = 1e-5;
	relaxType   = "SPAI";
	relaxParam  = 1.0;
	relaxPre 	= 2;
	relaxPost   = 2;
	cycleType   ='W';
	coarseSolveType = "MUMPS";
	MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.0,0.0);
	shift = 0.15;
	Hparam = HelmholtzParam(Minv,zeros(0),zeros(0),0.0,true,true);
	Ainv = getShiftedLaplacianMultigridSolver(Hparam, MG,shift);
else   ######################## DIRECT SOLVER #################################################
	numCores 	= 16;
	BLAS.set_num_threads(numCores);
	Ainv = getMUMPSsolver([],0,0,2);
	# Ainv = getJuliaSolver();
	# Ainv = getPARsolver([],0,0,6);
end

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

println("omega*maximum(h): ",omega*maximum(Minv.h)*sqrt(maximum(1./(boundsLow.^2))));

# This is a list of workers for FWI. Ideally they should be on different machines.

println("The workers that we allocate for FWI are:");
println(workersFWI)

prepareFWIDataFiles(m,Minv,mref,boundsHigh,boundsLow,dataFilenamePrefix,omega,one(Complex128)*ones(size(omega)), pad,ABLPad,jumpSrc,
					offset,workersFWI,maxBatchSize,Ainv,useFilesForFields);

(Q,P,pMis,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow,resultsFilename) =
   setupFWI(m,dataFilenamePrefix,resultsFilename,plotting,workersFWI,maxBatchSize,Ainv,SSDFun,useFilesForFields);

########################################################################################################
# Setting up the inversion for slowness
########################################################################################################
function dump(mc,Dc,iter,pInv,PMis,resultsFilename)
	fullMc = slowSquaredToVelocity(reshape(Iact*pInv.modelfun(mc)[1] + sback,tuple((pInv.MInv.n+1)...)))[1];
	Temp = splitext(resultsFilename);
	if iter>0
		Temp = string(Temp[1],iter,Temp[2]);
	else
		Temp = resultsFilename;
	end
	if resultsFilename!=""
		writedlm(Temp,convert(Array{Float16},fullMc));
	end
	if plotting
		close(888);
		figure(888);
		plotModel(fullMc,true,false,[],0,[1.5,4.8],splitdir(Temp)[2]);
	end
end

mref 		= velocityToSlow(mref)[1];
t    		= copy(boundsLow);
boundsLow 	= velocityToSlow(boundsHigh)[1];
boundsHigh 	= velocityToSlow(t)[1]; t = 0;
modfun 		= slowToSlowSquared;

########################################################################################################
# Set up Inversion #################################################################################
########################################################################################################

maxStep=0.05*maximum(boundsHigh);

regparams = [1.0,1.0,1.0,1e-5];
regfun(m,mref,M) = wdiffusionRegNodal(m,mref,M,Iact=Iact,C=regparams);
cgit = 9;
alpha = 1e-10;
pcgTol = 1e-1;
maxit = 10;

HesPrec = getSSORCGRegularizationPreconditioner(1.0,1e-5,1000);

pInv = getInverseParam(Minv,modfun,regfun,alpha,mref[:],boundsLow,boundsHigh,
                         maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
						 minUpdate=1e-3, maxIter = maxit,HesPrec=HesPrec);
mc = copy(mref[:]);
mc,Dc,His = freqCont(mc, pInv, pMis,contDiv, 2, resultsFilename,dump,"Joint",1,1,"projGN");
# matwrite("lrSEG-res1.mat",Dict("mc"=>mc,"His"=>His))
mc,Dc,His = freqCont(mc, pInv, pMis,contDiv, 2, resultsFilename,dump,"Joint",2,2,"projGN");
# matwrite("lrSEG-res2.mat",Dict("mc"=>mc,"His"=>His))
mc,Dc,His = freqCont(mc, pInv, pMis,contDiv, 2, resultsFilename,dump,"Joint",3,3,"projGN");
# matwrite("lrSEG-res3.mat",Dict("mc"=>mc,"His"=>His))
##############################################################################################


text:
# Main archetype function
function buildarchetypes(project)
    stockdata = project.dat["master"]

    # Identify different plan shapes, based on vertex count
    L = length.(stockdata.dims)
    vertexcount(L,i) = findall(L.==i)
    ishapes = Dict()
    push!(ishapes,"4vert"=>vertexcount(L,4))
    push!(ishapes,"6vert"=>vertexcount(L,6))
    push!(ishapes,"8vert"=>vertexcount(L,8))

    # Prepare new columns for archetype fields
    L = nrow(stockdata)
    stockdata[!,:width] = fill!(Array{Union{Missing,Float64}}(undef,L),missing)
    stockdata[!,:depth] = fill!(Array{Union{Missing,Float64}}(undef,L),missing)
    stockdata[!,:tol] = fill!(Array{Union{Missing,Float64}}(undef,L),missing)
    stockdata[!,:tor] = fill!(Array{Union{Missing,Float64}}(undef,L),missing)
    stockdata[!,:tp] = fill!(Array{Union{Missing,Float64}}(undef,L),missing)
    stockdata[!,:bc1] = fill!(Array{Union{Missing,String}}(undef,L),missing)

    stockdata = assess4vertplans(stockdata,ishapes["4vert"])
    return project
end

# 4 vertex dwellings
function assess4vertplans(stockdata,ishape)
    elevpair(verts) = [round(mean(verts[[1,3]]),digits=3),round(mean(verts[[2,4]]),digits=3)]
    plans_dims = elevpair.(stockdata.dims[ishape])

    plans = DataFrame(
        :width => maximum.(plans_dims),
        :depth => minimum.(plans_dims)
    )

    adiab_odd(verts) = sum(iseven.(findall(verts.==true)))<sum(isodd.(findall(verts.==true)))
    adiab_even(verts) = sum(iseven.(findall(verts.==true)))>sum(isodd.(findall(verts.==true)))

    i_even = findall(adiab_even.(stockdata.adiabBool[ishape]))
    i_odd = findall(adiab_odd.(stockdata.adiabBool[ishape]))

    dim_i(dims,i) = dims[i]
    plans.width[i_even] .= dim_i.(plans_dims[i_even],1)
    plans.width[i_odd] .= dim_i.(plans_dims[i_odd],2)
    plans.depth[i_even] .= dim_i.(plans_dims[i_even],2)
    plans.depth[i_odd] .= dim_i.(plans_dims[i_odd],1)

    stockdata.width[ishape] .= plans.width
    stockdata.depth[ishape] .= plans.depth
    stockdata.tol[ishape] .= 0.
    stockdata.tor[ishape] .= 0.
    stockdata.tp[ishape] .= 0.

    function findadjacencies(adiab_array)
        i_all = findall(adiab_array)
        adjacencies = missing
        if isempty(i_all)
            adjacencies = "EEEE"
        elseif length(i_all) == 1
            adjacencies = "EEEA"
        elseif length(i_all) == 2
            if diff(i_all) == [2]
                adjacencies = "EEAA"
            else
                adjacencies = "EAEA"
            end
        elseif length(i_all) == 3
            adjacencies = "EAAA"
        end
        return adjacencies
    end
    stockdata.bc1[ishape] = findadjacencies.(stockdata.adiabBool[ishape])

    return stockdata
end

function abovebelowadj(stockdata)
    # Prep array
    stockdata[!,:bc2] = fill!(Array{Union{Missing,String}}(undef,nrow(stockdata)),missing)
    # Apply External (E) roof and Ground (G) adjacencies to all single-UPRN footprints
    # ... THIS OVERLOOKS FLATS ABOVE NON-RES,BUT THESE ARE PICKED UP BELOW
    stockdata.bc2[stockdata.nUPRNs_undersameroof.==1] .= "EG"

    findall_arr(val,arr) = findall(arr.==val)
    i_1to4UPRN = vcat(findall_arr.((stockdata.nUPRNs_undersameroof,),1:4)...)
    i_nonmissing = findall(ismissing.(stockdata.EPC_Dwelling_Type).==false)
    findoccursin_arr(val,arr) = i_nonmissing[findall(occursin.(val,arr))]
    j = findoccursin_arr.(("ground-floor ","mid-floor ","top-floor ","basement "),(lowercase.(stockdata.EPC_Dwelling_Type[i_nonmissing]),))
    stockdata.bc2[intersect(i_1to4UPRN,j[1])] .= "AG"
    stockdata.bc2[intersect(i_1to4UPRN,j[2])] .= "AA"
    stockdata.bc2[intersect(i_1to4UPRN,j[3])] .= "EA"
    stockdata.bc2[intersect(i_1to4UPRN,j[4])] .= "AG"
    return stockdata
end

function gen_archcode(stockdata,tol_geom)
    # Find intersecting non-missing rows (all data present in geom cols)
    ismissing_col(data,col) = findall(ismissing.(data[:,col]).==false)
    i_nonmisisngs = intersect(vcat(
        intersect(ismissing_col.((stockdata,),[:width,:depth,:tol,:tor,:tp]))...),
        intersect(vcat(ismissing_col.((stockdata,),[:RelH2,:RelHMax]))
    ...))
    stockdata_exrp = deepcopy(stockdata[i_nonmisisngs,:])

    # Apply archetyping geometry tolerence (arg, in metres)
    apply_geomtol(arr,tol,factor) = lpad(string(Int(factor*round(arr./tol).*tol)),3,"0")
    w = apply_geomtol.(stockdata_exrp.width,tol_geom,10)
    d = apply_geomtol.(stockdata_exrp.depth,tol_geom,10)
    tol = apply_geomtol.(stockdata_exrp.tol,tol_geom,10)
    tor = apply_geomtol.(stockdata_exrp.tor,tol_geom,10)
    tp = apply_geomtol.(stockdata_exrp.tp,tol_geom,10)
    he = apply_geomtol.(stockdata_exrp.RelH2,tol_geom/2,100)
    ha = apply_geomtol.(stockdata_exrp.RelHMax,tol_geom/2,100)

    # Compile archetypes
    compilearchetype(w,d,tol,tor,tp,he,ha,bc1,bc2) =
        "w$w-d$d-tol$tol-tor$tor-tp$tp-he$he-ha$ha-adj$bc1$bc2"
    archcode = compilearchetype.(w,d,tol,tor,tp,he,ha,
        stockdata_exrp.bc1,stockdata_exrp.bc2)

    return archcode
end


text:
module Tests

using NumericSuffixes
using Base.Test

@suffix I8 Int8
@suffix t16 x -> trunc(Int16, x)

@test 0I8 == Int8(0)
@test 8t16 == trunc(Int16, 8)

for (k, v) in NumericSuffixes.suffixes
    suffix = getfield(NumericSuffixes, k)
    totype = getfield(Base, v)
    @test 5suffix == totype(5)
end

@test_throws ErrorException NumericSuffixes.gentype(:(x + y))
@test_throws ErrorException NumericSuffixes.genfunc(:x, :(x + y))
@test_throws ErrorException NumericSuffixes.genfunc(:x, :((x, y) -> ()))
@test_throws ErrorException NumericSuffixes.genfunc(:x, 0)

end


text:
using NumericalMethodsforEngineers, Plots

old = pwd()
ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  f(x::Float64, y::Vector{Float64}) = (x + y[1])/x
  steps = 3
  h = 0.333

  x = 2.0
  y = [2.0]

  # y(3.0) = 4.2165

  r = Array{Float64,2}[]
  push!(r, euler(f, x, y, steps, h))
  push!(r, modified_euler(f, x, y, steps, h))
  push!(r, mid_point_euler(f, x, y, steps, h))
  push!(r, runga_kutta_4(f, x, y, steps, h))
  println()
  r |> display
  println()

  xs = x + Float64[(i-1)*h for i in 1:steps+1]
  titles = ["Euler", "Modified_Euler", "Mid_Point_Euler", "Runga_Kutta_4"]

  p1 = plot(xs, r[1][:, 2])
  plot!(p1, xs, r[2][:, 2])
  plot!(p1, xs, r[3][:, 2])
  plot!(p1, xs, r[4][:, 2])
  show(p1)
  savefig("Ex7.2.png")
  gui()

end


text:
using Sigma

## 1D Ising Model
## ==============
nbits = 10
bits = Sigma.RandArray([Sigma.flip() for i = 1:nbits])
noisyeq(a,b) = Sigma.flip(ifelse(a == b, 1.0, 0.2))

ising_condition = (&)(map(noisyeq, bits[2:end], bits[1:end-1])...)
rand(bits, ising_condition)


text:
##

using RadiativeTransfer, RadiativeTransfer.vSmartMOM

##

I_modeled_all = zeros(7, 16);
I_modeled_all_I = zeros(7, 16);
I_deltas_all = zeros(7, 16);

Q_modeled_all = zeros(7, 16);
Q_deltas_all = zeros(7, 16);

U_modeled_all = zeros(7, 16);
U_deltas_all = zeros(7, 16);

##

μ = [0.02, 0.06, 0.10, 0.16, 0.20, 0.28, 0.32, 0.40, 0.52, 0.64, 0.72, 0.84, 0.92, 0.96, 0.98, 1.00]
ϕs = collect(0.0:30.0:180.0)

acosd.(μ)

include("test/benchmarks/natraj_trues.jl")
parameters = parameters_from_yaml("test/benchmarks/natraj.yaml");

λ = 360.0
parameters.spec_bands = [1e7/λ (1e7/λ + 1)]
parameters.vza = acosd.(μ)
parameters.sza = acosd.(0.2)

τ = 0.5

for ϕ_i in 1:length(ϕs)
    parameters.vaz = repeat([ϕs[ϕ_i]], 16)
    model = model_from_parameters(parameters);
    model.τ_rayl[1] .= τ

    R = vSmartMOM.rt_run(model, i_band=1)
    @show size(R)

    I_modeled_all[ϕ_i,:] = R[:,1,1]
    Q_modeled_all[ϕ_i,:] = R[:,2,1]
    U_modeled_all[ϕ_i,:] = R[:,3,1]

    I_deltas_all[ϕ_i,:] = abs.(I_trues[:,ϕ_i] - I_modeled_all[ϕ_i,:]) ./ I_trues[:,ϕ_i]
    Q_deltas_all[ϕ_i,:] = abs.(Q_trues[:,ϕ_i] - Q_modeled_all[ϕ_i,:]) ./ Q_trues[:,ϕ_i]
    U_deltas_all[ϕ_i,:] = abs.(U_trues[:,ϕ_i] - U_modeled_all[ϕ_i,:]) ./ U_trues[:,ϕ_i]
end

parameters = parameters_from_yaml("test/benchmarks/natraj_I.yaml");

λ = 360.0
parameters.spec_bands = [1e7/λ (1e7/λ + 1)]
parameters.vza = acosd.(μ)
parameters.sza = acosd.(0.2)

τ = 0.5

for ϕ_i in 1:length(ϕs)
    parameters.vaz = repeat([ϕs[ϕ_i]], 16)
    model = model_from_parameters(parameters);
    model.τ_rayl[1] .= τ

    R = vSmartMOM.rt_run(model, i_band=1)
    @show size(R)

    I_modeled_all_I[ϕ_i,:] = R[:,1,1]

end

@show maximum(I_deltas_all)
println(maximum(Q_deltas_all[findall(i -> i >= 0.01, Q_modeled_all)]))
@show maximum(filter(!isnan, U_deltas_all[findall(i -> i >= 0.01, U_modeled_all)]))


text:
module Proj

using PROJ_jll
using CEnum
using CoordinateTransformations
using NetworkOptions: ca_roots

export PROJ_jll
export PJ_DIRECTION, PJ_FWD, PJ_IDENT, PJ_INV

"""
    struct Coord <: AbstractVector{Float64}

    Coord(x, y, [z = 0], [t = Inf])
    Coord(v::AbstractVector{<:Real})
    Coord(v)

Represent a coordinate with 4 Float64 values. Struct for interoperability with the PROJ C
API, specifically anything that expects a
[`PJ_COORD`](https://proj.org/development/reference/datatypes.html#c.PJ_COORD).

The constructors accept either 2 to 4 numbers, or 1 collection of 2 to 4 numbers, with the
first two required numbers representing x and y coordinates. Number 3 is the elevation
(defaults to 0), and number 4 is the time in years (defaults to Inf).

The four fields are called x, y, z and t, but note that depending on the [axis
ordering](https://proj.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent) the
field x may represent y and vice versa. You can use getindex `coord[2]` to avoid confusion.

The user generally does not need to create a `Coord` directly, since the input coordinates
are automatically converted to `Coord` by ccall. Therefore the constuctors need to be as
flexible as the call methods on the `Transformation`.

# Examples
```julia
julia> Proj.Coord(1.0, 2.0)
4-element Proj.Coord:
  1.0
  2.0
  0.0
 Inf

julia> Proj.Coord((1.0, 2.0, 3.0))
4-element Proj.Coord:
  1.0
  2.0
  3.0
 Inf
```
"""
struct Coord <: AbstractVector{Float64}
    x::Float64
    y::Float64
    z::Float64
    t::Float64
    Coord(x, y, z = 0.0, t = Inf) = new(x, y, z, t)
end

# this shields a StackOverflow from the splatting constructor
Coord(::Real) = error("Proj.Coord takes 2 to 4 numbers, one given")
Coord(v) = Coord(v...)

function Coord(v::AbstractVector{<:Real})
    n = length(v)
    return if n == 2
        Coord(v[begin], v[begin+1])
    elseif n == 3
        Coord(v[begin], v[begin+1], v[begin+2])
    elseif n == 4
        Coord(v[begin], v[begin+1], v[begin+2], v[begin+3])
    else
        error("Proj.Coord takes 2 to 4 numbers")
    end
end

Base.convert(::Type{Coord}, x) = Coord(x)
Base.convert(::Type{Coord}, x::Coord) = x

Base.length(::Coord) = 4
Base.size(::Coord) = (4,)
Base.getindex(coord::Coord, i::Int) = getfield(coord, i)
Base.IndexStyle(::Type{Coord}) = IndexLinear()
Base.eltype(::Coord) = Float64

# type aliases
const NTuple234 = Union{NTuple{2,Float64},NTuple{3,Float64},NTuple{4,Float64}}
const PROJ_COMPUTE_VERSION = VersionNumber
const GEODESIC_VERSION_NUM = VersionNumber

include("libproj.jl")
include("coord.jl")
include("error.jl")

"""
    unsafe_loadstringlist(ptr::Ptr{Cstring})

Load a null-terminated list of strings.

It takes a `PROJ_STRING_LIST`, which is a `Ptr{Cstring}`, and returns a `Vector{String}`.
"""
function unsafe_loadstringlist(ptr::Ptr{Cstring})
    strings = String[]
    (ptr == C_NULL) && return strings
    i = 1
    cstring = unsafe_load(ptr, i)
    while cstring != C_NULL
        push!(strings, unsafe_string(cstring))
        i += 1
        cstring = unsafe_load(ptr, i)
    end
    proj_string_list_destroy(ptr)
    return strings
end

"Prevent an error converting a null pointer to a string, returns `nothing` instead"
aftercare(x::Cstring) = x == C_NULL ? nothing : unsafe_string(x)
aftercare(x::Ptr{Cstring}) = unsafe_loadstringlist(x)

const PROJ_LIB = Ref{String}()

"Module initialization function"
function __init__()
    # register custom error handler
    funcptr = @cfunction(log_func, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Cstring))
    proj_log_func(C_NULL, funcptr)

    # set path to CA certificates
    ca_path = ca_roots()
    if ca_path !== nothing
        proj_context_set_ca_bundle_path(ca_path)
    end

    # point to the location of the provided shared resources
    PROJ_LIB[] = joinpath(PROJ_jll.artifact_dir, "share", "proj")
    proj_context_set_search_paths(1, [PROJ_LIB[]])
end

end # module


text:
export
    set_func_callback,
    set_grad_callback,
    set_hess_callback,
    # set_newpt_callback,
    # set_ms_process_callback,
    set_mip_node_callback #,
    # set_ms_initpt_callback,
    # set_puts_callback,

const callback_params = (Cint, Cint, Cint, Cint, Cint, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Nothing})

"""
Set the callback function that evaluates `obj` and `c` at `x`.

It may also evaluate `objGrad` and `jac` if EVALFC and EVALGA are
combined into a single call.
"""
function set_func_callback(kp::KnitroProblem, f::Function)
    cb = quote
        @cfunction($f, Cint, (Cint, Cint, Cint, Cint, Cint, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Nothing}))
    end
    return_code = @ktr_ccall(set_func_callback, Int32, (Ptr{Nothing},
                                                        Ptr{Nothing}), kp.env, eval(cb))
    if return_code != 0
        error("KNITRO: Error setting function callback")
    end
end

# /**
#  *  Do not modify "hessian" or "hessVector".
#  */
# int  KNITRO_API KTR_set_grad_callback (KTR_context_ptr       kc,
#                                        KTR_callback * const  fnPtr);

"""
Set the callback function that evaluates `objGrad` and `jac` at `x`.

It may do nothing if EVALFC and EVALGA are combined into a single call.
"""
function set_grad_callback(kp::KnitroProblem, f::Function)
    cb = quote
        @cfunction($f, Cint, (Cint, Cint, Cint, Cint, Cint, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Nothing}))
    end
    return_code = @ktr_ccall(set_grad_callback, Int32,(Ptr{Nothing},
                                                       Ptr{Nothing}), kp.env, eval(cb))
    if return_code != 0
        error("KNITRO: Error setting gradient callback")
    end
end

"""
Set the callback function that evaluates second derivatives at (x, lambda).

If `evalRequestCode` equals KTR_RC_EVALH, then the function must
return nonzeroes in `hess`.  If it equals KTR_RC_EVALHV, then the
function multiplies second derivatives by `hessVector` and returns the
product in `hessVector`.
"""
function set_hess_callback(kp::KnitroProblem, f::Function)
    cb = quote
        @cfunction($f, Cint, (Cint, Cint, Cint, Cint, Cint, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Nothing}))
    end
    return_code = @ktr_ccall(set_hess_callback, Int32, (Ptr{Nothing},
                                                        Ptr{Nothing}), kp.env, eval(cb))
    if return_code != 0
        error("KNITRO: Error setting hessian callback")
    end
end

"""
Set the callback function that is invoked after KNITRO computes a
new estimate of the solution point (i.e., after every major iteration).
"""
function set_newpt_callback(kp::KnitroProblem, f::Function)
    cb = quote
        @cfunction($f, Cint, (Ptr{Nothing}, Cint, Cint, Cint, Ptr{Cdouble},
                   Ptr{Cdouble}, Cdouble, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Nothing}))
    end
    return_code = @ktr_ccall(set_newpt_callback, Int32, (Ptr{Nothing},
                                                         Ptr{Nothing}), kp.env, eval(cb))
    if return_code != 0
        error("KNITRO: Error setting newpoint callback")
    end
end

"""
This callback function is for multistart (MS) problems only.

Set the callback function that is invoked after KNITRO finishes
processing a multistart solve.  The function should not modify any
KNITRO arguments.  Arguments `x` and `lambda` contain the solution from
the last solve. Arguments `obj` and `c` contain objective and constraint
values at `x`.  First and second derivative arguments are not currently
defined and should not be examined.
"""
function set_ms_process_callback(kp::KnitroProblem, f::Function)
    cb = quote
        @cfunction($f, Cint, (Cint, Cint, Cint, Cint, Cint, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Nothing}))
    end
    return_code = @ktr_ccall(set_ms_process_callback, Int32, (Ptr{Nothing},
                             Ptr{Nothing}), kp.env, eval(cb))
    if return_code != 0
        error("KNITRO: Error setting multi-start callback")
    end
end

"""
This callback function is for mixed integer (MIP) problems only.

Set the callback function that is invoked after KNITRO finishes
processing a node on the branch-and-bound tree.
"""
function set_mip_node_callback(kp::KnitroProblem, f::Function)
    cb = quote
        @cfunction($f, Cint, (Cint, Cint, Cint, Cint, Cint, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Nothing}))
    end
    return_code = @ktr_ccall(set_mip_node_callback, Int32, (Ptr{Nothing},
                                                            Ptr{Nothing}), kp.env, eval(cb))
    if return_code != 0
        error("KNITRO: Error setting MIP node callback")
    end
end

"""
This callback function is for multistart (MS) problems only.

Set the callback that allows applications to specify an initial point
before each local solve in the multistart procedure.
"""
function set_ms_initpt_callback(kp::KnitroProblem, f::Function)
    cb = quote
        @cfunction($f, Cint, (Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble},
                   Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Nothing}))
    end
    return_code = @ktr_ccall(set_ms_initpt_callback, Int32, (Ptr{Nothing},
                                                             Ptr{Nothing}), kp.env, eval(cb))
    if return_code < 0
        error("KNITRO: Error setting multi-start initial-point callback")
    end
end

"""
Set a put-string callback function to handle output generated by the KNITRO solver.

The callback function should return the number of characters that were printed.

By default KNITRO prints to stdout or a file named `knitro.log`,
as determined by KTR_PARAM_OUTMODE.
"""
function set_puts_callback(kp::KnitroProblem, f::Function)
    cb = quote
        @cfunction($f, Cint, (Ptr{Cchar}, Ptr{Nothing}))
    end
    return_code = @ktr_ccall(set_puts_callback, Int32, (Ptr{Nothing},
                             Ptr{Nothing}), kp.env, eval(cb))
    if return_code < 0
        error("KNITRO: Error setting put-string callback")
    end
end


text:
## Julia program for Hierarchical Modeling of means and variances
## author: weiya <szcfweiya@gmail.com>
## date: 27 August, 2018

using Distributions
using SpecialFunctions
using StatsBase
using DelimitedFiles

function higibbs(Y, T, mu0 = 50.0, gamma20 = 25.0, nu0 = 1.0, sigma20 = 100.0, eta0 = 1.0, tau20 = 100.0, a = 1.0, b = 1/100.0, alpha = 1.0, NUMAX = 40)
    m = size(unique(Y[:,1]), 1)
    # starting value
    ybar = ones(m)
    sv = ones(m)
    n = ones(m)
    for j = 1:m
        yj = Y[ [Y[i,1] == j for i = 1:end], 2]
        ybar[j] = mean(yj)
        sv[j] = var(yj)
        n[j] = size(yj, 1)
    end
    theta = ybar
    sigma2 = copy(sv)
#    sigma20 = 1 / mean(sigma2)
#    nu0 = 2 * mean(sigma2)^2 / var(sigma2)
    mu = mean(theta)
    tau2 = var(theta)

    THETA = ones(T, m)
    SIGMA2 = ones(T, m)
    # mu tau2 sigma20 nu0
    MTSN = ones(T, 4)

    for t = 1:T
        # sample mu
        varmu = 1 / (m / tau2 + 1 / gamma20)
        meanmu = varmu * (m * mean(theta) / tau2 + mu0 / gamma20)
        rnorm = Normal(meanmu, sqrt(varmu))
        mu = rand(rnorm, 1)[1]

        # sample 1/tau2
        shapetau = (eta0 + m) / 2
        ratetau = ( eta0 * tau20 + sum((theta .- mu).^2) ) / 2
        rgamma = Gamma(shapetau, 1/ratetau)
        tau2 = 1 / rand(rgamma, 1)[1]

        # sample theta
        for j = 1:m
            vartheta = 1 / (n[j] / sigma2[j] + 1 / tau2)
            meantheta = vartheta * ( n[j] * mean(Y[ [Y[i,1] == j for i = 1:end], 2]) / sigma2[j] + mu / tau2)
            rnorm = Normal(meantheta, sqrt(vartheta))
            theta[j] = rand(rnorm, 1)[1]
        end
        THETA[t, :] .= theta

        # sample sigma2
        for j = 1:m
            shapesig = (nu0 + n[j])/2
            yj = Y[ [Y[i,1] == j for i = 1:end], 2]
            ratesig = ( nu0*sigma20 + sum( (yj .- theta[j]).^2 ) )/2
            rgamma = Gamma(shapesig, 1/ratesig)
            sigma2[j] = 1 / rand(rgamma, 1)[1]
        end
        SIGMA2[t, :] .= sigma2

        # sample sigma20
        shapesig = a + 0.5 * m * nu0
        ratesig = b + 0.5 * nu0 * sum(1 ./ sigma2)
        rgamma = Gamma(shapesig, 1/ratesig)
        sigma20 = rand(rgamma, 1)[1]

        # sample nu0
        x = 1:NUMAX
        lpnu0 = ones(NUMAX)
        lpnu0 .= m * ( .5 * x .* log.(sigma20 * x / 2) .- lgamma.(x/2) ) .+ (x / 2 .+ 1) * sum(log.(1 ./ sigma2)) .- x .* (alpha + .5 * sigma20 * sum(1 ./ sigma2))
        #println(lpnu0)
        nu0 = sample(x, pweights(exp.(lpnu0 .- maximum(lpnu0))))
#        println(pweights(exp.(lpnu0 .- maximum(lpnu0))))

        # store results
        MTSN[t, :] .= [mu, tau2, sigma20, nu0]
    end
    return THETA, SIGMA2, MTSN, sv, n
end

# run
Y = readdlm("math-score-Y.csv")
THETA, SIGMA2, MTSN, sv, n = higibbs(Y, 5000)

using PyPlot
# histogram
plt[:subplot](221)
plt[:hist](MTSN[:,2])
ylabel(L"$\mu$")
plt[:subplot](222)
plt[:hist](MTSN[:,2])
ylabel(L"$\tau^2$")
plt[:subplot](223)
plt[:hist](MTSN[:,3])
ylabel(L"$\nu_0$")
plt[:subplot](224)
plt[:hist](MTSN[:,4])
ylabel(L"$\sigma_0^2$")
plt[:tight_layout]()
show()

# shrinkage
f, (ax1, ax2) = plt[:subplots](1, 2)
ax1[:scatter](sv, SIGMA2[end,:])
ax1[:plot](sv, sv)
ax1[:set_xlabel](L"$s^2$")
ax1[:set_ylabel](L"$\hat \sigma^2$")
ax2[:scatter](n, sv-SIGMA2[end,:])
ax2[:plot](n, zeros(size(n, 1)))
ax2[:set_xlabel]("sample size")
ax2[:set_ylabel](L"$s^2-\hat \sigma^2$")
plt[:tight_layout]()
show()

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
using Test
using OnlineStats: Series, fit!
using LyapunovExponents: VecMean, VecVariance

@time @testset "OnlineStats" begin
    max_dim = 10
    len = 5
    rng = MersenneTwister(0)
    @testset "VecMean dim=$dim" for dim in 1:max_dim
        stat = VecMean(dim)
        os = Series(stat)
        xs = randn(rng, dim, len)
        for i in 1:len
            fit!(os, xs[:, i])
            @test mean(stat) ≈ mean(xs[:, 1:i], 2)[:, 1]
        end
    end
    @testset "VecVariance dim=$dim" for dim in 1:max_dim
        stat = VecVariance(dim)
        os = Series(stat)
        xs = randn(rng, dim, len)
        for i in 1:len
            fit!(os, xs[:, i])
            @test mean(stat) ≈ mean(xs[:, 1:i], 2)[:, 1]
            @test var(stat) ≈ var(xs[:, 1:i], 2)[:, 1] nans=true
        end
    end
end


text:
function logisticsoftmax(x::AbstractVector{<:Real})
    σs = logistic.(x)
    return σs / sum(σs)
end

struct LogisticSoftMaxLink{Tθ<:AbstractVector} <: AbstractLink
    logθ::Tθ
end

LogisticSoftMaxLink(nclass::Integer) = LogisticSoftMaxLink(zeros(nclass))

function _get_const(l::BijectiveSimplexLink{<:LogisticSoftMaxLink})
    return exp(last(l.link.logθ)) * logistic(0)
end

_sum_θ(l::LogisticSoftMaxLink) = exp(logsumexp(l.logθ))

function _sum_θ(l::BijectiveSimplexLink{<:LogisticSoftMaxLink})
    return _get_const(l) + exp(logsumexp(l.link.logθ[1:(end - 1)]))
end

function _scale_σf(l::LogisticSoftMaxLink, f::AbstractVector{<:Real})
    return exp.(l.logθ) .* logistic.(f)
end

function _scale_σf(
    l::BijectiveSimplexLink{<:LogisticSoftMaxLink}, f::AbstractVector{<:Real}
)
    return exp.(l.link.logθ[1:(end - 1)]) .* logistic.(f)
end

function (l::LogisticSoftMaxLink)(f::AbstractVector{<:Real})
    σs = exp.(l.logθ) .* logistic.(f)
    return σs ./ sum(σs)
end

# Augmentations are possible for both options
const BijectiveLogisticSoftMaxLikelihood = CategoricalLikelihood{
    <:BijectiveSimplexLink{<:LogisticSoftMaxLink}
}
const LogisticSoftMaxLikelihood = CategoricalLikelihood{<:LogisticSoftMaxLink}
const LogisticSoftMaxLikelihoods = Union{
    BijectiveLogisticSoftMaxLikelihood,LogisticSoftMaxLikelihood
}

nlatent(l::BijectiveLogisticSoftMaxLikelihood) = length(l.invlink.link.logθ) - 1
nlatent(l::LogisticSoftMaxLikelihood) = length(l.invlink.logθ)

aux_field(::LogisticSoftMaxLikelihoods, Ω::NamedTuple) = values(Ω)
aux_field(::LogisticSoftMaxLikelihoods, Ω::TupleVector) = zip(Ω.ω, Ω.n)

function init_aux_variables(rng::AbstractRNG, l::LogisticSoftMaxLikelihoods, ndata::Int)
    return TupleVector(;
        ω=nestedview(rand(rng, PolyaGamma(1, 0), nlatent(l), ndata)),
        n=nestedview(rand(rng, Poisson(), nlatent(l), ndata)),
    )
end

function init_aux_posterior(T::DataType, lik::LogisticSoftMaxLikelihoods, n::Int)
    nl = nlatent(lik) # This needs # https://github.com/JuliaGaussianProcesses/GPLikelihoods.jl/pull/68
    return For(
        TupleVector(;
            y=nestedview(falses(nl, n)),
            c=nestedview(zeros(T, nl, n)),
            p=nestedview(zeros(T, nl, n)),
        ),
    ) do q
        PolyaGammaNegativeMultinomial(q.y, q.c, q.p)
    end
end

function aux_full_conditional(
    lik::LogisticSoftMaxLikelihoods, y::AbstractVector{<:Bool}, f::AbstractVector{<:Real}
)
    return PolyaGammaNegativeMultinomial(
        y, abs.(f), _scale_σf(lik.invlink, f) / _sum_θ(lik.invlink)
    )
end

function aux_posterior!(
    qΩ,
    lik::BijectiveLogisticSoftMaxLikelihood,
    y::AbstractVector{<:AbstractVector{<:Bool}},
    qf::AbstractVector{<:AbstractVector{<:Normal}},
)
    φ = only(qΩ.inds)
    for (i, φᵢ) in enumerate(φ)
        φᵢ.c .= sqrt.(second_moment.(qf[i]))
        φᵢ.y .= y[i]
        φᵢ.p .=
            approx_expected_logistic.(-mean.(qf[i]), φᵢ.c) /
            (_get_const(lik.invlink) + nlatent(lik))
    end
    return qΩ
end

function aux_posterior!(
    qΩ,
    lik::LogisticSoftMaxLikelihood,
    y::AbstractVector{<:AbstractVector},
    qf::AbstractVector{<:AbstractVector{<:Normal}},
)
    φ = only(qΩ.inds)
    for (i, φᵢ) in enumerate(φ)
        @. φᵢ.c = sqrt(second_moment(qf[i]))
        φᵢ.y .= y[i]
        φᵢ.p .= approx_expected_logistic.(-mean.(qf[i]), φᵢ.c) ./ nlatent(lik)
    end
    return qΩ
end

function auglik_potential(::LogisticSoftMaxLikelihoods, Ω, y::ArrayOfSimilarArrays{<:Bool})
    return nestedview(((flatview(y) - flatview(Ω.n)) / 2)')
    # We want to have a Tuple of vector of the same size as the number of classes
end

function auglik_precision(::LogisticSoftMaxLikelihoods, Ω, ::AbstractVector)
    return transpose_nested(Ω.ω)
end

function expected_auglik_potential(
    ::LogisticSoftMaxLikelihoods, qΩ, y::ArrayOfSimilarArrays{<:Bool}
)
    return nestedview(((flatview(y) - flatview(tvmean(qΩ).n)) / 2)')
end

function expected_auglik_precision(::LogisticSoftMaxLikelihoods, qΩ, ::AbstractVector)
    return transpose_nested(tvmean(qΩ).ω)
end

function expected_auglik_potential_and_precision(
    ::LogisticSoftMaxLikelihoods, qΩ, y::AbstractVector
)
    θ = tvmean(qΩ)
    return nestedview(((flatview(y) - flatview(θ.n)) / 2)'), transpose_nested(θ.ω)
end

function logtilt(
    ::LogisticSoftMaxLikelihoods,
    (ω, n)::Tuple{<:AbstractVector{<:Real},<:AbstractVector{<:Integer}},
    y::AbstractVector{<:Integer},
    f::AbstractVector{<:Real},
)
    return -sum(y + n) * logtwo + sum((y - n) .* f - abs2.(f) .* ω) / 2
end

function aux_prior(lik::LogisticSoftMaxLikelihoods, y::ArrayOfSimilarArrays{<:Bool})
    return For(y) do yᵢ
        aux_prior(lik, yᵢ)
    end
end

function aux_prior(lik::BijectiveLogisticSoftMaxLikelihood, y::AbstractVector{<:Integer})
    return PolyaGammaNegativeMultinomial(
        y, zeros(Int, length(y)), fill(inv(_sum_θ(lik.invlink)), nlatent(lik))
    )
end
function aux_prior(lik::LogisticSoftMaxLikelihood, y::AbstractVector{<:Integer})
    # There is no proper prior for this so we just pretend this will work
    return PolyaGammaNegativeMultinomial(
        y, zeros(Int, length(y)), fill(1 / nlatent(lik), nlatent(lik))
    )
end

function expected_logtilt(
    ::BijectiveLogisticSoftMaxLikelihood, qω, y, qf::AbstractVector{<:Normal}
)
    θ = ntmean(qω)
    return -sum(y + θ.n) * logtwo +
           mapreduce(+, y, θ.n, θ.ω, mean.(qf), var.(qf)) do yᵢ, nᵢ, ωᵢ, mᵢ, vᵢ
        ((yᵢ - nᵢ) * mᵢ - (abs2(mᵢ) + vᵢ) * ωᵢ) / 2
    end
end


text:
struct WangLi3SMil_AConstantCache <: StochasticDiffEqConstantCache end
@cache struct WangLi3SMil_ACache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  k₁::rateType
  tmp::uType
end

alg_cache(alg::WangLi3SMil_A,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = WangLi3SMil_AConstantCache()

function alg_cache(alg::WangLi3SMil_A,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  k₁ = zero(rate_prototype)
  tmp = zero(u)
  WangLi3SMil_ACache(u,uprev,k,k₁,tmp)
end

struct WangLi3SMil_BConstantCache <: StochasticDiffEqConstantCache end
@cache struct WangLi3SMil_BCache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  k₁::rateType
  tmp::uType
end

alg_cache(alg::WangLi3SMil_B,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = WangLi3SMil_BConstantCache()

function alg_cache(alg::WangLi3SMil_B,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  k₁ = zero(rate_prototype)
  tmp = zero(u)
  WangLi3SMil_BCache(u,uprev,k,k₁,tmp)
end

struct WangLi3SMil_CConstantCache <: StochasticDiffEqConstantCache end
@cache struct WangLi3SMil_CCache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  k₁::rateType
  tmp::uType
end

alg_cache(alg::WangLi3SMil_C,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = WangLi3SMil_CConstantCache()

function alg_cache(alg::WangLi3SMil_C,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  k₁ = zero(rate_prototype)
  tmp = zero(u)
  WangLi3SMil_CCache(u,uprev,k,k₁,tmp)
end

struct WangLi3SMil_DConstantCache <: StochasticDiffEqConstantCache end
@cache struct WangLi3SMil_DCache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  k₁::rateType
  tmp::uType
end

alg_cache(alg::WangLi3SMil_D,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = WangLi3SMil_DConstantCache()

function alg_cache(alg::WangLi3SMil_D,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  k₁ = zero(rate_prototype)
  tmp = zero(u)
  WangLi3SMil_DCache(u,uprev,k,k₁,tmp)
end

struct WangLi3SMil_EConstantCache <: StochasticDiffEqConstantCache end
@cache struct WangLi3SMil_ECache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  k₁::rateType
  tmp::uType
end

alg_cache(alg::WangLi3SMil_E,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = WangLi3SMil_EConstantCache()

function alg_cache(alg::WangLi3SMil_E,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  k₁ = zero(rate_prototype)
  tmp = zero(u)
  WangLi3SMil_ECache(u,uprev,k,k₁,tmp)
end

struct WangLi3SMil_FConstantCache <: StochasticDiffEqConstantCache end
@cache struct WangLi3SMil_FCache{uType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  k₁::rateType
  tmp::uType
end

alg_cache(alg::WangLi3SMil_F,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = WangLi3SMil_FConstantCache()

function alg_cache(alg::WangLi3SMil_F,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  k₁ = zero(rate_prototype)
  tmp = zero(u)
  WangLi3SMil_FCache(u,uprev,k,k₁,tmp)
end