text:
function _entity_name(s::String)
    if s === "acq"
        return s, true
    elseif s === "label"
        return s, true
    elseif s === "space"
        return s, true
    elseif s === "flip"
        return s, true
    elseif s === "rec"
        return s, true
    elseif s === "sub"
        return s, true
    elseif s === "sample"
        return s, true
    elseif s === "proc"
        return s, true
    elseif s === "part"
        return s, true
    elseif s === "split"
        return s, true
    elseif s === "desc"
        return s, true
    elseif s === "trc"
        return s, true
    elseif s === "recording"
        return s, true
    elseif s === "res"
        return s, true
    elseif s === "dir"
        return s, true
    elseif s === "ses"
        return s, true
    elseif s === "task"
        return s, true
    elseif s === "mt"
        return s, true
    elseif s === "run"
        return s, true
    elseif s === "mod"
        return s, true
    elseif s === "den"
        return s, true
    elseif s === "inv"
        return s, true
    elseif s === "ce"
        return s, true
    elseif s === "echo"
        return s, true
    else
        return s, false
    end
end

struct Entity
    entity::String
    label::String

    global _unsafe_entity(e::String, l::String) = new(e, l)

    function Entity(e::String, l::String)
        e, check = _entity_name(e)
        check || error("$e is not a valid entity.")
        return _unsafe_entity(e, l)
    end
    Entity(e, l) = Entity(String(e), String(l))
    function Entity(s::AbstractString)
        ss = split(s, "-")
        length(ss) !== 2 && error("$s cannot be parsed to an entity.")
        return Entity(ss[1], ss[2])
    end
end

struct EntityBaseName
    entities::Vector{Entity}
    suffix::String
    extension::String

    EntityBaseName(e::Vector{Entity}, s::AbstractString, ex::AbstractString) = new(e, s, ex)
    function EntityBaseName(s::AbstractString)
        filename, extension = splitext(s)
        filename_split = split(filename, "_")
        N = length(filename_split)
        # FIXME should probably for no negativity here
        entities = Vector{Entity}(undef, N-1)
        @inbounds for i in 1:(N-1)
            entities[i] = Entity(filename_split[i])
        end
        return EntityBaseName(entities, @inbounds(filename_split[end]), extension)
    end
end


text:
using Format
using Test
using Printf
using Random

@testset "cformat"    begin include( "cformat.jl" ) end
@testset "fmtspec"    begin include( "fmtspec.jl" ) end
@testset "formatexpr" begin include( "formatexpr.jl" ) end
@testset "fmt"        begin include( "fmt.jl" ) end

ext: jl


text:
#
# G(x) := g(x) + ind[S](c(x))
#
# prox_{a G}(x0) = argmin { G(x) + 1/(2a) ||x-x0||^2 }
#                    x
#                = argmin { g(x) + 1/(2a) ||x-x0||^2 | c(x) in S }
#                    x
#
# minimize         1/(2a) ||x-x0||^2 + g(x)
# subject to       c(x) in S
#
# ProxOptiModel( prob, x0, a, [y0] )
#   f(x) = 1/(2a) ||x-x0||^2
#   g(x) = g(x)
#   c(x) = c(x)
#   S    = S

import LinearAlgebra.norm

export ProxOptiModel

# Type
mutable struct ProxOptiModel <: AbstractOptiModel
  meta::OptiModelMeta
  base::AbstractOptiModel
  a::Real
  x0::AbstractVector

  function ProxOptiModel( meta::OptiModelMeta, prob::TP, a::Real, x0::Tx ) where {TP <: AbstractOptiModel, Tx <: AbstractVector}
    new( meta, prob, a, x0 )
  end
end

# Constructor
function ProxOptiModel( prob::TP, a::Real; x0::Tx=prob.meta.x0, y0::Tx=prob.meta.y0, name=prob.meta.name * "-prox" ) where {TP <: AbstractOptiModel, Tx <: AbstractVector}
  @assert a > 0
  @lencheck prob.meta.nvar x0
  @lencheck prob.meta.ncon y0
  meta = OptiModelMeta( prob.meta.nvar, prob.meta.ncon, x0=x0, y0=y0, minimize=true, name=name )
  pr = ProxOptiModel(meta, prob, a, x0)
  finalizer(prob -> finalize(prob.base), pr)
  return pr
end

##########################################################
# Methods
##########################################################
# obj, grad!, cons!, jprod!, jtprod!, proj!, prox!, objprox!

function obj(pr::ProxOptiModel, x::AbstractVector)
  @lencheck pr.meta.nvar x
  return (0.5 / pr.a) * norm( x - pr.x0, 2 )^2
end

function grad!(pr::ProxOptiModel, x::AbstractVector, dfx::AbstractVector)
  @lencheck pr.meta.nvar x dfx
  dfx = (x - pr.x0) ./ pr.a
  return nothing
end

function cons!(pr::ProxOptiModel, x::AbstractVector, cx::AbstractVector)
  @lencheck pr.meta.nvar x
  @lencheck pr.meta.ncon cx
  cons!(pr.base, x, cx)
  return nothing
end

function jprod!(pr::ProxOptiModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector)
  @lencheck prob.meta.nvar x v
  @lencheck pr.meta.ncon Jv
  jprod!(pr.base, x, v, Jv)
  return nothing
end

function jtprod!(pr::ProxOptiModel, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector)
  @lencheck pr.meta.nvar x Jtv
  @lencheck pr.meta.ncon v
  jtprod!(pr.base, x, v, Jtv)
  return nothing
end

function proj!(pr::ProxOptiModel, cx::AbstractVector, px::AbstractVector)
  @lencheck pr.meta.ncon cx px
  proj!(pr.base, cx, px )
  return nothing
end

function prox!(pr::ProxOptiModel, x::AbstractVector, a::Real, z::AbstractVector)
  @lencheck pr.meta.nvar x z
  prox!(pr.base, x, a, z )
  return nothing
end

function objprox!(pr::ProxOptiModel, x::AbstractVector, a::Real, z::AbstractVector)
  @lencheck pr.meta.nvar x z
  gz = objprox!( pr.base, x, a, z )
  return gz
end


text:
export CplexSolver

type CplexMathProgModel <: AbstractLinearQuadraticModel
    inner::Model
    lazycb
    cutcb
    heuristiccb
    branchcb
    incumbentcb
    infocb
    solvetime::Float64
end

function CplexMathProgModel(;options...)
    env = Env()
    # set_param!(env, "CPX_PARAM_MIPCBREDLP", 0) # access variables in original problem, not presolved
    # set_param!(env, "CPX_PARAM_PRELINEAR", 0) # MAY NOT BE NECESSARY, only performs linear presolving so can recover original variables
    set_param!(env, "CPX_PARAM_SCRIND", 1) # output logs to stdout by default
    for (name,value) in options
        set_param!(env, string(name), value)
    end

    m = CplexMathProgModel(_Model(env), nothing, nothing, nothing, nothing, nothing, nothing, NaN)
    return m
end

immutable CplexSolver <: AbstractMathProgSolver
    options
end
CplexSolver(;kwargs...) = CplexSolver(kwargs)
LinearQuadraticModel(s::CplexSolver) = CplexMathProgModel(;s.options...)

ConicModel(s::CplexSolver) = LPQPtoConicBridge(LinearQuadraticModel(s))
supportedcones(::CplexSolver) = [:Free,:Zero,:NonNeg,:NonPos,:SOC]

function loadproblem!(m::CplexMathProgModel, filename::AbstractString)
   read_model(m.inner, filename)
   prob_type = get_prob_type(m.inner)
   if prob_type in [:MILP,:MIQP, :MIQCP]
      m.inner.has_int = true
   end
   if prob_type in [:QP, :MIQP, :QCP, :MIQCP]
      m.inner.has_qc = true
   end
end

function loadproblem!(m::CplexMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
  add_vars!(m.inner, float(obj), float(collb), float(colub))

  neginf = typemin(eltype(rowlb))
  posinf = typemax(eltype(rowub))

  rangeconstrs = any((rowlb .!= rowub) & (rowlb .> neginf) & (rowub .< posinf))
  if rangeconstrs
    warn("Julia Cplex interface doesn't properly support range (two-sided) constraints.")
    add_rangeconstrs!(m.inner, float(A), float(rowlb), float(rowub))
  else
    b = Array(Float64,length(rowlb))
    senses = Array(Cchar,length(rowlb))
    for i in 1:length(rowlb)
      if rowlb[i] == rowub[i]
        senses[i] = 'E'
        b[i] = rowlb[i]
      elseif rowlb[i] > neginf
        senses[i] = 'G'
        b[i] = rowlb[i]
      else
        @assert rowub[i] < posinf
        senses[i] = 'L'
        b[i] = rowub[i]
      end
    end
    add_constrs!(m.inner, float(A), senses, b)
  end

  set_sense!(m.inner, sense)
end

writeproblem(m::CplexMathProgModel, filename::AbstractString) = write_model(m.inner, filename)

getvarLB(m::CplexMathProgModel) = get_varLB(m.inner)
setvarLB!(m::CplexMathProgModel, l) = set_varLB!(m.inner, l)
getvarUB(m::CplexMathProgModel) = get_varUB(m.inner)
setvarUB!(m::CplexMathProgModel, u) = set_varUB!(m.inner, u)

# CPXchgcoef
getconstrLB(m::CplexMathProgModel) = get_constrLB(m.inner)
setconstrLB!(m::CplexMathProgModel, lb) = set_constrLB!(m.inner, lb)
getconstrUB(m::CplexMathProgModel) = get_constrUB(m.inner)
setconstrUB!(m::CplexMathProgModel, ub) = set_constrUB!(m.inner, ub)

getobj(m::CplexMathProgModel) = get_obj(m.inner)
setobj!(m::CplexMathProgModel, c) = set_obj!(m.inner, c)

addvar!(m::CplexMathProgModel, l, u, coeff) = add_var!(m.inner, [], [], l, u, coeff)
addvar!(m::CplexMathProgModel, constridx, constrcoef, l, u, coeff) = add_var!(m.inner, constridx, constrcoef, l, u, coeff)

function addconstr!(m::CplexMathProgModel, varidx, coef, lb, ub)
  neginf = typemin(eltype(lb))
  posinf = typemax(eltype(ub))

  rangeconstrs = any((lb .!= ub) & (lb .> neginf) & (ub .< posinf))
  if rangeconstrs
    warn("Julia Cplex interface doesn't properly support range (two-sided) constraints.")
    add_rangeconstrs!(m.inner, [0], varidx, float(coef), float(lb), float(ub))
  else
    if lb == ub
      rel = 'E'
      rhs = lb
    elseif lb > neginf
      rel = 'G'
      rhs = lb
    else
      @assert ub < posinf
      rel = 'L'
      rhs = ub
    end
    add_constrs!(m.inner, ivec([1]), ivec(varidx), fvec(coef), convert(Vector{Cchar},[rel]), fvec(vec(collect(rhs))))
  end
end

getconstrmatrix(m::CplexMathProgModel) = get_constr_matrix(m.inner)

setsense!(m::CplexMathProgModel, sense) = set_sense!(m.inner, sense)

getsense(m::CplexMathProgModel) = get_sense(m.inner)

numvar(m::CplexMathProgModel) = num_var(m.inner)
numconstr(m::CplexMathProgModel) = num_constr(m.inner) + num_qconstr(m.inner)
numlinconstr(m::CplexMathProgModel) = num_constr(m.inner)
numquadconstr(m::CplexMathProgModel) = num_qconstr(m.inner)

# optimize!(m::CplexMathProgModel) = optimize!(m.inner)

function optimize!(m::CplexMathProgModel)
    # set callbacks if present
    if m.lazycb != nothing
        setmathproglazycallback!(m)
    end
    if m.cutcb != nothing
        setmathprogcutcallback!(m)
    end
    if m.heuristiccb != nothing
        setmathprogheuristiccallback!(m)
    end
    if m.branchcb != nothing
        setmathprogbranchcallback!(m)
    end
    if m.incumbentcb != nothing
        setmathprogincumbentcallback!(m)
    end
    if m.infocb != nothing
        setmathproginfocallback!(m)
    end
    start = time()
    optimize!(m.inner)
    m.solvetime = time() - start
end

function status(m::CplexMathProgModel)
    ret = get_status(m.inner)
    return (if ret in [:CPX_STAT_OPTIMAL, :CPXMIP_OPTIMAL, :CPXMIP_OPTIMAL_TOL]
        :Optimal
    elseif ret in [:CPX_STAT_UNBOUNDED, :CPXMIP_UNBOUNDED]
        :Unbounded
    elseif ret in [:CPX_STAT_INFEASIBLE, :CPXMIP_INFEASIBLE]
        :Infeasible
    elseif ret in [:CPX_STAT_INForUNBD, :CPXMIP_INForUNBD]
        Base.warn_once("CPLEX reported infeasible or unbounded. Set CPX_PARAM_REDUCE=1 to check
                        infeasibility or CPX_PARAM_REDUCE=2 to check unboundedness.")
        :InfeasibleOrUnbounded
    elseif contains(string(ret), "TIME_LIM")
        :UserLimit
    else
        ret
    end)
end

getobjval(m::CplexMathProgModel)   = get_objval(m.inner)
getobjbound(m::CplexMathProgModel) = get_best_bound(m.inner)
getsolution(m::CplexMathProgModel) = get_solution(m.inner)
getconstrsolution(m::CplexMathProgModel) = get_constr_solution(m.inner)
getreducedcosts(m::CplexMathProgModel) = get_reduced_costs(m.inner)
getconstrduals(m::CplexMathProgModel) = get_constr_duals(m.inner)
getrawsolver(m::CplexMathProgModel) = m.inner
getnodecount(m::CplexMathProgModel) = get_node_count(m.inner)

const var_type_map = Dict(
  'C' => :Cont,
  'B' => :Bin,
  'I' => :Int,
  'S' => :SemiCont,
  'N' => :SemiInt
)

const rev_var_type_map = Dict(
  :Cont => 'C',
  :Bin => 'B',
  :Int => 'I',
  :SemiCont => 'S',
  :SemiInt => 'N'
)

function setvartype!(m::CplexMathProgModel, v::Vector{Symbol})
  target_int = all(x->isequal(x,:Cont), v)
  prob_type = get_prob_type(m.inner)
  if target_int && prob_type in [:LP,:QP,:QCP]
    return nothing
  else
    set_vartype!(m.inner, map(x->rev_var_type_map[x], v))
    return nothing
  end
end

function getvartype(m::CplexMathProgModel)
  if m.inner.has_int
    return map(x->var_type_map[x], get_vartype(m.inner))
  else
    return fill(:Cont, num_var(m.inner))
  end
end

function getsolvetime(m::CplexMathProgModel)
    return m.solvetime
end

getinfeasibilityray(m::CplexMathProgModel) = get_infeasibility_ray(m.inner)
getunboundedray(m::CplexMathProgModel) = get_unbounded_ray(m.inner)

getbasis(m::CplexMathProgModel) = get_basis(m.inner)

setwarmstart!(m::CplexMathProgModel, v) = set_warm_start!(m.inner, v)

addsos1!(m::CplexMathProgModel, idx, weight) = add_sos!(m.inner, :SOS1, idx, weight)
addsos2!(m::CplexMathProgModel, idx, weight) = add_sos!(m.inner, :SOS2, idx, weight)

######
# QCQP
######
addquadconstr!(m::CplexMathProgModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs) =
    add_qconstr!(m.inner,linearidx,linearval,quadrowidx,quadcolidx,quadval,sense,rhs)
setquadobj!(m::CplexMathProgModel,rowidx,colidx,quadval) = add_qpterms!(m.inner,rowidx,colidx,quadval)

######
# Data
######
function getdettime(m::CplexMathProgModel)
    tim = Array(Cdouble,1)
    stat = @cpx_ccall(getdettime, Cint, (Ptr{Void},Ptr{Cdouble}), m.inner.env.ptr, tim)
    if stat != 0
        error(CplexError(m.inner.env, stat).msg)
    end
    return tim[1]
end

getobjgap(m::CplexMathProgModel) = get_rel_gap(m.inner)

###########
# Callbacks
###########
export cbaddboundbranchup!,
       cbaddboundbranchdown!,
       setmathprogbranchcallback!,
       cbgetnodelb,
       cbgetnodeub,
       cbgetnodeobjval,
       cbgetnodesleft,
       cbgetmipiterations,
       cbgetfeasibility,
       cbgetgap,
       cbgetstarttime,
       cbgetdetstarttime,
       cbgettimestamp,
       cbgetdettimestamp,
       cbgetintfeas

abstract CplexCallbackData <: MathProgCallbackData

# set to nothing to clear callback
setlazycallback!(m::CplexMathProgModel,f) = (m.lazycb = f)
setcutcallback!(m::CplexMathProgModel,f) = (m.cutcb = f)
setheuristiccallback!(m::CplexMathProgModel,f) = (m.heuristiccb = f)
setbranchcallback!(m::CplexMathProgModel,f) = (m.branchcb = f)
setincumbentcallback!(m::CplexMathProgModel,f) = (m.incumbentcb = f)
setinfocallback!(m::CplexMathProgModel,f) = (m.infocb = f)

function cbgetmipsolution(d::CplexCallbackData)
    @assert d.state == :MIPSol
    n = num_var(d.cbdata.model)
    sol = Array(Cdouble, n)
    stat = @cpx_ccall(getcallbacknodex, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Cdouble},Cint,Cint),
                      d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, sol, 0, n-1)
    if stat != 0
        error(CplexError(d.cbdata.model.env, stat).msg)
    end
    return sol
end

function cbgetmipsolution(d::CplexCallbackData, sol::Vector{Cdouble})
    @assert d.state == :MIPSol
    stat = @cpx_ccall(getcallbacknodex, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Cdouble},Cint,Cint),
                      d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, sol, 0, length(sol)-1)
    if stat != 0
        error(CplexError(d.cbdata.model.env, stat).msg)
    end
    return nothing
end

function cbgetlpsolution(d::CplexCallbackData)
    @assert d.state == :MIPNode || d.state == :MIPBranch
    n = num_var(d.cbdata.model)
    sol = Array(Cdouble, n)
    stat = @cpx_ccall(getcallbacknodex, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Cdouble},Cint,Cint),
                      d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, sol, 0, n-1)
    if stat != 0
        error(CplexError(d.cbdata.model.env, stat).msg)
    end
    return sol
end

function cbgetlpsolution(d::CplexCallbackData, sol::Vector{Cdouble})
    @assert d.state == :MIPNode || d.state == :MIPIncumbent || d.state == :MIPBranch
    stat = @cpx_ccall(getcallbacknodex, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Cdouble},Cint,Cint),
                      d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, sol, 0, length(sol)-1)
    if stat != 0
        error(CplexError(d.cbdata.model.env, stat).msg)
    end
    return nothing
end

for (func,param,typ) in ((:cbgetexplorednodes,CPX_CALLBACK_INFO_NODE_COUNT_LONG,:Int64),
                         (:cbgetnodesleft,CPX_CALLBACK_INFO_NODES_LEFT_LONG,:Int64),
                         (:cbgetmipiterations,CPX_CALLBACK_INFO_MIP_ITERATIONS_LONG,:Int64),
                         (:cbgetbestbound,CPX_CALLBACK_INFO_BEST_REMAINING,:Cdouble),
                         (:cbgetobj,CPX_CALLBACK_INFO_BEST_INTEGER,:Cdouble),
                         (:cbgetgap,CPX_CALLBACK_INFO_MIP_REL_GAP,:Cdouble),
                         (:cbgetfeasibility,CPX_CALLBACK_INFO_MIP_FEAS,:Cint),
                         (:cbgetstarttime,CPX_CALLBACK_INFO_STARTTIME,:Cdouble),
                         (:cbgetdetstarttime,CPX_CALLBACK_INFO_STARTDETTIME,:Cdouble),
                         (:cbgettimestamp,CPX_CALLBACK_INFO_ENDTIME,:Cdouble),
                         (:cbgetdettimestamp,CPX_CALLBACK_INFO_ENDDETTIME,:Cdouble))
    @eval begin
        function $(func)(d::CplexCallbackData)
            val = Array($(typ),1)
            ret = @cpx_ccall(getcallbackinfo, Cint, (Ptr{Void},Ptr{Void},Cint,Cint,Ptr{Void}),
                              d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, $(convert(Cint,param)), val)
            if ret != 0
                error(CplexError(d.cbdata.model.env, stat).msg)
            end
            return val[1]
        end
    end
end

# returns :MIPNode :MIPSol :Intermediate
cbgetstate(d::CplexCallbackData) = d.state

#const sensemap = Dict('=' => 'E', '<' => 'L', '>' => 'G')
function cbaddcut!(d::CplexCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode
    cbcut(d.cbdata, d.where, convert(Vector{Cint}, varidx), float(varcoef), sensemap[sense], float(rhs))
    unsafe_store!(d.userinteraction_p, convert(Cint,CPX_CALLBACK_ABORT_CUT_LOOP), 1)
end

function cbaddcutlocal!(d::CplexCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode
    cbcutlocal(d.cbdata, d.where, convert(Vector{Cint}, varidx), float(varcoef), sensemap[sense], float(rhs))
    unsafe_store!(d.userinteraction_p, convert(Cint,CPX_CALLBACK_ABORT_CUT_LOOP), 1)
end

function cbaddlazy!(d::CplexCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode || d.state == :MIPSol
    cblazy(d.cbdata, d.where, convert(Vector{Cint}, varidx), float(varcoef), sensemap[sense], float(rhs))
end

function cbaddlazylocal!(d::CplexCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode || d.state == :MIPSol
    cblazylocal(d.cbdata, d.where, convert(Vector{Cint}, varidx), float(varcoef), sensemap[sense], float(rhs))
end

function cbaddsolution!(d::CplexCallbackData)
    val = pointer_to_array(d.userinteraction_p, 1)
    if val[1] == CPX_CALLBACK_SET
        error("CPLEX only allows one heuristic solution for each call to the callback")
    end
    unsafe_store!(d.userinteraction_p, convert(Cint,CPX_CALLBACK_SET), 1)
end

function cbsetsolutionvalue!(d::CplexCallbackData,varidx,value)
    @assert 1 <= varidx <= num_var(d.cbdata.model)
    d.heur_x[varidx] = value
end

function cbaddboundbranchup!(d::CplexCallbackData,idx,bd,nodeest)
    seqnum = cbbranch(d.cbdata, d.where,convert(Cint,idx-1),convert(Cchar,'L'),bd,nodeest)
    unsafe_store!(d.userinteraction_p, convert(Cint,CPX_CALLBACK_SET), 1)
    seqnum
end

function cbaddboundbranchdown!(d::CplexCallbackData,idx,bd,nodeest)
    seqnum = cbbranch(d.cbdata, d.where,convert(Cint,idx-1),convert(Cchar,'U'),bd,nodeest)
    unsafe_store!(d.userinteraction_p, convert(Cint,CPX_CALLBACK_SET), 1)
    seqnum
end

function cbaddconstrbranch!(d::CplexCallbackData, indices, coeffs, rhs, sense, nodeest)
    seqnum = cbbranchconstr(d.cbdata,
                   d.where,
                   Cint[idx-1 for idx in indices],
                   Cdouble[c for c in coeffs],
                   convert(Cdouble, rhs),
                   convert(Cchar, sense),
                   nodeest)
    unsafe_store!(d.userinteraction_p, convert(Cint,CPX_CALLBACK_SET), 1)
    seqnum
end

function cbprocessincumbent!(d::CplexCallbackData,accept::Bool)
    if accept
        unsafe_store!(d.isfeas_p, convert(Cint, 1), 1)
    else
        unsafe_store!(d.isfeas_p, convert(Cint, 0), 1)
    end
    nothing
end

type CplexLazyCallbackData <: CplexCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    userinteraction_p::Ptr{Cint}
end

type CplexCutCallbackData <: CplexCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    userinteraction_p::Ptr{Cint}
end

terminate(model::CplexMathProgModel) = terminate(model.inner)

# breaking abstraction, define our low-level callback to eliminate
# a level of indirection
function mastercallback(env::Ptr{Void}, cbdata::Ptr{Void}, wherefrom::Cint, userdata::Ptr{Void}, userinteraction_p::Ptr{Cint})
    model = unsafe_pointer_to_objref(userdata)::CplexMathProgModel
    cpxrawcb = CallbackData(cbdata, model.inner)
    if wherefrom == CPX_CALLBACK_MIP_CUT_FEAS || wherefrom == CPX_CALLBACK_MIP_CUT_UNBD
        state = :MIPSol
    # elseif wherefrom == CPX_CALLBACK_MIP_CUT_LOOP || wherefrom == CPX_CALLBACK_MIP_CUT_LAST
    elseif wherefrom == CPX_CALLBACK_MIP_CUT_LAST
        state = :MIPNode
    else
        state = :Intermediate
    end

    if model.infocb != nothing
        cpxcb = CplexInfoCallbackData(cpxrawcb, state, wherefrom)
        stat = model.infocb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    if model.lazycb != nothing && state == :MIPSol
        cpxcb = CplexLazyCallbackData(cpxrawcb, state, wherefrom, userinteraction_p)
        stat = model.lazycb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    if model.cutcb != nothing && state == :MIPNode
        cpxcb = CplexCutCallbackData(cpxrawcb, state, wherefrom, userinteraction_p)
        stat = model.cutcb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    return convert(Cint, 0)
end

type CplexHeuristicCallbackData <: CplexCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    sol::Vector{Float64}
    heur_x::Vector{Float64}
    isfeas_p::Ptr{Cint}
    userinteraction_p::Ptr{Cint}
end

function masterheuristiccallback(env::Ptr{Void},
                                 cbdata::Ptr{Void},
                                 wherefrom::Cint,
                                 userdata::Ptr{Void},
                                 objval_p::Ptr{Cdouble},
                                 xx::Ptr{Cdouble},
                                 isfeas_p::Ptr{Cint},
                                 userinteraction_p::Ptr{Cint})
    model = unsafe_pointer_to_objref(userdata)::CplexMathProgModel
    cpxrawcb = CallbackData(cbdata, model.inner)
    if wherefrom == CPX_CALLBACK_MIP_HEURISTIC
        state = :MIPNode
    else
        state = :Intermediate
    end

    if model.infocb != nothing
        cpxcb = CplexInfoCallbackData(cpxrawcb, state, wherefrom)
        stat = model.infocb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    if model.heuristiccb != nothing && state == :MIPNode
        sol = pointer_to_array(xx, numvar(model))
        cpxcb = CplexHeuristicCallbackData(cpxrawcb, state, wherefrom, sol, fill(NaN, numvar(model)), isfeas_p, userinteraction_p)
        stat = model.heuristiccb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
        if any(x->!isnan(x), cpxcb.heur_x) # we filled in some solution values
            unsafe_store!(objval_p, dot(get_obj(model.inner), cpxcb.heur_x), 1)
            for i in 1:numvar(model)
                unsafe_store!(xx, cpxcb.heur_x[i], i)
            end
            if any(x->isnan(x), cpxcb.heur_x) # we have a partial solution
                unsafe_store!(isfeas_p, convert(Cint,CPX_ON),  1)
            else
                unsafe_store!(isfeas_p, convert(Cint,CPX_OFF), 1)
            end
        end
    end
    return convert(Cint, 0)
end

function setmathproglazycallback!(model::CplexMathProgModel)
    set_param!(model.inner.env, "CPX_PARAM_MIPCBREDLP", 0)
    set_param!(model.inner.env, "CPX_PARAM_PRELINEAR", 0)
    set_param!(model.inner.env, "CPX_PARAM_REDUCE", CPX_PREREDUCE_PRIMALONLY)
    cpxcallback = cfunction(mastercallback, Cint, (Ptr{Void}, Ptr{Void}, Cint, Ptr{Void}, Ptr{Cint}))
    stat = @cpx_ccall(setlazyconstraintcallbackfunc, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    nothing
end

function setmathprogcutcallback!(model::CplexMathProgModel)
    set_param!(model.inner.env, "CPX_PARAM_MIPCBREDLP", 0)
    set_param!(model.inner.env, "CPX_PARAM_PRELINEAR", 0)
    cpxcallback = cfunction(mastercallback, Cint, (Ptr{Void}, Ptr{Void}, Cint, Ptr{Void}, Ptr{Cint}))
    stat = @cpx_ccall(setusercutcallbackfunc, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    nothing
end

function setmathprogheuristiccallback!(model::CplexMathProgModel)
    set_param!(model.inner.env, "CPX_PARAM_MIPCBREDLP", 0)
    set_param!(model.inner.env, "CPX_PARAM_PRELINEAR", 0)
    cpxcallback = cfunction(masterheuristiccallback, Cint, (Ptr{Void}, Ptr{Void}, Cint, Ptr{Void}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}))
    stat = @cpx_ccall(setheuristiccallbackfunc, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    nothing
end

immutable BranchingChoice
    indices::Vector{Cint}
    bounds::Vector{Cdouble}
    lu::Vector{Cchar}
end
export BranchingChoice

type CplexBranchCallbackData <: CplexCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    userinteraction_p::Ptr{Cint}
    nodes::Vector{BranchingChoice}
end

function masterbranchcallback(env::Ptr{Void},
                              cbdata::Ptr{Void},
                              wherefrom::Cint,
                              userdata::Ptr{Void},
                              typ::Cint,
                              sos::Cint,
                              nodecnt::Cint,
                              bdcnt::Cint,
                              nodebeg::Ptr{Cint},
                              indices::Ptr{Cint},
                              lu::Ptr{Cchar},
                              bd::Ptr{Cdouble},
                              nodeest::Ptr{Cdouble},
                              userinteraction_p::Ptr{Cint})
    model = unsafe_pointer_to_objref(userdata)::CplexMathProgModel
    cpxrawcb = CallbackData(cbdata, model.inner)
    if wherefrom == CPX_CALLBACK_MIP_BRANCH
        @assert 0 <= nodecnt <= 2
        state = :MIPBranch
    else
        state = :Intermediate
    end

    if model.infocb != nothing
        cpxcb = CplexInfoCallbackData(cpxrawcb, state, wherefrom)
        stat = model.infocb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    if model.branchcb != nothing && state == :MIPBranch
        numbranchingvars = pointer_to_array(nodebeg, convert(Cint,nodecnt))::Vector{Cint} + 1
        idxs = pointer_to_array(indices, sum(numbranchingvars))::Vector{Cint}
        vals = pointer_to_array(bd, sum(numbranchingvars))::Vector{Cdouble}
        dirs = pointer_to_array(bd, sum(numbranchingvars))::Vector{Cdouble}
        nodes = Array(BranchingChoice, nodecnt)
        if nodecnt >= 1
            subidx = 1 : (numbranchingvars[1])
            nodes[1] = BranchingChoice(idxs[subidx], vals[subidx], dirs[subidx])
        end
        if nodecnt == 2
            subidx = (numbranchingvars[1]+1) : (numbranchingvars[2])
            nodes[2] = BranchingChoice(idxs[subidx], vals[subidx], dirs[subidx])
        end
        cpxcb = CplexBranchCallbackData(cpxrawcb, state, wherefrom, userinteraction_p, nodes)
        stat = model.branchcb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    return convert(Cint, 0)
end

function setmathprogbranchcallback!(model::CplexMathProgModel)
    set_param!(model.inner.env, "CPX_PARAM_MIPCBREDLP", 0)
    set_param!(model.inner.env, "CPX_PARAM_PRELINEAR", 0)
    set_param!(model.inner.env, "CPX_PARAM_REDUCE", CPX_PREREDUCE_PRIMALONLY)
    cpxcallback = cfunction(masterbranchcallback, Cint, (Ptr{Void},
                                                         Ptr{Void},
                                                         Cint,
                                                         Ptr{Void},
                                                         Cint,
                                                         Cint,
                                                         Cint,
                                                         Cint,
                                                         Ptr{Cint},
                                                         Ptr{Cint},
                                                         Ptr{Cchar},
                                                         Ptr{Cdouble},
                                                         Ptr{Cdouble},
                                                         Ptr{Cint}))
    stat = @cpx_ccall(setbranchcallbackfunc, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    nothing
end

type CplexIncumbentCallbackData <: CplexCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    sol::Vector{Float64}
    isfeas_p::Ptr{Cint}
    userinteraction_p::Ptr{Cint}
    nodes::Vector{BranchingChoice}
end

function masterincumbentcallback(env::Ptr{Void},
                                 cbdata::Ptr{Void},
                                 wherefrom::Cint,
                                 userdata::Ptr{Void},
                                 objval::Cdouble,
                                 xx::Ptr{Cdouble},
                                 isfeas_p::Ptr{Cint},
                                 useraction_p::Ptr{Cint})
    model = unsafe_pointer_to_objref(userdata)::CplexMathProgModel
    cpxrawcb = CallbackData(cbdata, model.inner)
    if wherefrom == CPX_CALLBACK_MIP_INCUMBENT_NODESOLN ||
       wherefrom == CPX_CALLBACK_MIP_INCUMBENT_HEURSOLN ||
       wherefrom == CPX_CALLBACK_MIP_INCUMBENT_USERSOLN
        state = :MIPIncumbent
    else
        state = :Intermediate
    end

    if model.infocb != nothing
        cpxcb = CplexInfoCallbackData(cpxrawcb, state, wherefrom)
        stat = model.infocb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    if model.incumbentcb != nothing && state == :MIPIncumbent
        sol = pointer_to_array(xx, numvar(model))
        cpxcb = CplexIncumbentCallbackData(cpxrawcb, state, wherefrom, sol, isfeas_p, useraction_p, BranchingChoice[])
        stat = model.incumbentcb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    return convert(Cint, 0)
end

function setmathprogincumbentcallback!(model::CplexMathProgModel)
    set_param!(model.inner.env, "CPX_PARAM_MIPCBREDLP", 0)
    set_param!(model.inner.env, "CPX_PARAM_PRELINEAR", 0)
    set_param!(model.inner.env, "CPX_PARAM_REDUCE", CPX_PREREDUCE_PRIMALONLY)
    cpxcallback = cfunction(masterincumbentcallback, Cint, (Ptr{Void},
                                                            Ptr{Void},
                                                            Cint,
                                                            Ptr{Void},
                                                            Cdouble,
                                                            Ptr{Cdouble},
                                                            Ptr{Cint},
                                                            Ptr{Cint}))
    stat = @cpx_ccall(setincumbentcallbackfunc, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    nothing
end

type CplexInfoCallbackData <: CplexCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
end

function masterinfocallback(env::Ptr{Void},
                            cbdata::Ptr{Void},
                            wherefrom::Cint,
                            userdata::Ptr{Void})
    model = unsafe_pointer_to_objref(userdata)::CplexMathProgModel
    if model.infocb != nothing
        state = :Intermediate
        cpxrawcb = CallbackData(cbdata, model.inner)
        cpxcb = CplexInfoCallbackData(cpxrawcb, state, wherefrom)
        stat = model.infocb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
    end
    return convert(Cint, 0)
end

function setmathproginfocallback!(model::CplexMathProgModel)
    cpxcallback = cfunction(masterinfocallback, Cint, (Ptr{Void}, Ptr{Void}, Cint, Ptr{Void}))
    stat = @cpx_ccall(setinfocallbackfunc, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    nothing
end

function cbgetnodelb(d::CplexCallbackData)
    n = num_var(d.cbdata.model)
    lb = Array(Cdouble,n)
    stat = @cpx_ccall(getcallbacknodelb, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Cdouble},Cint,Cint),
                      d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, lb, 0, n-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return lb
end

function cbgetnodeub(d::CplexCallbackData)
    n = num_var(d.cbdata.model)
    ub = Array(Cdouble,n)
    stat = @cpx_ccall(getcallbacknodeub, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Cdouble},Cint,Cint),
                      d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, ub, 0, n-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return ub
end

function cbgetnodeobjval(d::CplexCallbackData)
    val = Array(Cdouble,1)
    stat = @cpx_ccall(getcallbacknodeobjval, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Cdouble}),
                      d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, val)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return val[1]
end

function cbgetintfeas(d::CplexCallbackData)
    n = num_var(d.cbdata.model)
    feas = Array(Cint,n)
    stat = @cpx_ccall(getcallbacknodeintfeas, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Cint},Cint,Cint),
                      d.cbdata.model.env.ptr, d.cbdata.cbdata, d.where, feas, 0, n-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return convert(Vector{Int64},feas)
end


text:
import Aswing

# Read and input .asw file
inputasw = "hawk.asw"
asw = Aswing.readasw(inputasw)
Aswing.setgeom(asw)

# Read and input .pnt file
inputpnt = "hawk.pnt"
pnt = Aswing.readpnt(inputpnt)
Aswing.setpnt(pnt)

# Save .asw and .pnt files, can be used with original ASWING
outputasw = "out.asw"
Aswing.aswsav(outputasw)
outputpnt = "out.pnt"
Aswing.pntsav(outputpnt)

# ASWING files can also be read using the following functions
newasw = Aswing.readasw(outputasw)
newpnt = Aswing.readpnt(outputpnt)

# solves for first operating point
steady_solution = Aswing.solvesteady(1)

# gets stability derivatives of first operating point
stabderivs = Aswing.getstabilityderivatives(1)

# you can also get and/or set the paneling between break points
newintervals = [Aswing.getintervals(i) for i  in 1:4]*2
Aswing.nodeset(newintervals)
steady_solution = Aswing.solvesteady(1)

# quasi-steady and rigid static margins can be calculated, but the
# geometry and operating point conditions will need to be re-input afterwards
sm_qs, sm_rigid, fail = Aswing.getstaticmargin(1)

Aswing.setgeom(asw)
Aswing.setpnt(pnt)
steady_solution = Aswing.solvesteady(1)

# note that when you input the geometry again the paneling is reset to the defaults
intervals = [Aswing.getintervals(i) for i  in 1:4]

# you can also do unsteady solves
Aswing.retainparam()
Aswing.setcons(Aswing.Constraints("free"))
unsteady_solution = Aswing.solveunsteady(0.1, 10)

# and you can get eigenvalues
asweigs = Aswing.geteigs(16, -0.1+0.0im, 1)


text:
using Test

include("../src/main_model.jl")
m = getpage("NDCs", false, false)
run(m)

while m[:Discontinuity,:occurdis_occurrencedummy] != [0.,0.,0.,0.,0.,0.,0.,1.,1.,1.]
    run(m)
end

# Climate component
temp = m[:ClimateTemperature, :rt_g_globaltemperature]
temp_compare = readpagedata(m,"test/validationdata/noperm/rt_g_globaltemperature.csv")
@test temp ≈ temp_compare rtol=1e-4

# Abatement Costs
abatement = m[:TotalAbatementCosts, :tct_totalcosts]
abatement_compare = readpagedata(m,"test/validationdata/noperm/tct_totalcosts.csv")
@test abatement ≈ abatement_compare rtol=1e-2

te = m[:EquityWeighting, :te_totaleffect]
te_compare = 9.462179758012009e8
@test te ≈ te_compare rtol=1e4

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
function scent(f::Function,r::Real)
    if f(r)
        i = st_index[1]
        if i <= nstackframes
            ST[:,i] = stacktrace()[2:4]    # 1 is the line where scent is called
            if st_index[1] < nstackframes
                st_index[1] += 1
            else
                st_index[1] = 1
                @warn "DrWatson reached $nstackframes, the maximum number of stackframes."
            end
        end
    end
end

empty_symbol = Symbol("")
const unknown_st = Base.StackTraces.StackFrame(empty_symbol,empty_symbol,-1,nothing,true,false,0)

const nstackframes = 10000      # quite arbitrary for the moment
const stackframe_depth = 3

const ST = [unknown_st for _ in 1:stackframe_depth, _ in 1:nstackframes]
const st_index = [1]

get_stacktrace(i::Int=1) = ST[:,i]
get_stacktraces() = ST[:,1:st_index[1]-1]
get_st_index() = st_index[1]


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

using LinearAlgebra, Symbolics, DifferentialEquations, JLD2

struct DoublePendulum{T}
	m1::T
	m2::T
	l1::T
	l2::T
end

n = 2 # number of states
model = DoublePendulum(1.0, 1.0, 1.0, 1.0) # model

# kinematics
function kinematics_1(model::DoublePendulum, q)
	θ1, θ2 = q

	[0.5 * model.l1 * sin(θ1);
	 -0.5 * model.l1 * cos(θ1)]
end

function kinematics_2(model::DoublePendulum, q)
	θ1, θ2 = q

	[model.l1 * sin(θ1) + 0.5 * model.l2 * sin(θ1 + θ2);
	 -model.l1 * cos(θ1) - 0.5 * model.l2 * cos(θ1 + θ2)]
end

# fast kinematics functions
@variables q[1:n]
@variables q̇[1:n]

k1 = kinematics_1(model, q)
k2 = kinematics_2(model, q)

k1_exp = Symbolics.build_function(k1, q)
k2_exp = Symbolics.build_function(k2, q)

k1_func = eval(k1_exp[1])
k2_func = eval(k2_exp[1])

# kinematics Jacobians
j1 = Symbolics.jacobian(k1, q, simplify = true)
j2 = Symbolics.jacobian(k2, q, simplify = true)

j1_exp = Symbolics.build_function(j1, q)
j2_exp = Symbolics.build_function(j2, q)

j1_func = eval(j1_exp[1])
j2_func = eval(j2_exp[1])

# Lagrangian
function lagrangian(model, q, q̇)
	L = 0.0

	# mass 1
	v1 = j1_func(q) * q̇
	L += 0.5 * model.m1 * transpose(v1) * v1 		# kinetic energy
	L -= model.m1 * 9.81 * k1_func(q)[2]            # potential energy

	# mass 2
	v2 = j2_func(q) * q̇
	L += 0.5 * model.m2 * transpose(v2) * v2
	L -= model.m2 * 9.81 * k2_func(q)[2]

	return L
end

# fast Lagrangian
L = lagrangian(model, q, q̇)

dLq = Symbolics.gradient(L, q, simplify = true)
dLq̇ = Symbolics.gradient(L, q̇, simplify = true)
ddL = Symbolics.hessian(L, [q; q̇], simplify = true)

# mass matrix
M = ddL[n .+ (1:n), n .+ (1:n)]
M = simplify.(M)

# dynamics bias
C = ddL[n .+ (1:n), 1:n] * q̇ - dLq
C = simplify.(C)

# dynamics
# ẋ = [q̇; M \ (-1.0 * C)]
ẋ = [q̇; M \ (-0.5 * q̇ -1.0 * C)]
ẋ = simplify.(ẋ)

ẋ_exp = Symbolics.build_function(ẋ, q, q̇)
dynamics = eval(ẋ_exp[1])

# save dynamics function
path = joinpath(pwd(), "dynamics.jld2")
# @save path ẋ_exp
# @load path ẋ_exp

# DifferentialEquations.jl
function dynamics!(ẋ, x, p, t)
	ẋ .= dynamics(view(x, 1:n), view(x, n .+ (1:n)))
end

# simulate
x0 = [0.5 * π; 0.0; 0.0; 0.0]
tspan = (0.0, 10.0)
dt = 0.01
prob = ODEProblem(dynamics!, x0, tspan)
sol = solve(prob, Tsit5(), adaptive = false, dt = dt)

# MeshCat.jl
include(joinpath(pwd(), "visuals.jl"))
vis = Visualizer()
render(vis)
visualize_double_pendulum!(vis, model, sol.u, Δt = dt)


text:
@testset "Compressor/Decompressor" begin
    address(x) = UInt(Base.unsafe_convert(Ptr{Nothing}, x))

    for T in [Decompressor, Compressor]
        a = T()
        b = T()
        c = T()

        @test pointer_from_objref(a) != pointer_from_objref(c)
        @test pointer_from_objref(a) != pointer_from_objref(b)
        @test pointer_from_objref(b) != pointer_from_objref(c)
        @test address(a) != address(c)
        @test address(a) != address(b)
        @test address(b) != address(c)
    end

    c = Compressor()
    @test c.level == Compressor(6).level
end

@testset "Errors" begin
    # No space for decompression
    v = Vector{UInt8}("Hello, there!")
    c = Compressor()
    d = Decompressor()
    @test_throws LibDeflate.LibDeflateError(LibDeflate.OUTPUT_BUFFER_TOO_SMALL) compress!(c, zeros(UInt8, 16), v)

    # Not compressed data
    @test_throws LibDeflate.LibDeflateError(LibDeflate.BAD_DATA) decompress!(d, zeros(UInt8, 512), rand(UInt8, 32))

    # Decompressed data too short
    v = zeros(UInt8, 256)
    bytes = compress!(c, v, Vector{UInt8}("ABC"^51))
    compressed = v[1:bytes]
    @test_throws LibDeflate.LibDeflateError(LibDeflate.INSUFFICIENT_SPACE) decompress!(d, zeros(UInt8, 1024), compressed, 150)

    # Decompressed data too long
    @test_throws LibDeflate.LibDeflateError(LibDeflate.INSUFFICIENT_SPACE) decompress!(d, zeros(UInt8, 32), compressed)
    @test_throws LibDeflate.LibDeflateError(LibDeflate.SHORT_INPUT) decompress!(d, zeros(UInt8, 1024), compressed, 160)
end

@testset "Compression" begin
    COMPRESSIBLE = [
        vcat(rand(UInt8, 412), zeros(UInt8, 100)),
        rand(1:1000, 100),
        join(rand((['A', 'C', 'G', 'T']), 500)),
        ("Na " * "na " ^ 15 * "Batman! ")^2
    ]
    outbuffer = zeros(UInt8, 512)
    for i in COMPRESSIBLE
        v = unsafe_wrap(Array, Ptr{UInt8}(pointer(i)), sizeof(i))
        bytes = compress!(Compressor(), outbuffer, v)
        @test bytes < length(v)
    end
end

# Unsafe CRC is implicitly tested by decompressing gzip with
# codeczlib. So we can just compare it to the unsafe one
@testset "Safe CRC" begin
    for testdata in ["", "foo", "abracadabra!"]
        @test crc32(collect(codeunits(testdata))) == unsafe_crc32(pointer(testdata), ncodeunits(testdata))
    end
end

@testset "Round trip" begin
    INPUT_DATA = [
        "",
        "Abracadabra!",
        "A man, a plan, a canal, Panama!",
        "No, no, no, no, no, no, no, no, no, no, no!",
        "sXXbYltTe]EDP`kRNUoEPVRnkq]gS^cquEv^BVTwAhtjFGGQBC",
        rand(UInt8, 2048),
    ]
    outbuffer = Vector{UInt8}(undef, 4096)
    unsafe_outbuffer = similar(outbuffer)
    backbuffer1 = similar(outbuffer)
    backbuffer2 = similar(outbuffer)
    unsafe_backbuffer1 = similar(outbuffer)
    unsafe_backbuffer2 = similar(outbuffer)

    compressor = Compressor()
    decompressor = Decompressor()

    for i in INPUT_DATA
        v = Vector{UInt8}(i)

        c_bytes_unsafe = unsafe_compress!(compressor, pointer(unsafe_outbuffer),
                                          length(outbuffer), pointer(v), length(v))
        c_bytes_safe = compress!(compressor, outbuffer, v)

        @test c_bytes_unsafe == c_bytes_safe
        @test unsafe_outbuffer[1:c_bytes_unsafe] == outbuffer[1:c_bytes_safe]

        d_bytes_unsafe1 = unsafe_decompress!(Base.HasLength(), decompressor,
                                            pointer(unsafe_backbuffer1), length(v),
                                            pointer(unsafe_outbuffer), c_bytes_unsafe)

        d_bytes_unsafe2 = unsafe_decompress!(Base.SizeUnknown(), decompressor,
                                             pointer(unsafe_backbuffer2), length(unsafe_backbuffer2),
                                             pointer(unsafe_outbuffer), c_bytes_unsafe)

        d_bytes_safe1 = decompress!(decompressor, backbuffer1, outbuffer, length(v))
        d_bytes_safe2 = decompress!(decompressor, backbuffer2, outbuffer)

        @test d_bytes_safe1 == d_bytes_safe2 == d_bytes_unsafe1 == d_bytes_unsafe2

        @test v == backbuffer1[1:d_bytes_safe1] == backbuffer2[1:d_bytes_safe1] ==
              unsafe_backbuffer1[1:d_bytes_safe1] == unsafe_backbuffer2[1:d_bytes_safe1]
    end
end


text:
function rescaling(x, lb, ub)
    @assert all([x <= 1, x >= 0])
    lb + x * (ub-lb)
end

function create_parameters(number_holes, number_training_data, δ, lb, ub)
    ps_vec = zeros(number_holes)
    for it_training = 1:number_training_data
        cur_ps = rescaling.(rand(number_holes), lb, ub)
        ps_vec = vcat(ps_vec, cur_ps)
        for it_hole = 1:number_holes
            perturb = zeros(number_holes)
            perturb[it_hole]+= δ
            ps_vec = vcat(ps_vec, cur_ps .+ perturb)
        end
    end
    return transpose(reshape(ps_vec[number_holes+1:end], (number_holes, number_training_data*6)))
end

function finite_difference_deriv(curv, vperturb, δ)
    return (vperturb-curv)/δ
end

function get_val_gradient_from_data(data_array; δ=1e-2)
    n_variable = size(data_array)[2]-2
    n_data_val_gradient = size(data_array)[1]÷(1+n_variable)

    real_val_gradient = zeros(Float64, (n_data_val_gradient, 2*n_variable+1))
    imag_val_gradient = zeros(Float64, (n_data_val_gradient, 2*n_variable+1))

    for i=1:n_data_val_gradient
        real_vector_data = data_array[(i-1)*(1+n_variable)+1, 1:n_variable] #push the data input
        real_curval = data_array[(i-1)*(1+n_variable)+1, n_variable+1]
        push!(real_vector_data, real_curval)

        imag_vector_data = data_array[(i-1)*(1+n_variable)+1, 1:n_variable]
        imag_curval = data_array[(i-1)*(1+n_variable)+1, n_variable+2]
        push!(imag_vector_data, imag_curval)

        for j=1:n_variable
            push!(real_vector_data, finite_difference_deriv(real_curval, data_array[(i-1)*(1+n_variable)+1+j, n_variable+1], δ))
            push!(imag_vector_data, finite_difference_deriv(imag_curval, data_array[(i-1)*(1+n_variable)+1+j, n_variable+2], δ))
        end
        real_val_gradient[i,:] = real_vector_data
        imag_val_gradient[i,:] = imag_vector_data
    end
    return real_val_gradient, imag_val_gradient
end


text:
using CSV, DataFrames, StatsBase
using StanSample, StatsPlots
using StatisticalRethinkingPlots, StatisticalRethinking

df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
df[!, :M] = zscore(df.Marriage)
df[!, :A] = zscore(df.MedianAgeMarriage)
df[!, :D] = zscore(df.Divorce)
data = (N=size(df, 1), D=df.D, A=df.A, M=df.M)

stan5_1 = "
data {
    int < lower = 1 > N; // Sample size
    vector[N] D; // Outcome
    vector[N] A; // Predictor
}
parameters {
    real a; // Intercept
    real bA; // Slope (regression coefficients)
    real < lower = 0 > sigma;    // Error SD
}
transformed parameters {
    vector[N] mu;               // mu is a vector
    for (i in 1:N)
        mu[i] = a + bA * A[i];
}
model {
    a ~ normal(0, 0.2);         //Priors
    bA ~ normal(0, 0.5);
    sigma ~ exponential(1);
    D ~ normal(mu , sigma);     // Likelihood
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

stan5_2 = "
data {
    int N;
    vector[N] D;
    vector[N] M;
}
parameters {
    real a;
    real bM;
    real<lower=0> sigma;
}
transformed parameters {
    vector[N] mu;
    for (i in 1:N)
        mu[i]= a + bM * M[i];

}
model {
    a ~ normal( 0 , 0.2 );
    bM ~ normal( 0 , 0.5 );
    sigma ~ exponential( 1 );
    D ~ normal( mu , sigma );
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

stan5_3 = "
data {
  int N;
  vector[N] D;
  vector[N] M;
  vector[N] A;
}
parameters {
  real a;
  real bA;
  real bM;
  real<lower=0> sigma;
}
transformed parameters {
    vector[N] mu;
    for (i in 1:N)
        mu[i] = a + bA * A[i] + bM * M[i];
}
model {
  a ~ normal( 0 , 0.2 );
  bA ~ normal( 0 , 0.5 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  D ~ normal( mu , sigma );
}
generated quantities{
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

fname = joinpath(@__DIR__, "coeftab_plot.png")
isfile(fname) && rm(fname)

m5_1s = SampleModel("m5.1s", stan5_1)
rc5_1s = stan_sample(m5_1s; data)

m5_2s = SampleModel("m5.2s", stan5_2)
rc5_2s = stan_sample(m5_2s; data)

m5_3s = SampleModel("m5.3s", stan5_3)
rc5_3s = stan_sample(m5_3s; data)
models = [m5_1s, m5_2s, m5_3s]

if success(rc5_1s) && success(rc5_2s) && success(rc5_3s)

    m5_1s_df = read_samples(m5_1s, :dataframe)
    m5_2s_df = read_samples(m5_2s, :dataframe)
    m5_3s_df = read_samples(m5_3s, :dataframe)

    if isinteractive()

        coeftab_plot(m5_1s_df, m5_2s_df, m5_3s_df; pars=(:bA, :bM),
            names=["m5.1", "m5.2", "m5.3"])
        savefig(joinpath(@__DIR__, "coeftab_plot.png"))

    end
end


text:
"""
    CustomCallback(f, TEvent, [TPhase = Phase, access = (;)])

A callback that runs `f(learner)` every time `TEvent` is triggered
in `TPhase`.

If `f` needs to access learner state, pass `access`, a named tuple
in the same form as [`stateaccess`](#).
"""
mutable struct CustomCallback{E<:Event,P<:Phase} <: Callback
    f::Any
    access
end

function CustomCallback(f, E::Type{<:Event}, P::Type{<:Phase}, access = (;))
    return CustomCallback{E, P}(f, access)
end

stateaccess(cc::CustomCallback) = cc.access

function on(::E, ::P, cb::CustomCallback{E,P}, learner) where {E<:Event,P<:Phase}
    cb.f(learner)
end


text:
module TestFlatLens

include("preamble.jl")

@testset begin
    l = FlatLens(2, 3, 4)
    @test get(((1, 2), (3, 4, 5), (6, 7, 8, 9)), l) ==
        (1, 2, 3, 4, 5, 6, 7, 8, 9)
    @test set(nothing, l, (1, 2, 3, 4, 5, 6, 7, 8, 9)) ==
        ((1, 2), (3, 4, 5), (6, 7, 8, 9))
end

end  # module


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), Apr 21, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames

# ------------------------------------------------------------------------------------------
# ## Possible pitfalls
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Know what is copied when creating a `DataFrame`
# ------------------------------------------------------------------------------------------

x = DataFrame(rand(3, 5))

y = DataFrame(x)
x === y # no copyinng performed

y = copy(x)
x === y # not the same object

all(x[i] === y[i] for i in ncol(x)) # but the columns are the same

x = 1:3; y = [1, 2, 3]; df = DataFrame(x=x,y=y) # the same when creating arrays or assigning columns, except ranges

y === df[:y] # the same object

typeof(x), typeof(df[:x]) # range is converted to a vector

# ------------------------------------------------------------------------------------------
# ### Do not modify the parent of `GroupedDataFrame`
# ------------------------------------------------------------------------------------------

x = DataFrame(id=repeat([1,2], outer=3), x=1:6)
g = groupby(x, :id)

x[1:3, 1]=[2,2,2]
g # well - it is wrong now, g is only a view

# ------------------------------------------------------------------------------------------
# ### Remember that you can filter columns of a `DataFrame` using booleans
# ------------------------------------------------------------------------------------------

srand(1)
x = DataFrame(rand(5, 5))

x[x[:x1] .< 0.25] # well - we have filtered columns not rows by accident as you can select columns using booleans

x[x[:x1] .< 0.25, :] # probably this is what we wanted

# ------------------------------------------------------------------------------------------
# ### Column selection for DataFrame creates aliases unless explicitly copied
# ------------------------------------------------------------------------------------------

x = DataFrame(a=1:3)
x[:b] = x[1] # alias
x[:c] = x[:, 1] # also alias
x[:d] = x[1][:] # copy
x[:e] = copy(x[1]) # explicit copy
display(x)
x[1,1] = 100
display(x)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# datatype and API for doing checkpointing

# Enum for checkpoint states
global const CheckpointFree = 1
global const CheckpointUsed = 2

# file name to write the checkpointer to
global const Checkpointer_fname = "checkpointer.dat"

# file name stem for solution files
global const Solution_fname = "qvec"

# file name for input_file
global const Input_fname = "input_vals_restart"

# file name step for user checkpoint data
global const CheckpointData_fname = "checkpointdata"

"""
  Every function that wants to checkpoint must implement an
  AbstractCheckpointData containing all the data that needs to be loaded
  when the code restarts.  The values in the AbstractSolutionData can be
  different for each MPI process.

  This type must contain only "julia" data, ie. objects that are managed
  by Julia and have deterministic values.  No pointers, MPI communicators etc.

  File IO is expensive, so include only as much data as needed in the
  AbstractCheckpointData.  Anything that can be recalculated should be.

  Generally, only the outermost function (ie. the time stepper or the
  nonlinear solver for steady problems) will need to checkpoint.

  There are no required fields for this type
"""
abstract type AbstractCheckpointData end

"""
  This type keeps track of which checkpoints are in use and has an API
  for loading and saving checkpoints

  The fields of the type should *never* be directly accessed.  Use the API
  instead.

  Every checkpoint must be in its own directory.

  Checkpointing has two uses: loading a previous state and restarting.
  The first one is rather easy, all that needs to happen is the solution
  variables get loaded.  Restarting is more involved because all the local
  data of the NonlinearSolver that was running needs to be stored and then
  loaded again in a new session.

  An input file called $(Input_fname) is created in the current directory
  that can be used to restart the code.

  **Fields**

   * ncheckpoints: total number of checkpoints
   * paths: absolute paths to the checkpoints (array of strings)
   * status: status of each checkpoint (used or free), (array of Ints)
   * history: list of checkpoints from most recently used to least recently
              used, unused entries set to -1

"""
mutable struct Checkpointer
  ncheckpoints::Int  # number of checkpoints
  paths::Array{String, 1}  # paths to the directories
  status::Array{Int, 1}  # is checkpoint free
  history::Array{Int, 1}  # list of checkpoints, from most recently used
                          # (first) to least recently used
                          # unused entries are set to -1
end

#------------------------------------------------------------------------------
# outer constructors

"""
  Constructs a Checkpointer with a given number of checkpoints.

  The checkpoints are put into directories called "checkpoint*", where the
  "*" is replaced by the index of the checkpoint (from 1 to ncheckpoints).
  This constructor does not error if the directories already exist, however
  any data in the checkpoints may be overwritten.

  **Inputs**

   * myrank: MPI rank of this process

   * ncheckpoints: number of checkpoints, defaults to 2.  Using fewer than
                   2 is not recommended (there will be a moment in time during
                   which the old checkpoint has been partially overwritten but
                   the new checkpoint is not complete yet)

   * prefix: added to directory names (with an underscore added in-between).
              Defaults to empty string. This is useful if there are multiple
              sets of checkpoints (and therefore multiple Checkpointer objects).

  **Outputs**

   * a Checkpointer object, fully initialized
"""
function Checkpointer(myrank::Integer, ncheckpoints::Integer=2,
                      prefix::String="")

  paths = Array{String}(ncheckpoints)
  status = Array{Int}(ncheckpoints)
  history = Array{Int}(ncheckpoints)
  fill!(history, -1)

  if prefix != ""
    prefix = string(prefix, "_")
  end

  # give initial values: all checkpoints free
  for i=1:ncheckpoints

    # create an absolute path in case the user calls cd() somewhere else in
    # the code
    paths[i] = joinpath_ascii(pwd(), "$(prefix)checkpoint$i")
    status[i] = CheckpointFree

    if !isdir(paths[i]) && myrank == 0
      mkdir(paths[i])
    end

  end

  checkpointer = Checkpointer(ncheckpoints, paths, status, history)

  # make sure the flag files don't exist (perhaps left over from a previous
  # run)

  if myrank == 0
    for i=1:ncheckpoints
      deleteFlagFile(checkpointer, i)
    end
  end

  return checkpointer
end

"""

  This constructor loads a Checkpointer object from a file.  This should
  be used for restarting.  Do *not* use the other constructor for restarting.
  It will load the most recently written checkpoint that is complete (an
  important distinction if the code was killed in the middle of writing
  a checkpoint)

  **Inputs**

   * opts: the options dictionary.

  **Outputs**

   * a Checkpointer object, fully initialized and aware of which checkpoints
     are in use and which are free

  This function only loads the Checkpointer from the checkpoint.
  See [`loadLastCheckpoint`](@ref) and [`readCheckpointData`](@ref) to
  load the rest of checkpoint data.
  The checkpoint that was loaded can be accessed via [`getLastCheckpoint`](@ref)

  **Implementation notes:**

  Uses 4 keys in the options dictionary:

   * "writing checkpoint": index of checkpoint that might be complete
   * "writing_checkpoint_path": absolute path of checkpoint
   * "most_recent_checkpoint": index of checkpoint that is definitely
                                complete, but might be older than the above
   * "most_recent_checkpoint_path": absolute path of checkpoint

  This system ensure it is possible to restart even if the code is
  killed in the middle of writing a checkpoint.

  The Checkpointer object is the same on all processes.
"""
function Checkpointer(opts::Dict, myrank::Integer)

  # figure out whether checkpoint that was in the process of being written
  # finished.
  # if not, use the checkpoint before that.
  checkpoint1_path = opts["writing_checkpoint_path"]
  checkpoint2_path = opts["most_recent_checkpoint_path"]

  if checkFlagFile(checkpoint1_path)
    if myrank == 0
      println(BSTDOUT, "Loading checkpoint most recent checkpoint")
    end
    checkpoint = opts["writing_checkpoint"]
    checkpoint_path = checkpoint1_path
  else
    @assert checkFlagFile(checkpoint2_path)
    if myrank == 0
      println(BSTDOUT, "Most recent checkpoint incomplete, loading older checkpoint")
    end
    checkpoint = opts["most_recent_checkpoint"]
    checkpoint_path = checkpoint2_path
  end

  if myrank == 0
    println(BSTDOUT, "checkpoint path = ", checkpoint_path)
  end
  fname = joinpath_ascii(checkpoint_path, Checkpointer_fname)
  f = open(fname, "r")
  checkpointer = deserialize(f)
  close(f)

  return checkpointer::Checkpointer
end

import Base.copy, Base.copy!

"""
  Recursively copies all fields to make a new Checkpointer object.
  Note that changes to one Checkpointer object will not affect the other
  (ie. the record of which checkpoints are used and which are not).
  This could easily lead to corrupting a checkpoint.
  For this reason, this function should rarely be used.
"""
function copy(chkpointer::Checkpointer)

  ncheckpoints = chkpointer.ncheckpoints
  paths = copy(chkpointer.paths)
  for i=1:ncheckpoints
    paths[i] = copy(chkpointer.paths[i])
  end
  status = copy(checkpointer.status)
  history = copy(checkpointer.history)

  return Checkpointer(ncheckpoints, paths, status, history)
end

"""
  2 argument version of copy().  See that function for details.
"""
function copy!(dest::Checkpointer, src::Checkpointer)

  dest.ncheckpoints = src.ncheckpoints
  dest.paths = copy(src.paths)
  for i=1:dest.ncheckpoints
    dest.paths[i] = copy(src.paths[i])
  end
  dest.status = copy(src.status)
  dest.history = copy(src.history)

  return nothing
end
#------------------------------------------------------------------------------
# Internal functions

"""
  Call on master process only

  Does *not* check if the checkpoint is already used (for reading the
  state back after restart, this checkpoint should already be marked as
  used)

"""
function writeCheckpointer(chkpointer::Checkpointer, checkpoint::Int)

  fname = joinpath_ascii(chkpointer.paths[checkpoint], Checkpointer_fname)
  f = open(fname, "w")
  serialize(f, chkpointer)
  close(f)

  return nothing
end

"""
  Save to a specified checkpoint.  Throw error if checkpoint is not free.
  Users should not generally call this function directly.  Instead, they
  should prefer [`saveNextFreeCheckpoint`](@ref).

  This function automatically saves eqn.q_vec and the checkpointer to
  a file.  Any additional data should be in checkpoint_data.

  Note: mesh adaptation is not compatable with checkpointing
  #TODO; add a field to the mesh to record that it has been modified

  **Inputs**

   * checkpointer: the CheckPointer
   * checkpoint: the index of the checkpoint
   * mesh: an AbstractMesh object
   * sbp: SBP operator
   * eqn: an AbstractSolutionData
   * checkpoint_data: an AbstractCheckpointData.  This is the random bag
                      of data the user needs saved.

   **Inputs/Outputs**

    * opts: options dictionary

  Implementation Notes:

  Uses options dictionary keys described by [`Checkpointer`](@ref)
  Note that the checkpoint is eagerly marked as used, before finishing writing
  the checkpoint.  Upon restart the code needs to check if this checkpoint
  is really finished.
"""
function saveCheckpoint(checkpointer::Checkpointer, checkpoint::Int,
                         mesh::AbstractMesh,
                         sbp::AbstractOperator, eqn::AbstractSolutionData,
                         opts, checkpoint_data::AbstractCheckpointData)

  @assert checkpoint > 0
  @assert checkpoint <= checkpointer.ncheckpoints

  if checkpointer.status[checkpoint] != CheckpointFree
    throw(ErrorException("Cannot save to non-free checkpoint number $checkpoint"))
  end

  # wait for all processes to catch up before overwriting whatever data
  # might be in the checkpoint
  MPI.Barrier(mesh.comm)
  if mesh.myrank == 0
    deleteFlagFile(checkpointer, checkpoint)
  end

  # record the state change now, because we have started to overwrite things
  # the options dictionary handles the case of the program getting killed
  # while the checkpoint is being written
  markCheckpointUsed(checkpointer, checkpoint)

  # record which checkpoint to restart from
  opts["writing_checkpoint"] = checkpoint
  opts["writing_checkpoint_path"] = checkpointer.paths[checkpoint]

  # write the solution to the file
  # note: the residual is not saved
  solution_fname = joinpath_ascii(checkpointer.paths[checkpoint], Solution_fname)
  writeSolutionFiles(mesh, sbp, eqn, opts, solution_fname)

  # write the user data
  writeCheckpointData(checkpointer, checkpoint, checkpoint_data, mesh.myrank)

  # write input file
  if mesh.myrank == 0
    # create a file suitable for restarting
    opts_restart = copy(opts)
    opts_restart["is_restart"] = true  # general flag that all parts of
                                       # the code will have to look for

    opts_fname = joinpath_ascii(checkpointer.paths[checkpoint], Input_fname)
    opts_fname = make_input(opts_restart, opts_fname)

    # symlink it into the current directory so we know which checkpoint to
    # load on restart
    try
      rm(Input_fname)
    end
    symlink(opts_fname, Input_fname)

    #TODO: deal with mesh adaptation
    # if only coordinates have changed, need to make sure that all the
    # solver-visible parts of the mesh object do not depend on the original
    # mesh coordinates (I think deciding elementL vs elementR for interfaces
    # does)
    # OR: write mesh coordinates to a file, load original mesh, then adapt
    # it.
    # it would be faster to have a flag in the mesh determining if it has been
    # adapted, to avoid needing to write the coordinates and read them back
    # in again

    # write checkpointer to file
    # this is needed in the case of restarting an unsteady adjoint computation

    writeCheckpointer(checkpointer, checkpoint)
  end

  # wait for all processes to finish writing
  MPI.Barrier(mesh.comm)

  # write flag file
  if mesh.myrank == 0
    writeFlagFile(checkpointer, checkpoint)
  end
  opts["most_recent_checkpoint"] = checkpoint
  opts["most_recent_checkpoint_path"] = checkpointer.paths[checkpoint]

  return nothing
end

"""
  Loads a specified checkpoint.  Only loads eqn.q_vec.  Does not load
  the checkpointer (because this is loading a checkpoint and not a restart).
  Also does not load the AbstractCheckpointData.
  See [`readCheckpointData`](@ref) for loading of the AbstractCheckpointData.
  Users should generally not call this function directly.  Users should
  prefer [`loadLastCheckpoint`](@ref) whenever possible.

  Note that loading a checkpoint does not mark it as free.  Users must
  explictly call [`freeCheckpoint`](@ref).

  **Inputs**

   * checkpointer: the Checkpointer object
   * checkpoint: the index of the checkpoint
   * mesh: an AbstractMesh object
   * sbp: the SBP operator
   * opts: the options dictionary

  **Inputs/Output**

   * eqn: an AbstractSolutionData, eqn.q_vec is overwritten
"""
function loadCheckpoint(checkpointer::Checkpointer, checkpoint::Int,
                        mesh::AbstractMesh,
                        sbp::AbstractOperator, eqn::AbstractSolutionData,
                        opts)

  @assert checkpoint > 0
  @assert checkpoint <= checkpointer.ncheckpoints

  if checkpointer.status[checkpoint] != CheckpointUsed
    throw(ErrorException("Cannot load unused checkpoint number $checkpoint"))
  end

  solution_fname = joinpath_ascii(checkpointer.paths[checkpoint], Solution_fname)
  readSolutionFiles(mesh, sbp, eqn, opts, solution_fname)

  # we can't load the AbstractCheckpointData here in a type stable way

  return nothing
end

# file name of flag file
global const FlagFile_fname = "flagfile.dat"

# value to write to flag file
global const FlagFile_val = Int64(42)

"""
  Write the flag file that signifies a checkpoint is complete
  Call on master process only!

  **Inputs**

   * checkpointer: the Checkpointer object
   * checkpoint: the checkpoint index
"""
function writeFlagFile(checkpointer::Checkpointer, checkpoint::Integer)

  fpath = joinpath_ascii(checkpointer.paths[checkpoint], FlagFile_fname)
  f = open(fpath, "w")
  write(f, FlagFile_val)
  close(f)

  return nothing
end

"""
  Returns true if the flag file exists and is consistent, returns false
  otherwise

  **Inputs**

   * checkpointer: the Checkpointer object
   * checkpoint: the checkpoint index
"""
function checkFlagFile(checkpointer::Checkpointer, checkpoint::Integer)

  fpath = checkpointer.paths[checkpoint]
  return checkFlagFile(fpath)
end

"""
  Sometimes need to check the flag file before the Checkpointer is available.
  See the other method ofr details.

  **Inputs**

   * fpath: path to the checkpoint directory
"""
function checkFlagFile(fpath::AbstractString)

  fpath = joinpath_ascii(fpath, FlagFile_fname)
  if !isfile(fpath)
    return false
  end

  f = open(fpath, "r")
  val = read(f, typeof(FlagFile_val))
  close(f)

  if val == FlagFile_val
    return true
  else
    return false
  end
end

"""
  Deletes the flag file if it exists.  Does not error if flag file does not
  exist.
  Call on master process only

  **Inputs**

   * checkpointer: the Checkpointer object
   * checkpoint: the checkpoint index
"""
function deleteFlagFile(checkpointer::Checkpointer, checkpoint::Integer)

  fpath = joinpath_ascii(checkpointer.paths[checkpoint], FlagFile_fname)
  if isfile(fpath)
    rm(fpath)
  end

  return nothing
end

"""
  Writes the AbstractCheckpointData to a file.  The file can be read by
  [`readCheckpointData`](@ref).

  **Inputs**

   * checkpoint: Checkpointer object
   * checkpoint: the checkpoint index
   * obj: the AbstractCheckpointData object
   * comm_rank: the MPI rank of this process
"""
function writeCheckpointData(chkpointer::Checkpointer, chkpoint::Integer,
                             obj::AbstractCheckpointData, comm_rank::Integer)

  fname = string(get_parallel_fname(CheckpointData_fname, comm_rank), ".dat")
  fpath = joinpath_ascii(chkpointer.paths[chkpoint], fname)

  f = open(fpath, "w")
  serialize(f, obj)
  close(f)

  return nothing
end

"""
  Marks a checkpoint as unused.

  **Inputs**

   * checkpointer: the Checkpointer object
   * checkpoint: the checkpoint index
"""
function markCheckpointUsed(chkpointer::Checkpointer, chkpoint::Integer)

  @assert chkpointer.status[chkpoint] == CheckpointFree

  chkpointer.status[chkpoint] = CheckpointUsed

  # shift everything right by one
  @assert chkpointer.history[end] == -1  # there has to be at least one unused
  last_idx = 0
  for i=1:chkpointer.ncheckpoints
    if chkpointer.history[i] == -1
      last_idx = i
      break
    end
  end

  for i=(last_idx-1):-1:1
    chkpointer.history[i+1] = chkpointer.history[i]
  end

  # add the checkpoint to the head of the line
  chkpointer.history[1] = chkpoint

  return nothing
end

#------------------------------------------------------------------------------
# API

"""
  This function reads an object of type AbstractSolutionData from a file
  and reuturns the object.  This is type unstable, so every AbstractCheckpointData
  implementation should create a constructor that calls this function and
  uses a type assertion to specify that the object must be of their concrete
  type.

  **Inputs**

   * chkpointer: a Checkpointer
   * chkpoint: the index of hte checkpoint
   * comm_rank: the MPI rank of this process

  **Outputs**

   * the [`AbstractCheckpointData`](@ref) loaded from the checkpoint
"""
function readCheckpointData(chkpointer::Checkpointer, chkpoint::Integer,
                            comm_rank::Integer)

  @assert chkpoint > 0
  @assert chkpoint <= chkpointer.ncheckpoints

  fname = string(get_parallel_fname(CheckpointData_fname, comm_rank), ".dat")
  fpath = joinpath_ascii(chkpointer.paths[chkpoint], fname)

  f = open(fpath, "r")
  obj = deserialize(f)
  close(f)

  return obj
end

"""
  Simple wrapper around [`readCheckpointData`](@ref) to loading the most
  recently saved checkpoint data

  **Inputs**

   * chkpointer: a Checkpointer
   * comm_rank: MPI rank of this process
"""
function readLastCheckpointData(chkpointer::Checkpointer, comm_rank::Integer)

  chkpoint = getLastCheckpoint(chkpointer)
  obj = readCheckpointData(chkpointer, chkpoint, comm_rank)

  return obj
end

"""
  This function saves to the next free checkpoint

  **Inputs**

   * checkpointer: the Checkpointer object
   * mesh: an AbstractMesh
   * sbp: an SBP operator
   * eqn: an AbstractSolutionData object
   * checkpoint_data: the users AbstractCheckpointData implementation

  **Inputs/Outputs**

   * opts: the options dictionary

  **Outputs**

   * checkpoint: the index of the checkpoint saved.
"""
function saveNextFreeCheckpoint(checkpointer::Checkpointer,
                                mesh::AbstractMesh,
                                sbp::AbstractOperator, eqn::AbstractSolutionData,
                                opts, checkpoint_data::AbstractCheckpointData)

  checkpoint = getNextFreeCheckpoint(checkpointer)
  saveCheckpoint(checkpointer, checkpoint, mesh, sbp, eqn, opts, checkpoint_data)

  return checkpoint

end

"""
  This function loads the most recently saved checkpoint.

  **Inputs**

   * checkpointer: the CheckPointer object
   * mesh: an AbstractMesh
   * sbp: an SBP operator
   * opts: the options dictionary

  **Inputs/Outputs**

   * eqn: an AbstractSolutionData object

  **Outputs**

   * the checkpoint index that was loaded
"""
function loadLastCheckpoint(checkpointer::Checkpointer,
                            mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts)

  checkpoint = getLastCheckpoint(checkpointer)
  loadCheckpoint(checkpointer, checkpoint, mesh, sbp, eqn, opts)

  return checkpoint
end

"""
  This function returns the number of free checkpoints

  **Inputs**

   * checkpointer: the Checkpoint object
"""
function countFreeCheckpoints(checkpointer::Checkpointer)

  cnt = 0
  for i=1:checkpointer.ncheckpoints
    if checkpointer.status[i] == CheckpointFree
      cnt += 1
    end
  end

  return cnt
end

"""
  This function returns the index of the next free checkpoint,  or 0 if
  there are no free checkpoints

  **Input**

   * checkpointer: a Checkpointer
"""
function getNextFreeCheckpoint(checkpointer::Checkpointer)

  for i=1:checkpointer.ncheckpoints
    if checkpointer.status[i] == CheckpointFree
      return i
    end
  end

  return 0  # no free checkpoints
end

"""
  This function returns the index of the most recently written checkpoint

  **Inputs**

   * checkpointer: the Checkpointer object

  **Outputs**

   * the index of the most recently saved checkpoint
"""
function getLastCheckpoint(checkpointer::Checkpointer)

  return checkpointer.history[1]
end

"""
  Get the least recently written checkpoint

  **Inputs**

   * checkpointer: the Checkpointer

  **Outputs**

   * the index of the least recently written checkpoint
"""
function getOldestCheckpoint(checkpointer::Checkpointer)

  last_idx = checkpointer.ncheckpoints
  for i=checkpointer.ncheckpoints:-1:1
    if checkpointer.history[i] != -1
      last_idx = i
      break
    end
  end

  return checkpointer.history[last_idx]
end

"""
  Frees the least recently written checkpoint. Calling this function when
  all checkpoints are free is allowed.

  **Inputs**

   * checkpointer: the Checkpointer

  **Outputs**

   * returns the checkpoint freed (0 all checkpoints were free on entry)
"""
function freeOldestCheckpoint(checkpointer::Checkpointer)

  checkpoint = getOldestCheckpoint(checkpointer)
  if checkpoint != -1
    freeCheckpoint(checkpointer, checkpoint)
  end

  return checkpoint
end

"""
  This function frees a checkpoint (marks it as available to be overwritten)
  Unlike certain algorithms that use free lists (cough, cough, malloc)
  freeing an already free checkpoint is allowed.

  **Inputs**

   * checkpointer: the Checkpointer object
   * checkpoint: the index of the checkpoint to free

  The user must explictly free checkpoints (loading a checkpoint does not
  free it).
  Note that a free checkpoint can still be loaded for restarting if it has
  not been saved to yet.
"""
function freeCheckpoint(checkpointer::Checkpointer, checkpoint::Integer)

  @assert checkpoint > 0
  @assert checkpoint <= checkpointer.ncheckpoints

  if checkpointer.status[checkpoint] == CheckpointFree
    # nothing to do
    return nothing
  end

  # remove the checkpoint from the history
  checkpoint_idx = 0
  for i=1:checkpointer.ncheckpoints
    if checkpointer.history[i] == checkpoint
      checkpoint_idx = i
      break
    end
  end

  last_idx = 0
  for i=1:checkpointer.ncheckpoints
    if checkpointer.history[i] == -1
      last_idx = i - 1
      break
    end
  end

  # shift all entries after the checkpoint to be removed left by one
  for i=checkpoint_idx:last_idx
    checkpointer.history[i] = checkpointer.history[i+1]
  end

  # if we are freeing the last checkpoint, last_idx == 0
  if last_idx != 0
    checkpointer.history[last_idx] = -1  # the last one is now unused
  else
    checkpointer.history[end] = -1
  end
  checkpointer.status[checkpoint] = CheckpointFree

  return nothing
end


text:
function align_ic(ic, sample_ids, grm_ids)
    leftjoin!(
        grm_ids,
        DataFrame(IC=ic, SAMPLE_ID=sample_ids),
        on=:SAMPLE_ID
    )
    aligned_ic = grm_ids.IC
    select!(grm_ids, Not(:IC))
    return coalesce.(aligned_ic, 0)
end

"""
    bit_distances(sample_grm, nτs)

Returns a matrix of shape (n_samples,nτs) where n_samples is the
size of sample_grm.
The sample_grm comes from the gcta software.
The process is as follows:
- Round between -1 and 1 as some elements may be beyond that limit
- Take 1 - this value to convert this quantity to a distance
- For each τ return if the distance between individuals is less or equal than τ
"""
function bit_distances(sample_grm, τs)
    distances = 1 .-  max.(min.(sample_grm, 1), -1)
    return convert(Matrix{Float32}, permutedims(distances) .<= τs)
end

default_τs(nτs;max_τ=2) = Float32[max_τ*(i-1)/(nτs-1) for i in 1:nτs]

function build_work_list(prefix, grm_ids; pval=0.05)
    dirname_, prefix_ = splitdir(prefix)
    dirname__ = dirname_ == "" ? "." : dirname_
    hdf5files = filter(
            x -> startswith(x, prefix_) && endswith(x, ".hdf5"),
            readdir(dirname__)
    )
    hdf5files = [joinpath(dirname_, x) for x in hdf5files]

    influence_curves = Vector{Float32}[]
    n_obs = Int[]
    file_tmlereport_pairs = Pair{String, String}[]
    for hdf5file in hdf5files
        jldopen(hdf5file) do io
            tmlereports = io["TMLEREPORTS"]
            for key in keys(tmlereports)
                tmlereport = tmlereports[key]
                if pvalue(ztest(tmlereport)) <= pval
                    phenotype = string(tmlereport.target_name)
                    sample_ids = parse.(Int, io["SAMPLE_IDS"][phenotype])
                    push!(influence_curves, align_ic(tmlereport.influence_curve, sample_ids, grm_ids))
                    push!(n_obs, size(sample_ids, 1))
                    push!(file_tmlereport_pairs, hdf5file => key)
                end
            end
        end
    end
    return reduce(vcat, transpose(influence_curves)), n_obs, file_tmlereport_pairs
end

"""
    normalize(variances, n_observations)

Divides the variance estimates by the effective number of observations
used for each phenotype at estimation time.
"""
normalize!(variances, n_observations) =
    variances ./= permutedims(n_observations)

"""
    aggregate_variances(influence_curves, indicator, sample)

This function computes the sum for a single index i, see also `compute_variances`.
As the GRM is symetric it is performed as :
    2 times off-diagonal elements with j < i + diagonal term
and this for all τs. Some diagonal elements are not equal to 1 while in theory they should be,
this is an artefact of the GRM. We thus always include diagonal elements in the aggregate
whatever the value of the indicator.
"""
function aggregate_variances(influence_curves, indicator, sample)
    @views begin
        D_off_diag = transpose(influence_curves[:, 1:sample-1])
        D_diag = transpose(influence_curves[:, sample])
        return D_diag .* (2indicator[:, 1:sample-1] * D_off_diag .+ D_diag)
    end
end

"""
    compute_variances(influence_curves, nτs, grm_files)

An overall variance estimate for a distance function d, a threshold τ
and influence curve D is given by:
            σ̂ = 1/n ∑ᵢ∑ⱼ1(dᵢⱼ ≤ τ)DᵢDⱼ
              = 1/n 2∑ᵢDᵢ∑ⱼ<ᵢ1(dᵢⱼ ≤ τ)DᵢDⱼ + 1(dᵢᵢ ≤ τ)DᵢDᵢ

This function computes those variance estimates at each τ for all phenotypes
and queries.

# Arguments:
- influence_curves: Array of size (n_samples, n_queries, n_phenotypes)
- grm: Vector containing the lower elements of the GRM ie: n(n+1)/2 elements
- τs: 1 row matrix containing the distance thresholds between individuals
- n_obs: Vector of size (n_phenotypes,), containing the number of effective
observations used during estimation

# Returns:
- variances: An Array of size (nτs, n_curves) where n_curves is the number of influence curves.
"""
function compute_variances(influence_curves, grm, τs, n_obs)
    n_curves, n_samples = size(influence_curves)
    variances = zeros(Float32, length(τs), n_curves)
    start_idx = 1
    for sample in 1:n_samples
        # lower diagonal of the GRM are stored in a single vector
        # that are accessed one row at a time
        end_idx = start_idx + sample - 1
        sample_grm = view(grm, start_idx:end_idx)
        indicator = bit_distances(sample_grm, τs)
        variances .+= aggregate_variances(influence_curves, indicator, sample)
        start_idx = end_idx + 1
    end
    normalize!(variances, n_obs)
    return variances
end

function grm_rows_bounds(n_samples)
    bounds = Pair{Int, Int}[]
    start_idx = 1
    for sample in 1:n_samples
        # lower diagonal of the GRM are stored in a single vector
        # that are accessed one row at a time
        end_idx = start_idx + sample - 1
        push!(bounds, start_idx => end_idx)
        start_idx = end_idx + 1
    end
    return bounds
end

function save_results(outfilename, τs, variances, std_errors, file_queryreport_pairs)
    jldopen(outfilename, "w") do io
        io["TAUS"] = τs
        io["VARIANCES"] = variances
        io["SOURCEFILE_REPORTID_PAIRS"] = file_queryreport_pairs
        io["STDERRORS"] = std_errors
    end
end

corrected_stderrors(variances, n_obs) =
    sqrt.(view(maximum(variances, dims=1), 1, :) ./ n_obs)

function sieve_variance_plateau(parsed_args)
    prefix = parsed_args["prefix"]
    pval = parsed_args["pval"]
    outfilename = parsed_args["out"]

    τs = default_τs(parsed_args["nb-estimators"];max_τ=parsed_args["max-tau"])
    grm, grm_ids = readGRM(parsed_args["grm-prefix"])

    influence_curves, n_obs, file_queryreport_pairs = build_work_list(prefix, grm_ids; pval=pval)
    variances = compute_variances(influence_curves, grm, τs, n_obs)
    std_errors = corrected_stderrors(variances, n_obs)

    save_results(outfilename, τs, variances, std_errors, file_queryreport_pairs)
end


text:
#################### Model-Based Posterior Statistics ####################
"""
    dic(mc::ModelChains)

Compute the Deviance Information Criterion (DIC) of Spiegelhalter et al. and Gelman et al. from MCMC sampler output.

Returns a `ChainSummary` type object with DIC results from the methods of Spiegelhalter and Gelman in the first and second rows of the `value` field, and the DIC value and effective numbers of parameters in the first and second columns; where

``\\text{DIC} = -2 \\mathcal{L}(\\bar{\\Theta}) + 2 p,``

such that ``\\mathcal{L}(\\bar{\\Theta})`` is the log-likelihood of model outputs given the expected values of model parameters ``\\Theta``, and ``p`` is the effective number of parameters. The latter is defined as ``p_D = -2 \\bar{\\mathcal{L}}(\\Theta) + 2 \\mathcal{L}(\\bar{\\Theta})`` for the method of Spiegelhalter and as ``p_V = \\frac{1}{2} \\operatorname{var}(-2 \\mathcal{L}(\\Theta))`` for the method of Gelman. Results are for all chains combined.

* `mc` : sampler output from a model fit with the `mcmc()` function.
"""
function dic(mc::ModelChains)
  nodekeys = keys(mc.model, :output)

  Dhat = -2.0 * logpdf(mc, mean, nodekeys)
  D = -2.0 * logpdf(mc, nodekeys).value
  p = [mean(D) - Dhat, 0.5 * var(D)]

  ChainSummary([Dhat .+ 2.0 .* p  p], ["pD", "pV"],
               ["DIC", "Effective Parameters"], header(mc))
end

"""
    logpdf(mc::ModelChains, f::Function, nodekeys::Vector{Symbol})

Compute the sum of log-densities at each iteration of MCMC output for stochastic nodes.

Returns a `ModelChains` object of resulting summed log-densities at each MCMC iteration of the supplied chain.

* `mc` : sampler output from a model fit with the `mcmc()`` function.

*  `nodekey/nodekeys` : stochastic model node(s) over which to sum densities (default: all).

* `f` : ??
"""
function logpdf(mc::ModelChains, f::Function, nodekeys::Vector{Symbol})
  m = mc.model

  relistkeys, updatekeys = getsimkeys(mc, nodekeys)
  relistkeys = union(relistkeys, intersect(nodekeys, keys(m, :block)))
  inds = names2inds(mc, relistkeys)
  for key in relistkeys
    isa(m[key], TreeStochastic) && throw("not possible with tree objects")
  end

  m[relistkeys] = relist(m, map(i -> f(mc.value[:, i, :]), inds), relistkeys)
  update!(m, updatekeys)
  mapreduce(key -> logpdf(m[key]), +, nodekeys)
end

logpdf(mc::ModelChains, nodekey::Symbol) = logpdf(mc, [nodekey])
"""
    logpdf(mc::ModelChains,
            nodekeys::Vector{Symbol}=keys(mc.model, :stochastic))

"""
function logpdf(mc::ModelChains,
                nodekeys::Vector{Symbol}=keys(mc.model, :stochastic))
  N = length(mc.range)
  K = size(mc, 3)

  relistkeys, updatekeys = getsimkeys(mc, nodekeys)
  relistkeys = union(relistkeys, intersect(nodekeys, keys(mc.model, :block)))
  inds = names2inds(mc, relistkeys)

  frame = ChainProgressFrame(
    "MCMC Processing of $N Iterations x $K Chain" * "s"^(K > 1), true
  )

  lsts = [
    Any[mc[:, :, [k]], nodekeys, relistkeys, inds, updatekeys,
        ChainProgress(frame, k, N)]
    for k in 1:K
  ]
  sims = pmap2(logpdf_modelchains_worker, lsts)

  ModelChains(cat(sims..., dims=3), mc.model)
end

function logpdf_modelchains_worker(args::Vector)
  mc, nodekeys, relistkeys, inds, updatekeys, meter = args
  m = mc.model

  sim = Chains(size(mc, 1), 1, start=first(mc), thin=step(mc), names=["logpdf"])

  for i in 1:size(mc.value, 1)
    m[relistkeys] = relist(m, mc.value[i, inds, 1], relistkeys)
    update!(m, updatekeys)
    sim.value[i, 1, 1] = mapreduce(key -> logpdf(m[key]), +, nodekeys)
    next!(meter)
  end

  sim
end

predict(mc::ModelChains, nodekey::Symbol) = predict(mc, [nodekey])
"""
    predict(mc::ModelChains,
             nodekeys::Vector{Symbol}=keys(mc.model, :output))

Generate MCMC draws from a posterior predictive distribution.

Returns a `ModelChains` object of draws simulated at each MCMC iteration of the supplied chain. For observed data node ``y``, simulation is from the posterior predictive distribution

``p(\\tilde{y} | y) = \\int p(\\tilde{y} | \\Theta) p(\\Theta | y) d\\Theta,``

where ``\\tilde{y}`` is an unknown observation on the node, ``p(\\tilde{y} | \\Theta)`` is the data likelihood, and ``p(\\Theta | y)`` is the posterior distribution of unobserved parameters ``\\Theta``.

* `mc` : sampler output from a model fit with the `mcmc()`` function.

* `nodekey/nodekeys` : observed Stochastic model node(s) for which to generate draws from the predictive distribution (default: all observed data nodes).
"""
function predict(mc::ModelChains,
                 nodekeys::Vector{Symbol}=keys(mc.model, :output))
  m = mc.model

  outputs = keys(m, :output)
  all(key -> key in outputs, nodekeys) ||
    throw(ArgumentError(string(
      "nodekeys are not all observed Stochastic nodess : ",
      join(map(string, outputs), ", ")
    )))

  nodenames = names(m, nodekeys)
  relistkeys, updatekeys = getsimkeys(mc, nodekeys)
  for key in relistkeys
    isa(m[key], TreeStochastic) && throw("not possible with tree objects")
  end

  inds = names2inds(mc, relistkeys)

  c = Chains(size(mc, 1), length(nodenames), chains=size(mc, 3),
             start=first(mc), thin=step(mc), names=nodenames)

  iters, _, chains = size(c.value)
  for k in 1:chains
    for i in 1:iters
      m[relistkeys] = relist(m, mc.value[i, inds, k], relistkeys)
      update!(m, updatekeys)
      f = key -> unlist(m[key], rand(m[key]))
      c.value[i, :, k] = vcat(map(f, nodekeys)...)
    end
  end

  ModelChains(c, m)
end

#################### Auxiliary Functions ####################

function getsimkeys(mc::ModelChains, nodekeys::Vector{Symbol})
  relistkeys = Symbol[]
  updatekeys = Symbol[]

  m = mc.model
  dag = ModelGraph(m)

  nodekeys = intersect(nodekeys, keys(m, :stochastic))
  blockkeys = keys(m, :block)
  dynamickeys = union(blockkeys, keys(m, :target, blockkeys))
  terminalkeys = union(keys(m, :stochastic), keys(mc, :dependent))

  for v in vertices(dag.graph)
    vkey = dag.keys[v]
    if vkey in dynamickeys
      if any(key -> key in nodekeys, gettargets(dag, v, terminalkeys))
        vkey in terminalkeys ?
          push!(relistkeys, vkey) :
          push!(updatekeys, vkey)
      end
    end
  end
  if !isempty(relistkeys) append!(updatekeys, nodekeys) end

  relistkeys, intersect(keys(m, :dependent), updatekeys)
end