text:
@testset "zygote_adjoints" begin
    rng = MersenneTwister(123456)
    x = rand(rng, 5)
    y = rand(rng, 5)
    r = rand(rng, 5)
    Q = Matrix(Cholesky(rand(rng, 5, 5), 'U', 0))
    @assert isposdef(Q)

    gzeucl = gradient(:Zygote, [x, y]) do xy
        evaluate(Euclidean(), xy[1], xy[2])
    end
    gzsqeucl = gradient(:Zygote, [x, y]) do xy
        evaluate(SqEuclidean(), xy[1], xy[2])
    end
    gzdotprod = gradient(:Zygote, [x, y]) do xy
        evaluate(KernelFunctions.DotProduct(), xy[1], xy[2])
    end
    gzdelta = gradient(:Zygote, [x, y]) do xy
        evaluate(KernelFunctions.Delta(), xy[1], xy[2])
    end
    gzsinus = gradient(:Zygote, [x, y]) do xy
        evaluate(KernelFunctions.Sinus(r), xy[1], xy[2])
    end
    gzsqmaha = gradient(:Zygote, [Q, x, y]) do xy
        evaluate(SqMahalanobis(xy[1]), xy[2], xy[3])
    end

    gfeucl = gradient(:FiniteDiff, [x, y]) do xy
        evaluate(Euclidean(), xy[1], xy[2])
    end
    gfsqeucl = gradient(:FiniteDiff, [x, y]) do xy
        evaluate(SqEuclidean(), xy[1], xy[2])
    end
    gfdotprod = gradient(:FiniteDiff, [x, y]) do xy
        evaluate(KernelFunctions.DotProduct(), xy[1], xy[2])
    end
    gfdelta = gradient(:FiniteDiff, [x, y]) do xy
        evaluate(KernelFunctions.Delta(), xy[1], xy[2])
    end
    gfsinus = gradient(:FiniteDiff, [x, y]) do xy
        evaluate(KernelFunctions.Sinus(r), xy[1], xy[2])
    end
    gfsqmaha = gradient(:FiniteDiff, [Q, x, y]) do xy
        evaluate(SqMahalanobis(xy[1]), xy[2], xy[3])
    end

    @test all(gzeucl .≈ gfeucl)
    @test all(gzsqeucl .≈ gfsqeucl)
    @test all(gzdotprod .≈ gfdotprod)
    @test all(gzdelta .≈ gfdelta)
    @test all(gzsinus .≈ gfsinus)
    @test all(gzsqmaha .≈ gfsqmaha)
end


text:
#function to obtain positional likelihoods for a sequence under a given BHMM.
function get_BGHMM_symbol_lh(seq::AbstractMatrix, hmm::BHMM)
    @assert size(seq)[1] == 1
    (seq=Array(transpose(seq))) # one sequence at a time only
    symbol_lhs = zeros(length(seq))
    length_mask = [length(seq)-1]

    lls = bw_llhs(hmm, seq) #obtain log likelihoods for sequences and states
    log_α = messages_forwards_log(hmm.a, hmm.A, lls, length_mask) #get forward messages
    log_β = messages_backwards_log(hmm.A, lls, length_mask) # get backwards messages

    #calculate observation probability and γ weights
    K,Tmaxplus1,Strand = size(lls) #the last T value is the 0 end marker of the longest T

    #transforms to cut down log_ξ, log_γ assignment times
    lls = permutedims(lls, [2,1,3]) # from (K,T) to (T,K)
    log_α = permutedims(log_α, [2,1,3])
    log_β = permutedims(log_β, [2,1,3])

    log_γ = fill(-Inf, Tmaxplus1,K)
    log_pobs = logsumexp(lps.(log_α[1,:], log_β[1,:]))

    @inbounds for i = 1:K, t = 1:length_mask[1]
            log_γ[t,i] = lps(log_α[t,i],log_β[t,i],-log_pobs)
    end

    for t in 1:length_mask[1]
        symbol_lh::AbstractFloat = -Inf #ie log(p=0)
        for k = 1:K #iterate over states
                state_symbol_lh::AbstractFloat = lps(log_γ[t,k], log(hmm.B[k].p[seq[t]])) #state symbol likelihood is the γ weight * the state symbol probability (log implementation)
                symbol_lh = logaddexp(symbol_lh, state_symbol_lh) #sum the probabilities over states
        end
        symbol_lhs[t] = symbol_lh
    end

    return symbol_lhs[1:end-1] #remove trailing index position
end

#function to calculate BGHMM from an observation set and a dict of BGHMMs
function BGHMM_likelihood_calc(observations::DataFrame, BGHMM_dict::Dict, code_partition_dict = get_partition_code_dict(false); symbol=:PadSeq)
    lh_matrix_size = ((findmax(length.(collect(values(observations[!, symbol]))))[1]), length(observations[!, symbol]))
    BGHMM_lh_matrix = zeros(lh_matrix_size) #T, Strand, O

    BGHMM_fragments = fragment_observations_by_BGHMM(observations[!, symbol], observations.MaskMatrix)

    #@showprogress 1 "Writing frags to matrix.."
    Threads.@threads for (jobid, frag) in BGHMM_fragments
        (frag_start, o, partition, strand) = jobid

        partition_BGHMM::BHMM = BGHMM_dict[code_partition_dict[partition]]
        no_symbols = length(partition_BGHMM.B[1].p)
        order = Int(log(4,no_symbols) - 1)

        order_seq = get_order_n_seqs([frag], order)
        coded_seq = code_seqs(order_seq)

        subseq_symbol_lh = get_BGHMM_symbol_lh(coded_seq, partition_BGHMM)

        if strand == -1
            subseq_symbol_lh = reverse(subseq_symbol_lh)
        end #positive, unstranded frags  are inserted as-is
        BGHMM_lh_matrix[frag_start:frag_start+length(frag)-1,o] = subseq_symbol_lh
    end

    return BGHMM_lh_matrix
end

function fragment_observations_by_BGHMM(seqs::AbstractVector, masks::AbstractVector)
    likelihood_jobs = Vector{Tuple{Tuple,LongSequence{DNAAlphabet{2}}}}()
    @showprogress 1 "Fragmenting observations by partition..." for (o, obs_seq) in enumerate(seqs)
        mask = masks[o]
        frags = Vector{LongSequence{DNAAlphabet{2}}}() #container for all subsequences in observation
        frag = LongSequence{DNAAlphabet{2}}()

        frag_end=0
        frag_start = 1

        while frag_start < length(obs_seq) # while we're not at the sequence end
            curr_partition = mask[frag_start,1] #get the partition code of the frag start
            curr_strand = mask[frag_start,2] #get the strand of the frag start

            #JOBID COMPOSED HERE
            jobid = (frag_start, o, curr_partition, curr_strand) #compose an identifying index for this frag

            findnext(!isequal(curr_partition),mask[:,1],frag_start) != nothing ? frag_end = findnext(!isequal(curr_partition),mask[:,1],frag_start) -1 : frag_end = length(obs_seq) #find the next position in the frag that has a different partition mask value from hte current one and set that position-1 to frag end, alternately frag end is end of the overall  sequence
            frag = obs_seq[frag_start:frag_end] #get the frag bases
            if curr_strand == -1 #if the fragment is reverse stranded
                reverse_complement!(frag) #use the reverse complement sequence
            end

            push!(likelihood_jobs,(jobid, frag)) #put the frag in the jobs vec
            frag_start = frag_end + 1 #move on
        end
    end
    return likelihood_jobs
end

ext: jl


text:
function removenode(C, A, C_age, h,r) #depends only on C operates on everything

row,col = findnz(C);

maxa = maximum(row);
for i = 1:maxa

    if isempty(row)||maximum(row)<maxa
        break
    end
    if isempty(find(row .== i))
        C = clipsimmat(C,i);
        A = clipA(A,i);
        C_age = clipsimmat(C_age,i);
        h = clipvect(h,i);
        r = r-1;
        row,col,() = findnz(C);
    end

end
return C, A, C_age, h,r

end

function clipsimmat(C,i)
   if i>size(C,2)
   error("out of bounds")
end

C = C[:, 1:size(C,2) .!= i];
C = C[1:size(C,1) .!= i,:];
ZORE = zeros(size(C,1),1);
C = [C ZORE];
ZERO = zeros(1,size(C,2));
C = [C;ZERO];
return C
end

function clipvect(V, i)
if i>size(V,2)
   error("out of bounds")
end

V = V[:,1:size(V,2) .!= i];
V = [V 0];
end

function clipA(A, i)
if i>size(A,2)
   error("out of bounds")
end
A = A[:,1:size(A,2) .!= i];
ZERO = zeros(size(A,1),1);
A = [A ZERO];
return A
end


text:
let
    P3 = Point3{Float64}
    # p = Translation(1.0,2.0,3.0)
    # r = zero(UnitQuaternion)
    # t = compose(p,LinearMap(r))
    x = P3(0.0,0.0,0.0)
    l = Line(zeros(P3),ones(P3))
    tr = Triangle(
        zeros(P3),
        ones(P3),
        0.5*ones(P3))
    q = GeometryBasics.Ngon{3,Float64,4,P3}(
        zeros(P3),
        ones(P3),
        -ones(P3),
        ones(P3).+P3(1,2,3))

    for t in [
            compose(Translation(0.0,0.0,0.0),LinearMap(zero(UnitQuaternion))),
            compose(Translation(1.0,2.0,3.0),LinearMap(zero(UnitQuaternion))),
            compose(Translation(1.0,2.0,3.0),LinearMap(RotZ(π/4))),
        ]
        for g in [l,tr,q]
            gt = t(g) # transform
            for (p,pt) in zip(g.points,gt.points)
                @test array_isapprox(t(p),pt)
            end
            n = LDrawParser.NgonElement(1,g)
            nt = t(n)
            for (p,pt) in zip(n.geom.points,nt.geom.points)
                @test array_isapprox(t(p),pt)
            end
        end
        n = LDrawParser.OptionalLineElement(1,l,l)
        nt = t(n)
        for (p,pt) in zip(n.geom.points,nt.geom.points)
            @test array_isapprox(t(p),pt)
        end
    end
end


text:
__precompile__()

module BioDataStructures

include("IntRangeSet.jl")
include("IntRangeDict.jl")

end


text:
include("traces.jl")
include("adam.jl")

export SGD, SGA
export Adam

mutable struct SGD{T, TF} <: AbstractOptimizer where {T<:Real,TF<:AbstractFuncApprox}
    α::T
    f::TF
    function SGD(fun::AbstractFuncApprox, α::T) where {T<:Real}
        tp = eltype(get_params(fun))
        new{tp, typeof(fun)}(eltype(get_params(fun))(α), fun)
    end
end

mutable struct SGA{T, TF} <: AbstractOptimizer where {T<:Real,TF<:AbstractFuncApprox}
    α::T
    f::TF
    function SGA(fun::AbstractFuncApprox, α::T) where {T<:Real}
        tp = eltype(get_params(fun))
        new{tp, typeof(fun)}(eltype(get_params(fun))(α), fun)
    end
end

function update!(opt::SGD{T}, g::Array{T}) where{T<:Real}
    add_to_params!(opt.f, -opt.α .* g)
end

function update!(opt::SGD{T}, delta::T, g::Array{T}) where{T<:Real}
    add_to_params!(opt.f, -opt.α * delta .* g)
end

function update!(opt::SGA{T}, g::Array{T}) where{T<:Real}
    add_to_params!(opt.f, opt.α .* g)
end

function update!(opt::SGA{T}, delta::T, g::Array{T}) where{T<:Real}
    add_to_params!(opt.f, opt.α * delta .* g)
end

function new_episode!(opt::AbstractOptimizer, rng=nothing)
    return nothing
end

function clone(opt::SGD)
    opt2 = SGD(opt.f, opt.α)
end

function clone(opt::SGD, f::AbstractFuncApprox)
    opt2 = SGD(f, opt.α)
end

function clone(opt::SGA)
    opt2 = SGA(opt.f, opt.α)
end

function clone(opt::SGA, f::AbstractFuncApprox)
    opt2 = SGA(f, opt.α)
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
e = exp(1.0)

function area(radius)
    a = pi * radius^2
    return a
end

A = area(1)

@show e area(1);

# More Recursion
# Turing Thesis

# Checking Types


text:
@discord_object struct Emoji
    id::Snowflake
    name::String
    roles::Vector{Snowflake}
    user::User
    require_colons::Bool
    managed::Bool
    animated::Bool
    available::Bool
end


text:
using MathProgBase
export CplexSolver

mutable struct CplexMathProgModel <: AbstractLinearQuadraticModel
    inner::Model
    lazycb
    cutcb
    heuristiccb
    branchcb
    incumbentcb
    infocb
    solvetime::Float64
    mipstart_effortlevel::Cint
    heuristic_buffer::Vector{Float64}
end

function CplexMathProgModel(;mipstart_effortlevel::Cint = CPX_MIPSTART_AUTO, 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, mipstart_effortlevel, [])
    return m
end

mutable struct CplexSolver <: AbstractMathProgSolver
    options
end
CplexSolver(;kwargs...) = CplexSolver(kwargs)
MathProgBase.LinearQuadraticModel(s::CplexSolver) = CplexMathProgModel(;s.options...)

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

function MathProgBase.setparameters!(s::CplexSolver; mpboptions...)
    opts = collect(Any,s.options)
    for (optname, optval) in mpboptions
        if optname == :TimeLimit
            push!(opts, (:CPX_PARAM_TILIM, optval))
        elseif optname == :Silent
            if optval == true
                push!(opts, (:CPX_PARAM_SCRIND, 0))
            end
        else
            error("Unrecognized parameter $optname")
        end
    end
    s.options = opts
    return
end

function MathProgBase.setparameters!(s::CplexMathProgModel; mpboptions...)
    for (optname, optval) in mpboptions
        if optname == :TimeLimit
            setparam!(m.inner, "CPX_PARAM_TILIM", optval)
        elseif optname == :Silent
            if optval == true
                setparam!(m.inner,"CPX_PARAM_SCRIND",0)
            end
        else
            error("Unrecognized parameter $optname")
        end
    end
end

function MathProgBase.loadproblem!(m::CplexMathProgModel, filename::String)
   read_model(m.inner, filename)
end

function MathProgBase.loadproblem!(m::CplexMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
  # throw away old model but keep env
  m.inner = Model(m.inner.env)
  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 = Vector{Float64}(undef, length(rowlb))
    senses = Vector{Cchar}(undef, 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

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

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

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

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

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

function MathProgBase.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

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

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

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

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

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

function MathProgBase.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 MathProgBase.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 occursin("TIME_LIM", string(ret)) || occursin("MIP_ABORT", string(ret))
        :UserLimit
    else
        ret
    end)
end

MathProgBase.getobjval(m::CplexMathProgModel)   = get_objval(m.inner)
MathProgBase.getobjbound(m::CplexMathProgModel) = get_best_bound(m.inner)
MathProgBase.getsolution(m::CplexMathProgModel) = get_solution(m.inner)
MathProgBase.getconstrsolution(m::CplexMathProgModel) = get_constr_solution(m.inner)
MathProgBase.getreducedcosts(m::CplexMathProgModel) = get_reduced_costs(m.inner)
MathProgBase.getconstrduals(m::CplexMathProgModel) = get_constr_duals(m.inner)
MathProgBase.getrawsolver(m::CplexMathProgModel) = m.inner
MathProgBase.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 MathProgBase.setvartype!(m::CplexMathProgModel, v::Vector{Symbol})
    target_int = all(x->isequal(x,:Cont), v)
    prob_type = get_prob_type(m.inner)
    if target_int
        if m.inner.has_sos # if it has sos we need to keep has_int==true and the MI(prob_type) version.
            set_vartype!(m.inner, map(x->rev_var_type_map[x], v))
        else
            m.inner.has_int = false
            if !(prob_type in [:LP,:QP,:QCP])
                toggleproblemtype!(m)
            end
        end
    else
        if prob_type in [:LP,:QP,:QCP]
            toggleproblemtype!(m)
        end
        set_vartype!(m.inner, map(x->rev_var_type_map[x], v))
    end
    return nothing
end

const prob_type_toggle_map = Dict(
    :LP    => :MILP,
    :MILP  => :LP,
    :QP    => :MIQP,
    :MIQP  => :QP,
    :QCP   => :MIQCP,
    :MIQCP => :QCP
)

function toggleproblemtype!(m::CplexMathProgModel)
    prob_type = get_prob_type(m.inner)
    set_prob_type!(m.inner, prob_type_toggle_map[prob_type])
end

function MathProgBase.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 MathProgBase.getsolvetime(m::CplexMathProgModel)
    return m.solvetime
end

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

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

function MathProgBase.setwarmstart!(m::CplexMathProgModel, v)
    # This means that warm starts are ignored if you haven't called setvartype! first
    if m.inner.has_int
        set_warm_start!(m.inner, v, m.mipstart_effortlevel)
    end
end

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

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

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

MathProgBase.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 type CplexCallbackData <: MathProgCallbackData end

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

function MathProgBase.cbgetmipsolution(d::CplexCallbackData)
    @assert d.state == :MIPSol || d.state == :MIPIncumbent
    n = num_var(d.cbdata.model)
    sol = Vector{Cdouble}(undef, n)
    stat = @cpx_ccall(getcallbacknodex, Cint, (Ptr{Cvoid},Ptr{Cvoid},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 MathProgBase.cbgetmipsolution(d::CplexCallbackData, sol::Vector{Cdouble})
    @assert d.state == :MIPSol
    stat = @cpx_ccall(getcallbacknodex, Cint, (Ptr{Cvoid},Ptr{Cvoid},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 MathProgBase.cbgetlpsolution(d::CplexCallbackData)
    @assert d.state == :MIPNode || d.state == :MIPBranch
    n = num_var(d.cbdata.model)
    sol = Vector{Cdouble}(undef, n)
    stat = @cpx_ccall(getcallbacknodex, Cint, (Ptr{Cvoid},Ptr{Cvoid},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 MathProgBase.cbgetlpsolution(d::CplexCallbackData, sol::Vector{Cdouble})
    @assert d.state == :MIPNode || d.state == :MIPIncumbent || d.state == :MIPBranch
    stat = @cpx_ccall(getcallbacknodex, Cint, (Ptr{Cvoid},Ptr{Cvoid},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),
                         (:cbgetobj,CPX_CALLBACK_INFO_BEST_INTEGER,:Cdouble),
                         (:cbgetbestbound,CPX_CALLBACK_INFO_BEST_REMAINING,:Cdouble))

    @eval begin
        function (MathProgBase.$func)(d::CplexCallbackData)
            val = Vector{$(typ)}(undef, 1)
            ret = @cpx_ccall(getcallbackinfo, Cint, (Ptr{Cvoid},Ptr{Cvoid},Cint,Cint,Ptr{Cvoid}),
                              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

for (func,param,typ) in ((:cbgetnodesleft,CPX_CALLBACK_INFO_NODES_LEFT_LONG,:Int64),
                         (:cbgetmipiterations,CPX_CALLBACK_INFO_MIP_ITERATIONS_LONG,:Int64),
                         (: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 = Vector{$(typ)}(undef, 1)
            ret = @cpx_ccall(getcallbackinfo, Cint, (Ptr{Cvoid},Ptr{Cvoid},Cint,Cint,Ptr{Cvoid}),
                              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
MathProgBase.cbgetstate(d::CplexCallbackData) = d.state

#const sensemap = Dict('=' => 'E', '<' => 'L', '>' => 'G')
function MathProgBase.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 MathProgBase.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 MathProgBase.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 MathProgBase.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 MathProgBase.cbaddsolution!(d::CplexCallbackData)
    val = unsafe_wrap(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 MathProgBase.cbsetsolutionvalue!(d::CplexCallbackData,varidx,value)
    @assert 1 <= varidx <= num_var(d.cbdata.model)
    d.heur_x[varidx] = value
    d.user_solution = true
    nothing
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

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

mutable struct 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{Cvoid}, cbdata::Ptr{Cvoid}, wherefrom::Cint, userdata::Ptr{Cvoid}, userinteraction_p::Ptr{Cint})
    model = MathProgBase.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

mutable struct CplexHeuristicCallbackData <: CplexCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    sol::Vector{Float64}
    heur_x::Vector{Float64}
    isfeas_p::Ptr{Cint}
    userinteraction_p::Ptr{Cint}
    user_solution::Bool # true if user has set any solution values
end

function MathProgBase.cbgetlpsolution(d::CplexHeuristicCallbackData)
    return copy(d.sol)
end

function MathProgBase.cbgetlpsolution(d::CplexHeuristicCallbackData, sol::Vector{Cdouble})
    copyto!(sol,d.sol)
end

function masterheuristiccallback(env::Ptr{Cvoid},
                                 cbdata::Ptr{Cvoid},
                                 wherefrom::Cint,
                                 userdata::Ptr{Cvoid},
                                 objval_p::Ptr{Cdouble},
                                 xx::Ptr{Cdouble},
                                 isfeas_p::Ptr{Cint},
                                 userinteraction_p::Ptr{Cint})
    model = MathProgBase.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 = unsafe_wrap(Array, xx, numvar(model))
        cpxcb = CplexHeuristicCallbackData(cpxrawcb, state, wherefrom, sol, model.heuristic_buffer, isfeas_p, userinteraction_p, false)
        stat = model.heuristiccb(cpxcb)
        if stat == :Exit
            terminate(model.inner)
        end
        if cpxcb.user_solution # 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
            fill!(model.heuristic_buffer, NaN)
        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{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cvoid}, Ptr{Cint}))
    stat = @cpx_ccall(setlazyconstraintcallbackfunc, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      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{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cvoid}, Ptr{Cint}))
    stat = @cpx_ccall(setusercutcallbackfunc, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      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{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}))
    stat = @cpx_ccall(setheuristiccallbackfunc, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    model.heuristic_buffer = fill(NaN, numvar(model))
    nothing
end

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

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

function masterbranchcallback(env::Ptr{Cvoid},
                              cbdata::Ptr{Cvoid},
                              wherefrom::Cint,
                              userdata::Ptr{Cvoid},
                              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 = unsafe_wrap(Array, nodebeg, convert(Cint,nodecnt))::Vector{Cint} .+ 1
        idxs = unsafe_wrap(Array, indices, sum(numbranchingvars))::Vector{Cint}
        vals = unsafe_wrap(Array, bd, sum(numbranchingvars))::Vector{Cdouble}
        dirs = unsafe_wrap(Array, lu, sum(numbranchingvars))::Vector{Cchar}
        nodes = Vector{BranchingChoice}(undef, 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{Cvoid},
                                                         Ptr{Cvoid},
                                                         Cint,
                                                         Ptr{Cvoid},
                                                         Cint,
                                                         Cint,
                                                         Cint,
                                                         Cint,
                                                         Ptr{Cint},
                                                         Ptr{Cint},
                                                         Ptr{Cchar},
                                                         Ptr{Cdouble},
                                                         Ptr{Cdouble},
                                                         Ptr{Cint}))
    stat = @cpx_ccall(setbranchcallbackfunc, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    nothing
end

mutable struct 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{Cvoid},
                                 cbdata::Ptr{Cvoid},
                                 wherefrom::Cint,
                                 userdata::Ptr{Cvoid},
                                 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 = unsafe_wrap(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{Cvoid},
                                                            Ptr{Cvoid},
                                                            Cint,
                                                            Ptr{Cvoid},
                                                            Cdouble,
                                                            Ptr{Cdouble},
                                                            Ptr{Cint},
                                                            Ptr{Cint}))
    stat = @cpx_ccall(setincumbentcallbackfunc, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Any,
                      ),
                      model.inner.env.ptr, cpxcallback, model)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    nothing
end

mutable struct CplexInfoCallbackData <: CplexCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
end

function masterinfocallback(env::Ptr{Cvoid},
                            cbdata::Ptr{Cvoid},
                            wherefrom::Cint,
                            userdata::Ptr{Cvoid})
    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{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cvoid}))
    stat = @cpx_ccall(setinfocallbackfunc, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      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 = Vector{Cdouble}(undef, n)
    stat = @cpx_ccall(getcallbacknodelb, Cint, (Ptr{Cvoid},Ptr{Cvoid},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 = Vector{Cdouble}(undef, n)
    stat = @cpx_ccall(getcallbacknodeub, Cint, (Ptr{Cvoid},Ptr{Cvoid},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 = Vector{Cdouble}(undef, 1)
    stat = @cpx_ccall(getcallbacknodeobjval, Cint, (Ptr{Cvoid},Ptr{Cvoid},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 = Vector{Cint}(undef, n)
    stat = @cpx_ccall(getcallbacknodeintfeas, Cint, (Ptr{Cvoid},Ptr{Cvoid},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:
# Compatability for Julia below 1.3.
if VERSION < v"1.3"
    allocatedinline(::Type{T}) where {T} = (Base.@_pure_meta; ccall(:jl_array_store_unboxed, Cint, (Any,), T) != Cint(0))
    function aligned_sizeof(T)
        Base.@_pure_meta
        if Base.isbitsunion(T)
            sz = Ref{Csize_t}(0)
            algn = Ref{Csize_t}(0)
            ccall(:jl_islayout_inline, Cint, (Any, Ptr{Csize_t}, Ptr{Csize_t}), T, sz, algn)
            al = algn[]
            return (sz[] + al - 1) & -al
        elseif allocatedinline(T)
            al = Base.datatype_alignment(T)
            return (Core.sizeof(T) + al - 1) & -al
        else
            return Core.sizeof(Ptr{Cvoid})
        end
    end
else
    const allocatedinline = Base.allocatedinline
    const aligned_sizeof = Base.aligned_sizeof
end

size_to_fstrides(elsz::Integer, sz::Integer...) =
    isempty(sz) ? () : (elsz, size_to_fstrides(elsz * sz[1], sz[2:end]...)...)

size_to_cstrides(elsz::Integer, sz::Integer...) =
    isempty(sz) ? () : (size_to_cstrides(elsz * sz[end], sz[1:end-1]...)..., elsz)

isfcontiguous(o::AbstractArray) = strides(o) == size_to_fstrides(1, size(o)...)
isccontiguous(o::AbstractArray) = strides(o) == size_to_cstrides(1, size(o)...)

# TODO: what is the best way?
ismutablearray(x::Array) = true
ismutablearray(x::AbstractArray) = begin
    p = parent(x)
    p === x ? false : ismutablearray(p)
end

pybufferformat(::Type{T}) where {T} =
    T == Int8 ? "=b" :
    T == UInt8 ? "=B" :
    T == Int16 ? "=h" :
    T == UInt16 ? "=H" :
    T == Int32 ? "=i" :
    T == UInt32 ? "=I" :
    T == Int64 ? "=q" :
    T == UInt64 ? "=Q" :
    T == Float16 ? "=e" :
    T == Float32 ? "=f" :
    T == Float64 ? "=d" :
    T == Complex{Float16} ? "=Ze" :
    T == Complex{Float32} ? "=Zf" :
    T == Complex{Float64} ? "=Zd" :
    T == Bool ? "?" :
    T == Ptr{Cvoid} ? "P" :
    T == C.PyObjectRef ? "O" :
    if isstructtype(T) && isconcretetype(T) && allocatedinline(T)
        n = fieldcount(T)
        flds = []
        for i = 1:n
            nm = fieldname(T, i)
            tp = fieldtype(T, i)
            push!(flds, string(pybufferformat(tp), nm isa Symbol ? ":$nm:" : ""))
            d =
                (i == n ? sizeof(T) : fieldoffset(T, i + 1)) -
                (fieldoffset(T, i) + sizeof(tp))
            @assert d ≥ 0
            d > 0 && push!(flds, "$(d)x")
        end
        string("T{", join(flds, " "), "}")
    else
        "$(sizeof(T))x"
    end

pybufferformat_to_type(fmt::AbstractString) =
    fmt == "b" ? Cchar :
    fmt == "B" ? Cuchar :
    fmt == "h" ? Cshort :
    fmt == "H" ? Cushort :
    fmt == "i" ? Cint :
    fmt == "I" ? Cuint :
    fmt == "l" ? Clong :
    fmt == "L" ? Culong :
    fmt == "q" ? Clonglong :
    fmt == "Q" ? Culonglong :
    fmt == "e" ? Float16 :
    fmt == "f" ? Cfloat :
    fmt == "d" ? Cdouble :
    fmt == "?" ? Bool :
    fmt == "P" ? Ptr{Cvoid} :
    fmt == "O" ? C.PyObjectRef :
    fmt == "=e" ? Float16 :
    fmt == "=f" ? Float32 : fmt == "=d" ? Float64 : error("not implemented: $(repr(fmt))")

islittleendian() =
    Base.ENDIAN_BOM == 0x04030201 ? true : Base.ENDIAN_BOM == 0x01020304 ? false : error()

pytypestrdescr(::Type{T}) where {T} = begin
    c = islittleendian() ? '<' : '>'
    T == Bool ? ("$(c)b$(sizeof(Bool))", nothing) :
    T == Int8 ? ("$(c)i1", nothing) :
    T == UInt8 ? ("$(c)u1", nothing) :
    T == Int16 ? ("$(c)i2", nothing) :
    T == UInt16 ? ("$(c)u2", nothing) :
    T == Int32 ? ("$(c)i4", nothing) :
    T == UInt32 ? ("$(c)u4", nothing) :
    T == Int64 ? ("$(c)i8", nothing) :
    T == UInt64 ? ("$(c)u8", nothing) :
    T == Float16 ? ("$(c)f2", nothing) :
    T == Float32 ? ("$(c)f4", nothing) :
    T == Float64 ? ("$(c)f8", nothing) :
    T == Complex{Float16} ? ("$(c)c4", nothing) :
    T == Complex{Float32} ? ("$(c)c8", nothing) :
    T == Complex{Float64} ? ("$(c)c16", nothing) :
    T == C.PyObjectRef ? ("|O", nothing) :
    if isstructtype(T) && isconcretetype(T) && allocatedinline(T)
        n = fieldcount(T)
        flds = []
        for i = 1:n
            nm = fieldname(T, i)
            tp = fieldtype(T, i)
            ts, ds = pytypestrdescr(tp)
            isempty(ts) && return ("", nothing)
            push!(
                flds,
                (nm isa Integer ? "f$(nm-1)" : string(nm), ds === nothing ? ts : ds),
            )
            d =
                (i == n ? sizeof(T) : fieldoffset(T, i + 1)) -
                (fieldoffset(T, i) + sizeof(tp))
            @assert d ≥ 0
            d > 0 && push!(flds, ("", "|V$(d)"))
        end
        ("|$(sizeof(T))V", flds)
    else
        ("", nothing)
    end
end

pytypestrdescr_to_type(ts::String, descr) = begin
    # byte swapped?
    bsc = ts[1]
    bs =
        bsc == '<' ? !islittleendian() :
        bsc == '>' ? islittleendian() :
        bsc == '|' ? false : error("endianness character not supported: $ts")
    bs && error("byte-swapping not supported: $ts")
    # element type
    etc = ts[2]
    if etc == 'b'
        sz = parse(Int, ts[3:end])
        sz == sizeof(Bool) && return Bool
        error("bool of this size not supported: $ts")
    elseif etc == 'i'
        sz = parse(Int, ts[3:end])
        sz == 1 && return Int8
        sz == 2 && return Int16
        sz == 4 && return Int32
        sz == 8 && return Int64
        sz == 16 && return Int128
        error("signed int of this size not supported: $ts")
    elseif etc == 'u'
        sz = parse(Int, ts[3:end])
        sz == 1 && return UInt8
        sz == 2 && return UInt16
        sz == 4 && return UInt32
        sz == 8 && return UInt64
        sz == 16 && return UInt128
        error("unsigned int of this size not supported: $ts")
    elseif etc == 'f'
        sz = parse(Int, ts[3:end])
        sz == 2 && return Float16
        sz == 4 && return Float32
        sz == 8 && return Float64
        error("float of this size not supported: $ts")
    elseif etc == 'c'
        sz = parse(Int, ts[3:end])
        sz == 4 && return Complex{Float16}
        sz == 8 && return Complex{Float32}
        sz == 16 && return Complex{Float64}
        error("complex of this size not supported: $ts")
    elseif etc == 'O'
        return C.PyObjectRef
    else
        error("type not supported: $ts")
    end
end

### TYPE UTILITIES

# Used to signal a Python error from functions that return general Julia objects
struct PYERR end

# Used to signal that conversion failed
struct NOTIMPLEMENTED end

# Somewhere to stash results
const RESULT = Ref{Any}(nothing)
putresult(x) = (RESULT[] = x; 1)
putresult(x::PYERR) = -1
putresult(x::NOTIMPLEMENTED) = 0
takeresult(::Type{T} = Any) where {T} = (r = RESULT[]::T; RESULT[] = nothing; r)

tryconvert(::Type{T}, x::PYERR) where {T} = PYERR()
tryconvert(::Type{T}, x::NOTIMPLEMENTED) where {T} = NOTIMPLEMENTED()
tryconvert(::Type{T}, x::T) where {T} = x
tryconvert(::Type{T}, x) where {T} =
    try
        convert(T, x)
    catch
        NOTIMPLEMENTED()
    end

CTryConvertRule_wrapref(o, ::Type{S}) where {S} = putresult(S(C.PyObjectRef(o)))
CTryConvertRule_trywrapref(o, ::Type{S}) where {S} =
    try
        putresult(S(C.PyObjectRef(o)))
    catch
        0
    end

@generated _typeintersect(::Type{T1}, ::Type{T2}) where {T1,T2} = typeintersect(T1, T2)

@generated function _type_flatten_tuple(::Type{T}) where {T<:Tuple}
    S = T
    vars = []
    while !isa(S, DataType)
        push!(vars, S.var)
        S = S.body
    end
    Tuple{[foldr(UnionAll, vars; init = P) for P in S.parameters]...}
end

function pointer_from_obj(o::T) where {T}
    if T.mutable
        c = o
        p = Base.pointer_from_objref(o)
    else
        c = Ref{Any}(o)
        p = unsafe_load(Ptr{Ptr{Cvoid}}(Base.pointer_from_objref(c)))
    end
    p, c
end

isnull(p::Ptr) = Ptr{Cvoid}(p) == C_NULL
isnull(p::UnsafePtr) = isnull(pointer(p))
ism1(x::T) where {T<:Number} = x == (zero(T) - one(T))

# Put something in here to keep it around forever
const CACHE = []


text:
import Printf: ini_dec, fix_dec, ini_hex, ini_HEX

if VERSION < v"1.1"
  fix_dec(out, d::BFloat16,               flags::String, width::Int, precision::Int, c::Char) = fix_dec(out, Float32(d),          flags, width, precision, c)
  ini_dec(out, d::BFloat16, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = ini_dec(out, Float32(d), ndigits, flags, width, precision, c)
  ini_hex(out, d::BFloat16, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = ini_hex(out, Float32(d), ndigits, flags, width, precision, c)
  ini_HEX(out, d::BFloat16, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = ini_HEX(out, Float32(d), ndigits, flags, width, precision, c)
  ini_hex(out, d::BFloat16,               flags::String, width::Int, precision::Int, c::Char) = ini_hex(out, Float32(d),          flags, width, precision, c)
  ini_HEX(out, d::BFloat16,               flags::String, width::Int, precision::Int, c::Char) = ini_HEX(out, Float32(d),          flags, width, precision, c)
else
  fix_dec(out, d::BFloat16,               flags::String, width::Int, precision::Int, c::Char, digits) = fix_dec(out, Float32(d),          flags, width, precision, c, digits)
  ini_dec(out, d::BFloat16, ndigits::Int, flags::String, width::Int, precision::Int, c::Char, digits) = ini_dec(out, Float32(d), ndigits, flags, width, precision, c, digits)
  ini_hex(out, d::BFloat16, ndigits::Int, flags::String, width::Int, precision::Int, c::Char, digits) = ini_hex(out, Float32(d), ndigits, flags, width, precision, c, digits)
  ini_HEX(out, d::BFloat16, ndigits::Int, flags::String, width::Int, precision::Int, c::Char, digits) = ini_HEX(out, Float32(d), ndigits, flags, width, precision, c, digits)
  ini_hex(out, d::BFloat16,               flags::String, width::Int, precision::Int, c::Char, digits) = ini_hex(out, Float32(d),          flags, width, precision, c, digits)
  ini_HEX(out, d::BFloat16,               flags::String, width::Int, precision::Int, c::Char, digits) = ini_HEX(out, Float32(d),          flags, width, precision, c, digits)
end


text:
# Julia wrapper for header: cone.h
# Automatically generated using Clang.jl

function bring2cone(C, r, s)
    return ccall(
        (:bring2cone, ecos),
        Cvoid,
        (Ptr{cone}, Ptr{pfloat}, Ptr{pfloat}),
        C,
        r,
        s,
    )
end

function unitInitialization(C, s, z, scaling)
    return ccall(
        (:unitInitialization, ecos),
        Cvoid,
        (Ptr{cone}, Ptr{pfloat}, Ptr{pfloat}, pfloat),
        C,
        s,
        z,
        scaling,
    )
end

function updateScalings(C, s, z, lambda, mu)
    return ccall(
        (:updateScalings, ecos),
        idxint,
        (Ptr{cone}, Ptr{pfloat}, Ptr{pfloat}, Ptr{pfloat}, pfloat),
        C,
        s,
        z,
        lambda,
        mu,
    )
end

function evalSymmetricBarrierValue(siter, ziter, tauIter, kapIter, C, D)
    return ccall(
        (:evalSymmetricBarrierValue, ecos),
        pfloat,
        (Ptr{pfloat}, Ptr{pfloat}, pfloat, pfloat, Ptr{cone}, pfloat),
        siter,
        ziter,
        tauIter,
        kapIter,
        C,
        D,
    )
end

function scale(z, C, lambda)
    return ccall(
        (:scale, ecos),
        Cvoid,
        (Ptr{pfloat}, Ptr{cone}, Ptr{pfloat}),
        z,
        C,
        lambda,
    )
end

function scale2add(x, y, C)
    return ccall(
        (:scale2add, ecos),
        Cvoid,
        (Ptr{pfloat}, Ptr{pfloat}, Ptr{cone}),
        x,
        y,
        C,
    )
end

function unscale(lambda, C, z)
    return ccall(
        (:unscale, ecos),
        Cvoid,
        (Ptr{pfloat}, Ptr{cone}, Ptr{pfloat}),
        lambda,
        C,
        z,
    )
end

function conicProduct(u, v, C, w)
    return ccall(
        (:conicProduct, ecos),
        pfloat,
        (Ptr{pfloat}, Ptr{pfloat}, Ptr{cone}, Ptr{pfloat}),
        u,
        v,
        C,
        w,
    )
end

function conicDivision(u, v, C, w)
    return ccall(
        (:conicDivision, ecos),
        Cvoid,
        (Ptr{pfloat}, Ptr{pfloat}, Ptr{cone}, Ptr{pfloat}),
        u,
        v,
        C,
        w,
    )
end

function getSOCDetails(soc, conesize, eta_square, d1, u0, u1, v1, q)
    return ccall(
        (:getSOCDetails, ecos),
        Cvoid,
        (
            Ptr{socone},
            Ptr{idxint},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{Ptr{pfloat}},
        ),
        soc,
        conesize,
        eta_square,
        d1,
        u0,
        u1,
        v1,
        q,
    )
end

function unstretch(n, p, C, Pinv, Px, dx, dy, dz)
    return ccall(
        (:unstretch, ecos),
        Cvoid,
        (
            idxint,
            idxint,
            Ptr{cone},
            Ptr{idxint},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
        ),
        n,
        p,
        C,
        Pinv,
        Px,
        dx,
        dy,
        dz,
    )
end
# Julia wrapper for header: ctrlc.h
# Automatically generated using Clang.jl

# Julia wrapper for header: data.h
# Automatically generated using Clang.jl

# Julia wrapper for header: ecos.h
# Automatically generated using Clang.jl

function ECOS_setup(
    n,
    m,
    p,
    l,
    ncones,
    q,
    nex,
    Gpr,
    Gjc,
    Gir,
    Apr,
    Ajc,
    Air,
    c,
    h,
    b,
)
    return ccall(
        (:ECOS_setup, ecos),
        Ptr{pwork},
        (
            idxint,
            idxint,
            idxint,
            idxint,
            idxint,
            Ptr{idxint},
            idxint,
            Ptr{pfloat},
            Ptr{idxint},
            Ptr{idxint},
            Ptr{pfloat},
            Ptr{idxint},
            Ptr{idxint},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
        ),
        n,
        m,
        p,
        l,
        ncones,
        q,
        nex,
        Gpr,
        Gjc,
        Gir,
        Apr,
        Ajc,
        Air,
        c,
        h,
        b,
    )
end

function expConeLineSearch(w, dtau, dkappa, affine)
    return ccall(
        (:expConeLineSearch, ecos),
        pfloat,
        (Ptr{pwork}, pfloat, pfloat, idxint),
        w,
        dtau,
        dkappa,
        affine,
    )
end

function ECOS_solve(w)
    return ccall((:ECOS_solve, ecos), idxint, (Ptr{pwork},), w)
end

function ECOS_cleanup(w, keepvars)
    return ccall(
        (:ECOS_cleanup, ecos),
        Cvoid,
        (Ptr{pwork}, idxint),
        w,
        keepvars,
    )
end

function ECOS_ver()
    return ccall((:ECOS_ver, ecos), Cstring, ())
end

function ecos_updateDataEntry_h(w, idx, value)
    return ccall(
        (:ecos_updateDataEntry_h, ecos),
        Cvoid,
        (Ptr{pwork}, idxint, pfloat),
        w,
        idx,
        value,
    )
end

function ecos_updateDataEntry_c(w, idx, value)
    return ccall(
        (:ecos_updateDataEntry_c, ecos),
        Cvoid,
        (Ptr{pwork}, idxint, pfloat),
        w,
        idx,
        value,
    )
end

function ECOS_updateData(w, Gpr, Apr, c, h, b)
    return ccall(
        (:ECOS_updateData, ecos),
        Cvoid,
        (
            Ptr{pwork},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
        ),
        w,
        Gpr,
        Apr,
        c,
        h,
        b,
    )
end
# Julia wrapper for header: ecos_bb.h
# Automatically generated using Clang.jl

function ECOS_BB_setup(
    n,
    m,
    p,
    l,
    ncones,
    q,
    nex,
    Gpr,
    Gjc,
    Gir,
    Apr,
    Ajc,
    Air,
    c,
    h,
    b,
    num_bool_vars,
    bool_vars_idx,
    num_int_vars,
    int_vars_idx,
    stgs,
)
    return ccall(
        (:ECOS_BB_setup, ecos),
        Ptr{ecos_bb_pwork},
        (
            idxint,
            idxint,
            idxint,
            idxint,
            idxint,
            Ptr{idxint},
            idxint,
            Ptr{pfloat},
            Ptr{idxint},
            Ptr{idxint},
            Ptr{pfloat},
            Ptr{idxint},
            Ptr{idxint},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            idxint,
            Ptr{idxint},
            idxint,
            Ptr{idxint},
            Ptr{settings_bb},
        ),
        n,
        m,
        p,
        l,
        ncones,
        q,
        nex,
        Gpr,
        Gjc,
        Gir,
        Apr,
        Ajc,
        Air,
        c,
        h,
        b,
        num_bool_vars,
        bool_vars_idx,
        num_int_vars,
        int_vars_idx,
        stgs,
    )
end

function ECOS_BB_solve(prob)
    return ccall((:ECOS_BB_solve, ecos), idxint, (Ptr{ecos_bb_pwork},), prob)
end

function ECOS_BB_cleanup(prob, num_vars_keep)
    return ccall(
        (:ECOS_BB_cleanup, ecos),
        Cvoid,
        (Ptr{ecos_bb_pwork}, idxint),
        prob,
        num_vars_keep,
    )
end

function updateDataEntry_h(w, idx, value)
    return ccall(
        (:updateDataEntry_h, ecos),
        Cvoid,
        (Ptr{ecos_bb_pwork}, idxint, pfloat),
        w,
        idx,
        value,
    )
end

function updateDataEntry_c(w, idx, value)
    return ccall(
        (:updateDataEntry_c, ecos),
        Cvoid,
        (Ptr{ecos_bb_pwork}, idxint, pfloat),
        w,
        idx,
        value,
    )
end

function get_default_ECOS_BB_settings()
    return ccall((:get_default_ECOS_BB_settings, ecos), Ptr{settings_bb}, ())
end

function get_bool_node_id(idx, prob)
    return ccall(
        (:get_bool_node_id, ecos),
        Cstring,
        (idxint, Ptr{ecos_bb_pwork}),
        idx,
        prob,
    )
end

function get_int_node_id(idx, prob)
    return ccall(
        (:get_int_node_id, ecos),
        Ptr{pfloat},
        (idxint, Ptr{ecos_bb_pwork}),
        idx,
        prob,
    )
end

function abs_2(number)
    return ccall((:abs_2, ecos), pfloat, (pfloat,), number)
end

function pfloat_round(number)
    return ccall((:pfloat_round, ecos), pfloat, (pfloat,), number)
end

function pfloat_ceil(number, integer_tol)
    return ccall(
        (:pfloat_ceil, ecos),
        pfloat,
        (pfloat, pfloat),
        number,
        integer_tol,
    )
end

function pfloat_floor(number, integer_tol)
    return ccall(
        (:pfloat_floor, ecos),
        pfloat,
        (pfloat, pfloat),
        number,
        integer_tol,
    )
end

function float_eqls(a, b, integer_tol)
    return ccall(
        (:float_eqls, ecos),
        idxint,
        (pfloat, pfloat, pfloat),
        a,
        b,
        integer_tol,
    )
end
# Julia wrapper for header: equil.h
# Automatically generated using Clang.jl

function set_equilibration(w)
    return ccall((:set_equilibration, ecos), Cvoid, (Ptr{pwork},), w)
end

function unset_equilibration(w)
    return ccall((:unset_equilibration, ecos), Cvoid, (Ptr{pwork},), w)
end
# Julia wrapper for header: expcone.h
# Automatically generated using Clang.jl

function evalExpHessian(w, v, mu)
    return ccall(
        (:evalExpHessian, ecos),
        Cvoid,
        (Ptr{pfloat}, Ptr{pfloat}, pfloat),
        w,
        v,
        mu,
    )
end

function evalExpGradient(w, g)
    return ccall(
        (:evalExpGradient, ecos),
        Cvoid,
        (Ptr{pfloat}, Ptr{pfloat}),
        w,
        g,
    )
end

function evalBarrierValue(siter, ziter, fc, nexc)
    return ccall(
        (:evalBarrierValue, ecos),
        pfloat,
        (Ptr{pfloat}, Ptr{pfloat}, idxint, idxint),
        siter,
        ziter,
        fc,
        nexc,
    )
end

function scaleToAddExpcone(y, x, expcones, nexc, fc)
    return ccall(
        (:scaleToAddExpcone, ecos),
        Cvoid,
        (Ptr{pfloat}, Ptr{pfloat}, Ptr{expcone}, idxint, idxint),
        y,
        x,
        expcones,
        nexc,
        fc,
    )
end

function evalExpPrimalFeas(s, nexc)
    return ccall(
        (:evalExpPrimalFeas, ecos),
        idxint,
        (Ptr{pfloat}, idxint),
        s,
        nexc,
    )
end

function evalExpDualFeas(s, nexc)
    return ccall(
        (:evalExpDualFeas, ecos),
        idxint,
        (Ptr{pfloat}, idxint),
        s,
        nexc,
    )
end
# Julia wrapper for header: glblopts.h
# Automatically generated using Clang.jl

# Julia wrapper for header: kkt.h
# Automatically generated using Clang.jl

function kkt_factor(KKT, eps, delta)
    return ccall(
        (:kkt_factor, ecos),
        idxint,
        (Ptr{kkt}, pfloat, pfloat),
        KKT,
        eps,
        delta,
    )
end

function kkt_solve(KKT, A, G, Pb, dx, dy, dz, n, p, m, C, isinit, nitref)
    return ccall(
        (:kkt_solve, ecos),
        idxint,
        (
            Ptr{kkt},
            Ptr{spmat},
            Ptr{spmat},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            Ptr{pfloat},
            idxint,
            idxint,
            idxint,
            Ptr{cone},
            idxint,
            idxint,
        ),
        KKT,
        A,
        G,
        Pb,
        dx,
        dy,
        dz,
        n,
        p,
        m,
        C,
        isinit,
        nitref,
    )
end

function kkt_update(PKP, P, C)
    return ccall(
        (:kkt_update, ecos),
        Cvoid,
        (Ptr{spmat}, Ptr{idxint}, Ptr{cone}),
        PKP,
        P,
        C,
    )
end

function kkt_init(PKP, P, C)
    return ccall(
        (:kkt_init, ecos),
        Cvoid,
        (Ptr{spmat}, Ptr{idxint}, Ptr{cone}),
        PKP,
        P,
        C,
    )
end
# Julia wrapper for header: spla.h
# Automatically generated using Clang.jl

function sparseMV(A, x, y, a, newVector)
    return ccall(
        (:sparseMV, ecos),
        Cvoid,
        (Ptr{spmat}, Ptr{pfloat}, Ptr{pfloat}, idxint, idxint),
        A,
        x,
        y,
        a,
        newVector,
    )
end

function sparseMtVm(A, x, y, newVector, skipDiagonal)
    return ccall(
        (:sparseMtVm, ecos),
        Cvoid,
        (Ptr{spmat}, Ptr{pfloat}, Ptr{pfloat}, idxint, idxint),
        A,
        x,
        y,
        newVector,
        skipDiagonal,
    )
end

function vadd(n, x, y)
    return ccall(
        (:vadd, ecos),
        Cvoid,
        (idxint, Ptr{pfloat}, Ptr{pfloat}),
        n,
        x,
        y,
    )
end

function vsubscale(n, a, x, y)
    return ccall(
        (:vsubscale, ecos),
        Cvoid,
        (idxint, pfloat, Ptr{pfloat}, Ptr{pfloat}),
        n,
        a,
        x,
        y,
    )
end

function norm2(v, n)
    return ccall((:norm2, ecos), pfloat, (Ptr{pfloat}, idxint), v, n)
end

function norminf(v, n)
    return ccall((:norminf, ecos), pfloat, (Ptr{pfloat}, idxint), v, n)
end

function eddot(n, x, y)
    return ccall(
        (:eddot, ecos),
        pfloat,
        (idxint, Ptr{pfloat}, Ptr{pfloat}),
        n,
        x,
        y,
    )
end
# Julia wrapper for header: splamm.h
# Automatically generated using Clang.jl

function ecoscreateSparseMatrix(m, n, nnz, jc, ir, pr)
    return ccall(
        (:ecoscreateSparseMatrix, ecos),
        Ptr{spmat},
        (idxint, idxint, idxint, Ptr{idxint}, Ptr{idxint}, Ptr{pfloat}),
        m,
        n,
        nnz,
        jc,
        ir,
        pr,
    )
end

function newSparseMatrix(m, n, nnz)
    return ccall(
        (:newSparseMatrix, ecos),
        Ptr{spmat},
        (idxint, idxint, idxint),
        m,
        n,
        nnz,
    )
end

function freeSparseMatrix(M)
    return ccall((:freeSparseMatrix, ecos), Cvoid, (Ptr{spmat},), M)
end

function transposeSparseMatrix(M, MtoMt)
    return ccall(
        (:transposeSparseMatrix, ecos),
        Ptr{spmat},
        (Ptr{spmat}, Ptr{idxint}),
        M,
        MtoMt,
    )
end

function permuteSparseSymmetricMatrix(A, pinv, C, PK)
    return ccall(
        (:permuteSparseSymmetricMatrix, ecos),
        Cvoid,
        (Ptr{spmat}, Ptr{idxint}, Ptr{spmat}, Ptr{idxint}),
        A,
        pinv,
        C,
        PK,
    )
end

function pinv(n, p, pinv)
    return ccall(
        (:pinv, ecos),
        Cvoid,
        (idxint, Ptr{idxint}, Ptr{idxint}),
        n,
        p,
        pinv,
    )
end

function copySparseMatrix(A)
    return ccall((:copySparseMatrix, ecos), Ptr{spmat}, (Ptr{spmat},), A)
end

function printDenseMatrix(M, dim1, dim2, name)
    return ccall(
        (:printDenseMatrix, ecos),
        Cvoid,
        (Ptr{pfloat}, idxint, idxint, Cstring),
        M,
        dim1,
        dim2,
        name,
    )
end

function printDenseMatrix_i(M, dim1, dim2, name)
    return ccall(
        (:printDenseMatrix_i, ecos),
        Cvoid,
        (Ptr{idxint}, idxint, idxint, Cstring),
        M,
        dim1,
        dim2,
        name,
    )
end

function printSparseMatrix(M)
    return ccall((:printSparseMatrix, ecos), Cvoid, (Ptr{spmat},), M)
end

function dumpSparseMatrix(M, fn)
    return ccall((:dumpSparseMatrix, ecos), Cvoid, (Ptr{spmat}, Cstring), M, fn)
end

function dumpDenseMatrix(M, dim1, dim2, fn)
    return ccall(
        (:dumpDenseMatrix, ecos),
        Cvoid,
        (Ptr{pfloat}, Cint, Cint, Cstring),
        M,
        dim1,
        dim2,
        fn,
    )
end

function dumpDenseMatrix_i(M, dim1, dim2, fn)
    return ccall(
        (:dumpDenseMatrix_i, ecos),
        Cvoid,
        (Ptr{idxint}, Cint, Cint, Cstring),
        M,
        dim1,
        dim2,
        fn,
    )
end
# Julia wrapper for header: wright_omega.h
# Automatically generated using Clang.jl

function wrightOmega(z)
    return ccall((:wrightOmega, ecos), pfloat, (pfloat,), z)
end


text:
function split_cmd(s::String)
	v=[];t="";quo=false
	l=length(s)
	for i in 1:l
		if s[i]==' '&&!quo
			if t!=""
				push!(v,t)
				t=""
			end
			continue
		end
		if s[i]=='\\'
			if s[i+1]=='\\'
				i+=1
				t*='\\'
			elseif s[i+1]=='\`'
				i+=1
				t*='\`'
			end
		elseif s[i]=='`'
			quo=!quo
		else
			t*=s[i]
		end
	end
	if t!=""
		push!(v,t)
	end
	return v
end
function clean_str(s::String)
	a=""
	for i in s
		if i!='\t'&&i!='\n'&&i!='\r'
			a*=i
		end
	end
	return a
end
function esc_str2(io::IO,s::AbstractString)
	for c::AbstractChar in s
		c=='"' ? print(io,"\\=") :
		c=='\'' ? print(io,"\\-") :
		c=='\0' ? print(io,"\\0") :
		c=='\e' ? print(io,"\\e") :
		c=='\\' ? print(io,"\\/") :
		'\a'<=c<='\r' ? print(io,'\\',"abtnvfr"[Int(c)-6]) :
		print(io,c)
	end
end
esc_str2(s::AbstractString)=sprint(esc_str2,s;sizehint=lastindex(s))
function unesc_str2(io::IO,s::AbstractString)
	a=Iterators.Stateful(s)
	for c in a
		if !isempty(a) && c=='\\'
			c=popfirst!(a)
			if c=='x'||c=='u'||c=='U'
				n=k=0
				m=c=='x' ? 2 : c=='u' ? 4 : 8
				while (k+=1)<=m && !isempty(a)
					nc=peek(a)
					n='0'<=nc<='9' ? n<<4+(nc-'0') : 'a'<=nc<='f' ? n<<4+(nc-'a'+10) : 'A'<=nc<='F' ? n<<4+(nc-'A'+10) : break
					popfirst!(a)
				end
				if k==1 || n>0x10ffff
					u=m==4 ? 'u' : 'U'
					throw(ArgumentError("invalid $(m==2 ? "hex (\\x)" : "unicode (\\$u)") escape sequence"))
				end
				if m==2 # \x
					write(io,UInt8(n))
				else
					print(io,Char(n))
				end
			elseif '0'<=c<='7'
				k=1
				n=c-'0'
				while (k+=1)<=3 && !isempty(a)
					c=peek(a)
					n=('0'<=c<='7') ? n<<3+c-'0' : break
					popfirst!(a)
				end
				if n>255 throw(ArgumentError("octal escape sequence out of range")) end
				write(io,UInt8(n))
			else
				print(io,c=='a' ? '\a' : c=='b' ? '\b' : c=='t' ? '\t' : c=='n' ? '\n' : c=='v' ? '\v' : c=='f' ? '\f' : c=='r' ? '\r' : c=='e' ? '\e' :
				c=='/' ? '\\' : c=='=' ? '"' : c=='-' ? '\'' : c=='|' ? '`' : throw(ArgumentError("invalid escape sequence \\$c")))
			end
		else
			print(io,c)
		end
	end
end
unesc_str2(s::AbstractString)=sprint(unesc_str2,s;sizehint=lastindex(s))


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to use MathOptInterface and
# Knitro to solve the following # simple nonlinear optimization problem.
# This model is test problem # HS15 from the Hock & Schittkowski collection.
#
# min   100(x1 - x0^2)^2 +(1 - x0)^2
# s.t.  x0 x1 >= 1
#       x0 + x1^2 >= 0
#       x0 <= 0.5
#
# The standard start point(-2, 1) usually converges to the standard
# minimum at(0.5, 2.0), with final objective = 306.5.
# Sometimes the solver converges to another local minimum
# at(-0.79212, -1.26243), with final objective = 360.4.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO, MathOptInterface
using Test
const MOI = MathOptInterface

struct HS15 <: MOI.AbstractNLPEvaluator
    enable_hessian::Bool
end

MOI.features_available(d::HS15) = [:Grad, :Hess]
MOI.initialize(d::HS15, features) = nothing

MOI.jacobian_structure(d::HS15) = []
MOI.hessian_lagrangian_structure(d::HS15) = Tuple{Int64,Int64}[(1,1), (1, 2), (2, 2)]

#*------------------------------------------------------------------*
#*     FUNCTION callbackEvalF                                       *
#*------------------------------------------------------------------*
function MOI.eval_objective(d::HS15, x)
    return 100*(x[2] - x[1]^2)^2 + (1 - x[1])^2
end

function MOI.eval_constraint(::HS15, g, x)
    nothing
end
function MOI.eval_constraint_jacobian(::HS15, jac_g, x)
    nothing
end

#*------------------------------------------------------------------*
#*     FUNCTION callbackEvalG                                       *
#*------------------------------------------------------------------*
function MOI.eval_objective_gradient(d::HS15, grad_f, x)
    # Evaluate gradient of nonlinear objective
    dTmp = x[2] - x[1] * x[1]
    grad_f[1] = (-400.0 * dTmp * x[1]) - (2.0 * (1.0 - x[1]))
    grad_f[2] = 200.0 * dTmp
end

#*------------------------------------------------------------------*
#*     FUNCTION callbackEvalH                                       *
#*------------------------------------------------------------------*
function MOI.eval_hessian_lagrangian(d::HS15, H, x, σ, μ)
    H[1] = σ * ((-400.0 * x[2]) + (1200.0 * x[1] * x[1]) + 2.0) #(0,0)
    H[2] = σ * (-400.0 * x[1])  #(0,1)
    H[3] = σ * 200.0            #(1,1)

    return 0
end

#*------------------------------------------------------------------*
#*     main                                                         *
#*------------------------------------------------------------------*

# Initialize Knitro.
options = joinpath(dirname(@__FILE__), "..", "examples", "knitro.opt")
# Define a Knitro solver instance through MOI.
solver = KNITRO.Optimizer(outlev=3, opttol=1e-8)
lb =[]; ub=[]
block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), HS15(false), true)

n = 2
v = MOI.add_variables(solver, n)

u = [0.5, Inf]
start = [-2., 1.]

for i in 1:2
    MOI.add_constraint(solver, MOI.SingleVariable(v[i]), MOI.LessThan(u[i]))
    MOI.set(solver, MOI.VariablePrimalStart(), v[i], start[i])
end

# Add the constraints and set their lower bounds.
m = 2
# First constraint: x0 * x1 >= 1.
cf1 = MOI.ScalarQuadraticFunction{Float64}(
        MOI.ScalarAffineTerm.(0.0, v),
        [MOI.ScalarQuadraticTerm(1., v[1], v[2])],
        0.)
c1 = MOI.add_constraint(solver, cf1, MOI.GreaterThan{Float64}(1.))

# Second constraint: x0 + x1^2 >= 0.
cf2 = MOI.ScalarQuadraticFunction(
                                 [MOI.ScalarAffineTerm(1.0, v[1])],
                                 [MOI.ScalarQuadraticTerm(1., v[2], v[2])],
                                 0.)
c2 = MOI.add_constraint(solver, cf2, MOI.GreaterThan(0.))

MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE)
# Define NLP structure.
MOI.set(solver, MOI.NLPBlock(), block_data)

# Solve the problem.
MOI.optimize!(solver)

# Free solver environment properly
MOI.empty!(solver)


text:
module OlfactoryTopographicBehavior

using IntervalSets: minimum
using DataFrames
using IntervalSets
using Unitful: Hz, s, uconvert, Quantity
using DSP
using Statistics
using Requires

export TrialData, collect_by_trialtype, is_lick, session_frame, trialdatas, trialrange

include("types.jl")
include("utils.jl")
include("api.jl")
include("sniffs.jl")

## __init__

function __init__()
    @require PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" include("plotting.jl")
end

end


text:
export TensorOp, DiagonalOp,
    IdentityOp, TensorOpProduct,
    TensorOpProductSum,
    TensorOpProductKernel,
    create, destroy, sigma, jz, jplus, jminus,
    share, dim, tensor,
    iidx, oidx, all_close, apply_serial!

import Base: (*), transpose, conj, ctranspose, kron, A_mul_B!, sparse, full,
              convert, promote_rule

abstract TensorOp

immutable DiagonalOp{T<:Number} <: TensorOp
    # T: data type
    # J: number of zeros before first non-zero
    # K: number of zeros after last non-zero
    # M: Shift offset
    # vals: diagonal elements between the non-zeros
    vals::Vector{T}
    jkm::NTuple{3,Int64}
    # expr::Expr
    # DiagonalOp{T}(vals::Vector{T}, jkm, expr = :()) = new{T}(vals, jkm, expr)
end

# function DiagonalOp(T, op::DiagonalOp)
#   TT = typeof(op)
#   J,K,M = TT.parameters[2:end]
#   DiagonalOp{T,J,K,M}(convert(Vector{T}, op.vals))
# end
#

type IdentityOp{N} <: TensorOp
  # Identity operator with given dimension N
end

promote_rule{T,S}(::Type{DiagonalOp{T}}, ::Type{DiagonalOp{S}}) = DiagonalOp{promote_type(T,S)}
convert{T}(::Type{DiagonalOp{T}}, op::DiagonalOp) = DiagonalOp(convert(Vector{T}, op.vals), op.jkm)
convert{T,N}(::Type{DiagonalOp{T}}, op::IdentityOp{N}) = DiagonalOp(ones(T, N), (0,0,0))

function all_close(op1::DiagonalOp, op2::DiagonalOp; tol=1e-8)
  dim(op1) == dim(op2) || error("Incompatible shapes")
  op1.jkm == op2.jkm && norm(op1.vals-op2.vals) < tol
end

all_close{T}(op1::IdentityOp, op2::DiagonalOp{T}; kw...) = all_close(convert(DiagonalOp{T},op1),op2; kw...)
all_close{N}(op1::IdentityOp{N}, op2::IdentityOp{N}; kw...) = true
all_close{T}(op1::DiagonalOp{T}, op2::IdentityOp; kw...) = all_close(convert(DiagonalOp{T},op2),op1; kw...)

type TensorOpProduct{T,K}
    # Product of single degree tensor ops
    ops::NTuple{K,Union{DiagonalOp{T},IdentityOp}}
    kernel
end
TensorOpProduct(ops::TensorOp...) = TensorOpProduct(promote_diags(ops...), nothing)

function promote_diags(ops::TensorOp...)
  opsv = collect(ops)
  nontrivial=[isa(op, DiagonalOp) for op=opsv]
  promoted_opsv = Array(Any,length(ops))
  promoted_opsv[nontrivial] = collect(promote(opsv[nontrivial]...))
  promoted_opsv[~nontrivial] = opsv[~nontrivial]
  tuple(promoted_opsv...)
end

function convert_diags{T}(S::Type{DiagonalOp{T}}, ops::TensorOp...)
  nops = Any[]
  for op in ops
    if isa(op, DiagonalOp)
      push!(nops, convert(S, op))
    else
      push!(nops, op)
    end
  end
  tuple(nops...)
end

promote_rule{T,S,K}(::Type{TensorOpProduct{T,K}}, ::Type{TensorOpProduct{S,K}}) = TensorOpProduct{promote_type(T,S),K}
promote_rule{T,S}(::Type{TensorOpProduct{T}}, ::Type{TensorOpProduct{S}}) = TensorOpProduct{promote_type(T,S)}
convert{T,S,K}(::Type{TensorOpProduct{T,K}}, top::TensorOpProduct{S,K}) = TensorOpProduct(convert_diags(DiagonalOp{T}, top.ops...)...)
convert{T}(::Type{TensorOpProduct{T}}, top::TensorOpProduct) = TensorOpProduct(convert_diags(DiagonalOp{T}, top.ops...)...)

all_close(top1::TensorOpProduct, top2::TensorOpProduct; tol=1e-8) = all(Bool[all_close(op1, op2; tol=tol) for (op1,op2) in zip(top1.ops, top2.ops)])

immutable TensorOpProductKernel{K}
    # Product of single degree tensor ops
    ops::NTuple{K,Union{DenseVector,IdentityOp}}
    iidx::NTuple{K,Range{Int64}}
    oidx::NTuple{K,Range{Int64}}
end

function TensorOpProductKernel(tp::TensorOpProduct)
    iidxs = iidx(tp)
    oidxs = oidx(tp)
    ops_red = [(isa(op, DiagonalOp)? op.vals : op) for op in tp.ops]
    M = length(tp.ops)
    k = TensorOpProductKernel{M}(tuple(ops_red...), tuple(iidxs...), tuple(oidxs...))
    tp.kernel = k
    k
end

type TensorOpProductSum{T,K,N}
    # Sum of product of single degree tensor ops
    tops::NTuple{N,NTuple{K,Union{IdentityOp,DiagonalOp{T}}}}
    kernel
end

function TensorOpProductSum(tops::TensorOpProduct...)
  tops = promote(tops...)
  topops = [top.ops for top in tops]
  TensorOpProductSum(tuple(topops...), nothing)
end

all_close{T,S,K,N}(tops1::TensorOpProductSum{T,K,N},
          tops2::TensorOpProductSum{S,K,N}; tol=1e-8) = all(
            Bool[all_close(tops1.tops[n][k], tops2.tops[n][k]; tol=tol) for n=1:N, k=1:K])

#
# function share(a::SharedArray)
#     a
# end
#
# function share(a::Array)
#     as = SharedArray(eltype(a),size(a)...)
#     as[:] = a[:]
#     as
# end
#
# function share(op::IdentityOp)
#     op
# end
#
#
# function share{T,J,K,M}(op::DiagonalOp{T,J,K,M})
#     DiagonalOp{T,J,K,M}(share(op.vals))
# end
#
#
# function share(tp::TensorOpProduct)
#     TensorOpProduct(map(share, tp.ops))
# end
#
# function share(tpk::TensorOpProductKernel)
#     TensorOpProductKernel{length(tpk.ops)}(tuple(map(share, tpk.ops)), tpk.iidx, tpk.oidx)
# end
#

function dim(d::DiagonalOp)
  J,K,M = d.jkm
  J+length(d.vals)+K
end

function dim{N}(::IdentityOp{N})
    N
end

# range of input indices that are processed in own dimension.
iidx(d::DiagonalOp) = (1:length(d.vals))+(d.jkm[1]-d.jkm[3])
iidx{N}(::IdentityOp{N}) = 1:N

# range of output indices that are processed in own dimension.
oidx(d::DiagonalOp) =  (1:length(d.vals))+d.jkm[1]
oidx{N}(::IdentityOp{N})=1:N

for fn in [:iidx, :oidx]
    @eval $fn(tp::TensorOpProduct) = map($fn, tp.ops)
end

# conversion to sparse matrices
sparse{N}(::IdentityOp{N})= speye(N)
function sparse{T}(d::DiagonalOp{T})
    N = dim(d)
    J,K,M = d.jkm
    spdiagm([zeros(T,J);d.vals; zeros(T,K)], 0,N,N) * spdiagm(ones(N-abs(M)),-M, N, N)
end
sparse(tp::TensorOpProduct) = reduce(kron, (SparseMatrixCSC)[sparse(op) for op in reverse(tp.ops)])

# conversion of full matrices
full(d::Union{TensorOp, TensorOpProduct}) = d |> sparse |> full

create(N::Int) = DiagonalOp{Float64}(sqrt(1:(N-1)),(1,0,1))
destroy(N::Int) = DiagonalOp{Float64}(sqrt(1:(N-1)),(0,1,-1))
sigma(N::Int,j::Int,k::Int) = DiagonalOp{Float64}([1.],(j-1,N-j,j-k))

function jz(J::Int,trunc=-1)
  if trunc > 0
    N = trunc
  else
    N = 2J
  end
  DiagonalOp{Float64}(1.*collect(-J:(-J+N)),(0,0,0))
end

function jpm(J::Int, trunc=-1, pm=:(+))
  if trunc > 0
    N = trunc
  else
    N = 2J
  end
  rng = -J:(-J+N-1)
  diag = sqrt(J*(J+1) - rng .* (rng+1))
  if pm === :(+)
    DiagonalOp{Float64}(diag, (1,0,1))
  elseif pm === :(-)
    DiagonalOp{Float64}(diag, (0,1,-1))
  else
    error("Can only create J_+ and J_-")
  end
end

jplus(J, trunc=-1) = jpm(J,trunc, :(+))
jminus(J, trunc=-1) = jpm(J,trunc, :(-))

transpose(ii::IdentityOp) = ii
transpose(d::DiagonalOp) =  DiagonalOp(d.vals,(d.jkm[1]-d.jkm[3],d.jkm[2]+d.jkm[3],-d.jkm[3]))
conj(ii::IdentityOp) = ii
conj(d::DiagonalOp) = DiagonalOp(conj(d.vals), d.jkm)
ctranspose(op::TensorOp)= transpose(conj(op))

for fn in [:transpose, :conj, :ctranspose]
    @eval function $fn(tp::TensorOpProduct)
        TensorOpProduct([$fn(op) for op in tp.ops]...)
    end
end

# Multiplication of diagonal tensors on same degree of freedom
function (*){N,M}(a::IdentityOp{N}, b::IdentityOp{M})
    N == M || error("Wrong shape for multiplication: ($N x $N) * ($M x $M)")
    b
end

function (*){N}(a::IdentityOp{N}, b::TensorOp)
    M = dim(b)
    N == M || error("Wrong shape for multiplication: ($N x $N) * ($M x $M)")
    b
end

function (*){M}(a::TensorOp, b::IdentityOp{M})
    N = dim(a)
    N == M || error("Wrong shape for multiplication: ($N x $N) * ($M x $M)")
    a
end

function (*){S,T}(a::DiagonalOp{S}, b::DiagonalOp{T})
    N, M = dim(a), dim(b)
    J1,K1,M1 = a.jkm
    J2,K2,M2 = b.jkm
    N == M || error("Wrong shape for multiplication: ($N x $N) * ($M x $M)")
    da = [zeros(S, J1); a.vals; zeros(S,K1)]
    db = [zeros(T, J2); b.vals; zeros(T,K2)]
    if M1 >= 0
        dbshifted = [zeros(T, M1); db[1:(M-M1)]]
    else
        dbshifted = [db[(abs(M1)+1):M];zeros(T, abs(M1))]
    end
    dab = da .* dbshifted
    q, r = find(abs(dab) .> 0)[[1,end]]
    J = q - 1
    K = N - r
    DiagonalOp{eltype(dab)}(dab[q:r],(J,K,M1+M2))
end

# construction of tensor products
function tensor(a::TensorOp...)
    TensorOpProduct(a...)
end

function tensor(a::TensorOp, b::TensorOpProduct)
    TensorOpProduct(a, b.ops...)
end

function tensor(a::TensorOpProduct, b::TensorOp)
    TensorOpProduct(a.ops..., b)
end

# construction of tensor products for states
function tensor{S,T,M,N}(A::Array{T,M}, B::Array{S,N})
    broadcast(*, reshape(A, size(A)...,ones(Int64,N)...),reshape(B,ones(Int64,M)...,size(B)...))
end

# In-place mat-vec product without coefficients
A_mul_B!(C::DenseArray, A::Union{TensorOpProduct,TensorOpProductKernel}, B::DenseArray) = A_mul_B!(1, A, B, 0, C)

# In-place mat-vec product with coefficients
function A_mul_B!(a, A::TensorOpProduct, B::DenseArray, c, C::DenseArray)
    k = A.kernel === nothing ? TensorOpProductKernel(A) : A.kernel
    A_mul_B!(a, k, B, c, C)
end

function A_mul_B!{K}(a, A::TensorOpProductKernel{K}, B::DenseArray, c, C::DenseArray)
    ndims(B) == K || error("Mismatch between tensor factors")
    size(B) === size(C) || error("Dimension mismatch between input and output")
    if c != 1
      if c == 0
        fill!(C, 0)
      else
        scale!(C, c)
      end
    end
    top_kernel!(sub(C, A.oidx...), sub(B, A.iidx...), a, A.ops...)
    C
end

# code for generating mat-vec loops
@generated function top_kernel!(sC, sB, a, Avs::Union{DenseVector,IdentityOp}...)
   code = top_kernel_code(sC, sB, a, Avs)
  #  println(code)
   code
end
function top_kernel_code(sC, sB, a, Avs)
    K = length(Avs)
    @assert K >= 1
    loop_vars = [symbol("k_$(j)") for j=1:K]
    sCref = Expr(:ref, :sC, loop_vars...)
    sBref = Expr(:ref, :sB, loop_vars...)
    Avrefs = []
    for kk=1:K
        if Avs[kk] <: DenseArray
            push!(Avrefs, Expr(:ref, Expr(:ref, :Avs, kk), loop_vars[kk]))
        else
            @assert Avs[kk] <: IdentityOp
        end
    end

    innerloop_prod = Expr(:call, :*, :a, sBref, Avrefs...)
    innerloop = quote
        @simd for $(loop_vars[1]) = 1:size(sC,1)
            @inbounds $sCref += $innerloop_prod
        end
    end
    loop = innerloop
    for kk=2:K
        loop = quote
            for $(loop_vars[kk]) in 1:size(sC,$kk)
                $loop
            end
        end
    end
    loop
end

function apply_serial!(coeffs, toks, B, C, c=0)
    A_mul_B!(coeffs[1], toks[1], B, c, C)
    for kk=2:length(toks)
        A_mul_B!(coeffs[kk], toks[kk], B, 1, C)
    end
    C
end

function apply_serial!(coeffs, toks::TensorOpProductSum, B::DenseArray, C::DenseArray, c=0)
    apply_topsum!(coeffs, toks.tops, B, c, C)
end

@generated function apply_topsum!{T,N,K}(ts, tops::NTuple{N,NTuple{K,Union{IdentityOp,DiagonalOp{T}}}}, B::DenseArray{T,K}, c, C::DenseArray{T,K})
  code = apply_topsum_code(ts, tops, B, c, C)
    # println(code)
    code
end

function apply_topsum_code(ts, tops, B, c, C)
    # function diag_op_code(has_identity, offsets::Matrix{Int64}, nontrivial_factors::Matrix{Bool})
    # N = number of nontrivial terms in sum
    # D = number of tensor dimensions.
    N = length(tops.parameters)
    N >= 1 || error("Trivial Operator sum passed.")
    D = length(tops.parameters[1].parameters)

    nontrivial_factors = Array(Bool, N, D)
    for n=1:N
        for d=1:D
            nontrivial_factors[n,d] = tops.parameters[n].parameters[d] <: DiagonalOp
        end
    end

    IIDXs = [symbol("IIDX_$(n)") for n=1:N]
    kds = [symbol("k_$d") for d=1:D]
    flags = [symbol("apply_t$(n)_d$(d)") for n=1:N, d=1:D]

    diags = [symbol("diag_t$(n)_d$(d)") for n=1:N, d=1:D]
    jnds = [symbol("j_$(n)_$(d)") for n=1:N, d=1:D]
    knds = [symbol("k_$(n)_$(d)") for n=1:N, d=1:D]
    mnds = [symbol("m_$(n)_$(d)") for n=1:N, d=1:D]
    kdsmmnds = [nontrivial_factors[n,d] ? :($(kds[d])-$(mnds[n,d])) : kds[d] for n=1:N, d=1:D]

    diag_assignments = quote
        # Assign DiagonalOp.vals and .jkm properties to variables
    end
    for n=1:N
        for d=1:D
            if nontrivial_factors[n,d]
                push!(diag_assignments.args, quote
                    $(diags[n,d]) = tops[$n][$d].vals
                    ($(jnds[n,d]),$(knds[n,d]),$(mnds[n,d])) = tops[$n][$d].jkm
#                     println($n, ", ", $d, ": ", tops[$n][$d].jkm)
                    end)
            end
        end
    end

    iidx_init = quote
        IIDX=1
    end

    iidx_inc = quote
        IIDX += 1
    end

    for n=1:N
        terms = Any[]
        for d=1:D
            if nontrivial_factors[n,d]
                push!(terms, :($(mnds[n,d]) * stride(C,$d)))
            end
        end
        if length(terms) > 1
            offset = :(1-$(Expr(:call, :+, terms...)))
        elseif length(terms) == 1
            offset = :(1-$(terms[1]))
        else
            offset = :(1)
        end
        push!(iidx_init.args, :($(IIDXs[n])= $offset))
        push!(iidx_inc.args, :($(IIDXs[n]) += 1))
    end

    flag_updates = Any[quote end for d=1:D]

    for n=1:N
        if nontrivial_factors[n,D]
            push!(flag_updates[D].args,:(
            $(flags[n,D]) = ($(kdsmmnds[n,D]) in 1:size(C,$D)) && $(kds[D]) in (1+$(jnds[n,D])):(size(C,$D)-$(knds[n,D]));
            ))
        else
            push!(flag_updates[D].args,:($(flags[n,D]) = true))
        end

    end
    for d=1:D-1
        for n=1:N
            if nontrivial_factors[n,d]
                push!(flag_updates[d].args,:(
                $(flags[n,d]) = $(flags[n,d+1]) && ($(kdsmmnds[n,d]) in 1:(size(C,$d))) && $(kds[d]) in (1+$(jnds[n,d])):(size(C,$d)-$(knds[n,d]));
                ))

            else
                push!(flag_updates[d].args,:($(flags[n,d]) = $(flags[n,d+1])))
            end
        end
    end

    innerloop_stmts = quote
        C[IIDX] = c * C[IIDX]
    end

    for n=1:N

        tensor_elements = Expr(:call, :(*), :(ts[$n]), Expr(:ref, :B, IIDXs[n]))

        for d=1:D
            if nontrivial_factors[n,d]
                push!(tensor_elements.args,:($(diags[n,d])[$(kds[d])-$(jnds[n,d])]))
            end
        end

        push!(innerloop_stmts.args,
            :(
            if $(flags[n,1])
            C[IIDX] += $tensor_elements
#             else
#             println($n, ":    ", $(kds[1]),", ", $(kds[2]))
        end
            )
        )

    end

    innerloop = quote
        @simd for $(kds[1])=1:size(C,1)
            $(flag_updates[1])
            @inbounds begin
                $innerloop_stmts
            end
            $iidx_inc

        end
    end
    loop = innerloop
    for d=2:D
        loop = quote
            for $(kds[d])=1:size(C,$d)
                $(flag_updates[d])
                $loop
            end
        end
    end

    quote
        $diag_assignments
        $iidx_init
        $loop
        C
    end
end


text:
module Reproject

using FITSIO, WCS, Interpolations, SkyCoords
using SkyCoords: lat,lon

include("parsers.jl")
include("utils.jl")
include("core.jl")

export
    reproject

end # module


text:
"""
    Const <: Kernel

Constant kernel
```math
k(x,x') = σ²
```
with signal standard deviation ``σ``.
"""
mutable struct Const <: Kernel
    "Signal variance"
    σ2::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        Const(lσ::Float64)

    Create `Const` with signal standard deviation `exp(lσ)`.
    """
    Const(lσ::Float64) = new(exp(2 * lσ), [])
end

function set_params!(cons::Const, hyp::VecF64)
    length(hyp) == 1 || throw(ArgumentError("Constant kernel only has one parameters"))
    cons.σ2 = exp(2.0*hyp[1])
end

get_params(cons::Const) = Float64[log(cons.σ2)/2.0]
get_param_names(cons::Const) = [:lσ]
num_params(cons::Const) = 1

Statistics.cov(cons::Const) = cons.σ2
function Statistics.cov(cons::Const, x::VecF64, y::VecF64)
    return cov(cons)
end

@inline dk_dlσ(cons::Const) = 2.0*cov(cons)
@inline function dKij_dθp(cons::Const, X::MatF64, i::Int, j::Int, p::Int, dim::Int)
    if p == 1
        return dk_dlσ(cons)
    else
        return NaN
    end
end
@inline function dKij_dθp(cons::Const, X::MatF64, data::EmptyData, i::Int, j::Int, p::Int, dim::Int)
    if p == 1
        return dKij_dθp(cons, X, i, j, p, dim)
    else
        return NaN
    end
end


text:
module LinearControl

    require("../test/test.jl")
end


text:
function walk_directory(fn, directory)
    for file in readdir(directory)
        path = joinpath(directory, file)
        if isdir(path)
            walk_directory(fn, path)
        else
            fn(path)
        end
    end
end