text:
# GurobiMathProgInterface
# Standardized MILP interface

export GurobiSolver

import MathProgBase
const MPB = MathProgBase

mutable struct GurobiMathProgModel <: MPB.AbstractLinearQuadraticModel
    inner::Model
    last_op_type::Symbol  # To support arbitrary order of addVar/addCon
                          # Two possibilities :Var :Con
    changed_constr_bounds::Bool # have we updated the bounds below?
    obj::Vector{Float64} # objective vector
    lb::Vector{Float64} # persistent bounds on constraints to maintain
    ub::Vector{Float64} # abstraction for lb ≤ Ax ≤ ub
    lazycb
    cutcb
    heuristiccb
    infocb
    options
end

function GurobiMathProgModel(env=nothing;options...)
   finalize_env = (env == nothing)
   if env == nothing
       env = Env()
   end
   m = GurobiMathProgModel(Model(env,""; finalize_env=finalize_env), :Con, false, Float64[], Float64[], Float64[], nothing, nothing, nothing, nothing, options)
   setparams!(m)
   return m
end

function Base.copy(m::GurobiMathProgModel)

    m.lazycb == nothing || @warn("Callbacks can't be copied, lazy callback ignored")
    m.cutcb == nothing || @warn("Callbacks can't be copied, cut callback ignored")
    m.heuristiccb == nothing || @warn("Callbacks can't be copied, heuristic callback ignored")
    m.infocb == nothing || @warn("Callbacks can't be copied, info callback ignored")

    return GurobiMathProgModel(copy(m.inner),
                               m.last_op_type,
                               m.changed_constr_bounds,
                               copy(m.obj),
                               copy(m.lb),
                               copy(m.ub),
                               nothing,
                               nothing,
                               nothing,
                               nothing,
                               deepcopy(m.options))
end

function setparams!(m::GurobiMathProgModel)
    # Helper to set the parameters on the model's copy of env rather than
    # modifying the global env (ref: http://www.gurobi.com/support/faqs#P)

    # Set `InfUnbdInfo` to 1 by default so infeasibility rays available
    setparam!(m.inner, "InfUnbdInfo", 1)
    for (name,value) in m.options
        setparam!(m.inner, string(name), value)
    end
end

mutable struct GurobiSolver <: MPB.AbstractMathProgSolver
    env
    options
end
GurobiSolver(env=nothing; kwargs...) = GurobiSolver(env, kwargs)
MPB.LinearQuadraticModel(s::GurobiSolver) = GurobiMathProgModel(s.env; s.options...)
MPB.ConicModel(s::GurobiSolver) = MPB.LPQPtoConicBridge(MPB.LinearQuadraticModel(s))

MPB.supportedcones(::GurobiSolver) = [:Free,:Zero,:NonNeg,:NonPos,:SOC,:SOCRotated]

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

function MPB.setparameters!(m::GurobiMathProgModel; mpboptions...)
    for (optname, optval) in mpboptions
        if optname == :TimeLimit
            setparam!(m.inner, "TimeLimit", optval)
        elseif optname == :Silent
            if optval == true
                setparam!(m.inner,"OutputFlag",0)
            end
        else
            error("Unrecognized parameter $optname")
        end
    end
end

function MPB.loadproblem!(m::GurobiMathProgModel, filename::AbstractString)
    read_model(m.inner, filename)
    m.obj = MPB.getobj(m)
    if checkvalue(m.obj, GRB_INFINITY)
        _objwarning(m.obj)
    end
    _truncateobj!(m.obj)
    m.lb = MPB.getconstrLB(m)
    m.ub = MPB.getconstrUB(m)
end

function MPB.loadproblem!(m::GurobiMathProgModel, A::LinearAlgebra.Adjoint{T, Array{T, 2}},
                        collb, colub, obj, rowlb, rowub, sense) where T
    MPB.loadproblem!(m, collect(A), collb, colub, obj, rowlb, rowub, sense)
end

function MPB.loadproblem!(m::GurobiMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
  # throw away old model but keep env and finalize_env
  env = m.inner.env
  finalize_env = m.inner.finalize_env
  m.inner.finalize_env = false
  free_model(m.inner)
  m.inner = Model(env, "", finalize_env=finalize_env)
  # Re-set options on new model's copy of env
  setparams!(m)

  add_cvars!(m.inner, float(copy(obj)), float(collb), float(colub))
  update_model!(m.inner)

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

  # check if we have any range constraints
  # to properly support these, we will need to keep track of the
  # slack variables automatically added by gurobi.
  rangeconstrs = sum((rowlb .!= rowub) .& (rowlb .> neginf) .& (rowub .< posinf))
  if rangeconstrs > 0
      @warn("Julia Gurobi interface doesn't properly support range " *
            "(two-sided) constraints. See Gurobi.jl issue #14")
      add_rangeconstrs!(m.inner, float(A), float(rowlb), float(rowub))
      # A work around for the additional slack variables introduced and the
      # reformulation bug warning
      append!(obj, zeros(rangeconstrs))
  else
      b = Array{Float64}(undef, length(rowlb))
      senses = Array{Cchar}(undef, length(rowlb))
      for i in 1:length(rowlb)
          if rowlb[i] == rowub[i]
              senses[i] = '='
              b[i] = rowlb[i]
          elseif rowlb[i] > neginf
              senses[i] = '>'
              b[i] = rowlb[i]
          else
              @assert rowub[i] < posinf
              senses[i] = '<'
              b[i] = rowub[i]
          end
      end
      add_constrs!(m.inner, float(A), senses, b)
  end

  m.lb, m.ub = rowlb, rowub
  m.obj = copy(obj)
  _truncateobj!(m.obj)
  update_model!(m.inner)
  MPB.setsense!(m,sense)
end

function _truncateobj(v::Real)
    # Gurobi truncates objective coefficients to +/-GRB_INFINITY
    # We need to do this to the m.obj vector for consistency
    #   see #94
    if v > GRB_INFINITY
        return GRB_INFINITY
    elseif v < -GRB_INFINITY
        return -GRB_INFINITY
    end
    return v
end
_truncateobj!(obj::Vector) = map!(_truncateobj, obj, obj)

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

MPB.getvarLB(m::GurobiMathProgModel)     = get_dblattrarray( m.inner, "LB", 1, num_vars(m.inner))
function MPB.setvarLB!(m::GurobiMathProgModel, l)
    if checkvalue(l, GRB_BOUNDMAX)
        _boundwarning(l, getvarUB(m))
    end
    set_dblattrarray!(m.inner, "LB", 1, num_vars(m.inner), l)
end

MPB.getvarUB(m::GurobiMathProgModel)     = get_dblattrarray( m.inner, "UB", 1, num_vars(m.inner))
function MPB.setvarUB!(m::GurobiMathProgModel, u)
    if checkvalue(u, GRB_BOUNDMAX)
        _boundwarning(getvarLB(m), u)
    end
    set_dblattrarray!(m.inner, "UB", 1, num_vars(m.inner), u)
end

function MPB.getconstrLB(m::GurobiMathProgModel)
    sense = get_charattrarray(m.inner, "Sense", 1, num_constrs(m.inner))
    ret   = get_dblattrarray(m.inner, "RHS", 1, num_constrs(m.inner))
    for i = 1:num_constrs(m.inner)
        if sense[i] == '>' || sense[i] == '='
            # Do nothing
        else
            # LEQ constraint, so LB is -Inf
            ret[i] = -Inf
        end
     end
     return ret
end

function MPB.getconstrUB(m::GurobiMathProgModel)
    sense = get_charattrarray(m.inner, "Sense", 1, num_constrs(m.inner))
    ret   = get_dblattrarray(m.inner, "RHS", 1, num_constrs(m.inner))
    for i = 1:num_constrs(m.inner)
        if sense[i] == '<' || sense[i] == '='
            # Do nothing
        else
            # GEQ constraint, so UB is +Inf
            ret[i] = +Inf
        end
    end
    return ret
end

# setconstrLB!(m::GurobiMathProgModel, lb) = (m.changed_constr_bounds = true; m.last_op_type = :Con; m.lb = copy(lb))
# setconstrUB!(m::GurobiMathProgModel, ub) = (m.changed_constr_bounds = true; m.last_op_type = :Con; m.ub = copy(ub))
MPB.setconstrLB!(m::GurobiMathProgModel, lb) = (m.changed_constr_bounds = true; m.lb = copy(lb))
MPB.setconstrUB!(m::GurobiMathProgModel, ub) = (m.changed_constr_bounds = true; m.ub = copy(ub))

MPB.getobj(m::GurobiMathProgModel, i::Int=1) = get_dblattrarray( m.inner, "Obj", i, num_vars(m.inner))
function MPB.setobj!(m::GurobiMathProgModel, c, i::Int=1)
    if checkvalue(c, GRB_INFINITY)
        _objwarning(c)
    end
    m.obj = copy(c)
    _truncateobj!(m.obj)
    set_dblattrarray!(m.inner, "Obj", i, num_vars(m.inner), c)
end

set_multiobj_n!(m::Gurobi.Model, n::Int) = Gurobi.set_intattr!(m, "NumObj", n)
get_multiobj_n(m::Gurobi.Model)          = Gurobi.get_intattr(m, "NumObj")

function set_multiobj_c!(m::Model, i::Int, c::AbstractVector{Float64})
    _chklen(c, num_vars(m))
    i0 = get_int_param(m, "ObjNumber")
    set_int_param!(m, "ObjNumber", i-1)
    set_dblattrarray!(m, "ObjN", 1, num_vars(m), c)
    set_int_param!(m, "ObjNumber", i0)
end

function get_multiobj_c(m::Model, i::Int)
    i0 = get_int_param(m, "ObjNumber")
    set_int_param!(m, "ObjNumber", i-1)
    c = get_dblattrarray(m, "ObjN", 1, num_vars(m))
    set_int_param!(m, "ObjNumber", i0)
    return c
end

function set_multiobj_priority!(m::Model, i::Int, priority::Int)
    i0 = get_int_param(m, "ObjNumber")
    set_int_param!(m, "ObjNumber", i-1)
    set_intattr!(m, "ObjNPriority", priority)
    set_int_param!(m, "ObjNumber", i0)
end

function get_multiobj_priority(m::Model, i::Int)::Int
    i0 = get_int_param(m, "ObjNumber")
    set_int_param!(m, "ObjNumber", i-1)
    p = get_intattr(m, "ObjNPriority")
    set_int_param!(m, "ObjNumber", i0)
    return p
end

function set_multiobj_weight!(m::Model, i::Int, w::Float64)
    i0 = get_int_param(m, "ObjNumber")
    set_int_param!(m, "ObjNumber", i-1)
    set_dblattr!(m, "ObjNWeight", w)
    set_int_param!(m, "ObjNumber", i0)
end

function get_multiobj_weight(m::Model, i::Int)::Float64
    i0 = get_int_param(m, "ObjNumber")
    set_int_param!(m, "ObjNumber", i-1)
    w::Float64 = get_dblattr(m, "ObjNWeight")
    set_int_param!(m, "ObjNumber", i0)
    return w
end

function set_multiobj!(m::Model, i::Int, c::AbstractVector{Float64}, p::Int, w::Float64)
    nobj = get_multiobj_n(m)
    if !(1 ≤ i ≤ nobj)
        error("Tried to set objective $i of $nobj. You must use the set_multiobj_n to specify the correct number of objective functions for your problem.")
    end
    set_multiobj_c!(m, i, c)
    set_multiobj_priority!(m, i, p)
    set_multiobj_weight!(m, i, w)
end

function MPB.addvar!(m::GurobiMathProgModel, constridx, constrcoef, l, u, objcoef)
    if m.last_op_type == :Con
        updatemodel!(m)
        m.last_op_type = :Var
    end
    push!(m.obj, _truncateobj(objcoef))
    add_var!(m.inner, length(constridx), constridx, float(constrcoef), float(objcoef), float(l), float(u), GRB_CONTINUOUS)
end
MPB.addvar!(m::GurobiMathProgModel, l, u, objcoef) = MPB.addvar!(m, Int[], Float64[], l, u, objcoef)

function MPB.delvars!(m::GurobiMathProgModel, idx)
    deleteat!(m.obj, idx)
    del_vars!(m.inner, idx)
end

function MPB.addconstr!(m::GurobiMathProgModel, varidx, coef, lb, ub)
    if m.last_op_type == :Var
        updatemodel!(m)
        m.last_op_type = :Con
    end
    if lb == -Inf
        # <= constraint
        add_constr!(m.inner, varidx, coef, '<', ub)
    elseif ub == +Inf
        # >= constraint
        add_constr!(m.inner, varidx, coef, '>', lb)
    elseif lb == ub
        # == constraint
        add_constr!(m.inner, varidx, coef, '=', lb)
    else
        # Range constraint
        error("Adding range constraints not supported yet.")
    end
end
MPB.delconstrs!(m::GurobiMathProgModel, idx) = del_constrs!(m.inner, idx)

MPB.changecoeffs!(m::GurobiMathProgModel, cidx, vidx, val) = chg_coeffs!(m.inner, cidx, vidx, val)

function updatemodel!(m::GurobiMathProgModel)
    update_model!(m.inner)
    if Gurobi.version < v"7.0" && m.obj != MPB.getobj(m)
        @warn("""
            You have encountered a known bug in Gurobi. Any information you query from the model may be incorrect.
            This bug has existed since the first version of Gurobi but is fixed in Gurobi v7.0.

            For more information go to https://github.com/JuliaOpt/Gurobi.jl/issues/60
            Please leave a comment stating that you encountered this bug! We would like to know how prevalent it is.
        """)
    end
    if m.changed_constr_bounds
    # update lower/upper bounds on linear constraints (if they're consistent...)
        sense = get_charattrarray(m.inner, "Sense", 1, num_constrs(m.inner))
        rhs   = get_dblattrarray( m.inner, "RHS",   1, num_constrs(m.inner))
        lb, ub = m.lb, m.ub
        @assert (n_rows = num_constrs(m.inner)) == length(lb) == length(ub) ==
                                                   length(sense) == length(rhs)
        for i = 1:n_rows
            if -Inf < lb[i] == ub[i] < Inf
                sense[i] = '='
                rhs[i]   = lb[i]
            elseif (-Inf < lb[i] && ub[i] < Inf)
                error("Gurobi.jl does not currently support range constraints")
            elseif (lb[i] == -Inf && ub[i] < Inf)
                sense[i] = '<'
                rhs[i]   = ub[i]
            elseif (lb[i] > -Inf && ub[i] == Inf)
                sense[i] = '>'
                rhs[i]   = lb[i]
            else
                error("Internal error.")
            end
        end
        set_charattrarray!(m.inner, "Sense", 1, num_constrs(m.inner), sense)
        set_dblattrarray!(m.inner, "RHS", 1, num_constrs(m.inner), rhs)
        m.changed_constr_bounds = false
        update_model!(m.inner)
    end
end

MPB.getconstrmatrix(m::GurobiMathProgModel) = get_constrmatrix(m.inner)

function MPB.setsense!(m::GurobiMathProgModel, sense)
  if sense == :Min
    set_sense!(m.inner, :minimize)
  elseif sense == :Max
    set_sense!(m.inner, :maximize)
  else
    error("Unrecognized objective sense $sense")
  end
end
function MPB.getsense(m::GurobiMathProgModel)
  v = get_intattr(m.inner, "ModelSense")
  if v == -1
    return :Max
  else
    return :Min
  end
end

MPB.numvar(m::GurobiMathProgModel)    = (updatemodel!(m); num_vars(m.inner))
MPB.numconstr(m::GurobiMathProgModel) = num_constrs(m.inner) + num_qconstrs(m.inner)
MPB.numlinconstr(m::GurobiMathProgModel) = num_constrs(m.inner)
MPB.numquadconstr(m::GurobiMathProgModel) = num_qconstrs(m.inner)

function MPB.optimize!(m::GurobiMathProgModel)
    # set callbacks if present
    if m.lazycb != nothing || m.cutcb != nothing || m.heuristiccb != nothing || m.infocb != nothing
        updatemodel!(m)
        if !is_mip(m.inner)
            @warn("Gurobi ignores branch-and-bound callbacks when no discrete elements are present in the model.")
        end
        setmathprogcallback!(m)
    end
    if m.lazycb != nothing
      setparam!(m.inner, "LazyConstraints", 1)
    end
    updatemodel!(m)
    optimize(m.inner)
end

function MPB.status(m::GurobiMathProgModel)
  s = get_status(m.inner)
  if s == :optimal
    return :Optimal
  elseif s == :infeasible
    return :Infeasible
  elseif s == :unbounded
    return :Unbounded
  elseif s == :inf_or_unbd
    @warn("Gurobi reported infeasible or unbounded. Set InfUnbdInfo=1 for more specific status.")
    return :InfeasibleOrUnbounded
  elseif s == :iteration_limit || s == :node_limit || s == :time_limit || s == :solution_limit
    return :UserLimit
  elseif s == :numeric
    return :Numeric
  elseif s == :suboptimal
    return :Suboptimal # not very useful status
  elseif s == :interrupted # ended by user?
    return :UserLimit
  elseif s == :inprogress
    return :InProgress
  elseif s == :user_obj_limit
    return :UserObjLimit
  elseif s == :cutoff
    return :Cutoff
  else
    error("Unrecognized solution status: $s")
  end
end

MPB.getobjval(m::GurobiMathProgModel)   = get_objval(m.inner)
MPB.getobjbound(m::GurobiMathProgModel) = get_objbound(m.inner)
MPB.getsolution(m::GurobiMathProgModel) = get_solution(m.inner)

function MPB.getconstrsolution(m::GurobiMathProgModel)
    sense = get_charattrarray(m.inner, "Sense", 1, num_constrs(m.inner))
    rhs   = get_dblattrarray( m.inner, "RHS",   1, num_constrs(m.inner))
    slack = get_dblattrarray( m.inner, "Slack", 1, num_constrs(m.inner))
    ret   = zeros(num_constrs(m.inner))
    for i = 1:num_constrs(m.inner)
        if sense[i] == '='
            # Must be equal to RHS if feasible
            ret[i] = rhs[i]
        else
            # Slack variable is non-negative for <= constraints, and non-positive for >= constraints
            ret[i] = rhs[i] - slack[i]
        end
    end
    return ret
end

function MPB.getreducedcosts(m::GurobiMathProgModel)
    if is_qcp(m.inner) && get_int_param(m.inner, "QCPDual") == 0
        return fill(NaN, num_vars(m.inner))
    else
        return get_dblattrarray(m.inner, "RC", 1, num_vars(m.inner))
    end
end

function MPB.getconstrduals(m::GurobiMathProgModel)
    if is_qcp(m.inner) && get_int_param(m.inner, "QCPDual") == 0
        return fill(NaN, num_constrs(m.inner))
    else
        return get_dblattrarray(m.inner, "Pi", 1, num_constrs(m.inner))
    end
end

function MPB.getquadconstrduals(m::GurobiMathProgModel)
    if is_qcp(m.inner) && get_int_param(m.inner, "QCPDual") == 0
        return fill(NaN, num_qconstrs(m.inner))
    else
        return get_dblattrarray(m.inner, "QCPi", 1, num_qconstrs(m.inner))
    end
end

MPB.getinfeasibilityray(m::GurobiMathProgModel) = -get_dblattrarray(m.inner, "FarkasDual", 1, num_constrs(m.inner)) # note sign is flipped
MPB.getunboundedray(m::GurobiMathProgModel) = get_dblattrarray(m.inner, "UnbdRay", 1, num_vars(m.inner))

MPB.getbasis(m::GurobiMathProgModel) = get_basis(m.inner)

MPB.getrawsolver(m::GurobiMathProgModel) = 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 MPB.setvartype!(m::GurobiMathProgModel, vartype::Vector{Symbol})
    # do this to make sure we deal with new columns
    updatemodel!(m)
    nvartype = map(x->rev_var_type_map[x], vartype)
    set_charattrarray!(m.inner, "VType", 1, length(nvartype), nvartype)
    updatemodel!(m) # otherwise getvartype! will return old values
end

function MPB.getvartype(m::GurobiMathProgModel)
    ret = get_charattrarray(m.inner, "VType", 1, num_vars(m.inner))
    map(x->var_type_map[x], ret)
end

function MPB.setwarmstart!(m::GurobiMathProgModel, v)
    for j = 1:length(v)
        if isnan(v[j])
            v[j] = 1e101  # GRB_UNDEFINED
        end
    end
    set_dblattrarray!(m.inner, "Start", 1, num_vars(m.inner), v)
end

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

# Callbacks

MPB.setlazycallback!(m::GurobiMathProgModel,f) = (m.lazycb = f)
MPB.setcutcallback!(m::GurobiMathProgModel,f) = (m.cutcb = f)
MPB.setheuristiccallback!(m::GurobiMathProgModel,f) = (m.heuristiccb = f)
MPB.setinfocallback!(m::GurobiMathProgModel,f) = (m.infocb = f)

mutable struct GurobiCallbackData <: MPB.MathProgCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    sol::Vector{Float64}  # Used for heuristic callbacks
#    model::GurobiMathProgModel # not needed?
end

function MPB.cbgetmipsolution(d::GurobiCallbackData)
    @assert d.state == :MIPSol
    return cbget_mipsol_sol(d.cbdata, d.where)
end

function MPB.cbgetmipsolution(d::GurobiCallbackData,output)
    @assert d.state == :MIPSol
    return cbget_mipsol_sol(d.cbdata, d.where, output)
end

function MPB.cbgetlpsolution(d::GurobiCallbackData)
    @assert d.state == :MIPNode
    return cbget_mipnode_rel(d.cbdata, d.where)
end

function MPB.cbgetlpsolution(d::GurobiCallbackData, output)
    @assert d.state == :MIPNode
    return cbget_mipnode_rel(d.cbdata, d.where, output)
end

# TODO: macro for these getters?
function MPB.cbgetobj(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_objbst(d.cbdata, d.where)
    elseif d.state == :Intermediate
        return cbget_mip_objbst(d.cbdata, d.where)
    elseif d.state == :MIPSol
        error("Gurobi does not implement cbgetobj when state == MIPSol")
        # https://groups.google.com/forum/#!topic/gurobi/Az_X6Ag-y6k
        # https://github.com/JuliaOpt/MathProgBase.jl/issues/23
        # A complex workaround is possible but hasn't been implemented.
        # return cbget_mipsol_objbst(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

function MPB.cbgetbestbound(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_objbnd(d.cbdata, d.where)
    elseif d.state == :MIPSol
        return cbget_mipsol_objbnd(d.cbdata, d.where)
    elseif d.state == :Intermediate
        return cbget_mip_objbnd(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

function MPB.cbgetexplorednodes(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_nodcnt(d.cbdata, d.where)
    elseif d.state == :MIPSol
        return cbget_mipsol_nodcnt(d.cbdata, d.where)
    elseif d.state == :Intermediate
        return cbget_mip_nodcnt(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

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

function MPB.cbaddcut!(d::GurobiCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode
    cbcut(d.cbdata, convert(Vector{Cint}, varidx), float(varcoef), sense, float(rhs))
end

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

function MPB.cbaddsolution!(d::GurobiCallbackData)
    # Gurobi doesn't support adding solutions on MIPSol.
    # TODO: support this anyway
    @assert d.state == :MIPNode
    cbsolution(d.cbdata, d.sol)
    # "Wipe" solution back to GRB_UNDEFINIED
    for i in 1:length(d.sol)
        d.sol[i] = 1e101  # GRB_UNDEFINED
    end
end

function MPB.cbsetsolutionvalue!(d::GurobiCallbackData,varidx,value)
    d.sol[varidx] = value
end

# breaking abstraction, define our low-level callback to eliminatate
# a level of indirection

function mastercallback(ptr_model::Ptr{Cvoid}, cbdata::Ptr{Cvoid}, where::Cint, userdata::Ptr{Cvoid})

    model = unsafe_pointer_to_objref(userdata)::GurobiMathProgModel
    grbrawcb = CallbackData(cbdata,model.inner)
    if where == CB_MIPSOL
        state = :MIPSol
    elseif where == CB_MIPNODE
        state = :MIPNode
        # skip callback if node is reported to be cut off or infeasible --
        # nothing to do.
        # TODO: users may want this information
        status = cbget_mipnode_status(grbrawcb, where)
        if status != 2
            return convert(Cint,0)
        end
    elseif where == CB_MIP
        state = :Intermediate
    else
        # State with no relevant callbacks
        return convert(Cint,0)
    end

    grbcb = GurobiCallbackData(grbrawcb, state, where, [0.0])
    if model.infocb != nothing
        ret = model.infocb(grbcb)
        if ret == :Exit
            terminate(model.inner)
        end
    end
    if model.cutcb != nothing && state == :MIPNode
        ret = model.cutcb(grbcb)
        if ret == :Exit
            terminate(model.inner)
        end
    end
    if model.heuristiccb != nothing && state == :MIPNode
        grbcb.sol = fill(1e101, MPB.numvar(model))  # GRB_UNDEFINED
        ret = model.heuristiccb(grbcb)
        if ret == :Exit
            terminate(model.inner)
        end
    end
    if model.lazycb != nothing && (state == :MIPSol || state == :MIPNode)
        ret = model.lazycb(grbcb)
        if ret == :Exit
            terminate(model.inner)
        end
    end

    return convert(Cint,0)
end

# User callback function should be of the form:
# callback(cbdata::MPB.MathProgCallbackData)
# return :Exit to indicate an error

function setmathprogcallback!(model::GurobiMathProgModel)
    if Sys.iswindows() && Sys.WORD_SIZE != 64
        error("Callbacks not currently supported on Win32. Use 64-bit Julia with 64-bit Gurobi.")
    end
    grbcallback = @cfunction(mastercallback, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ptr{Cvoid}))
    ret = @grb_ccall(setcallbackfunc, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Any), model.inner.ptr_model, grbcallback, model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

# QCQP

function MPB.setquadobj!(m::GurobiMathProgModel, rowidx, colidx, quadval)
    delq!(m.inner)
    update_model!(m.inner)
    scaledvals = similar(quadval)
    for i in 1:length(rowidx)
      if rowidx[i] == colidx[i]
        # rescale from matrix format to "terms" format
        scaledvals[i] = quadval[i] / 2
      else
        scaledvals[i] = quadval[i]
      end
    end
    add_qpterms!(m.inner, rowidx, colidx, scaledvals)
    update_model!(m.inner)
end

function MPB.addquadconstr!(m::GurobiMathProgModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    if m.last_op_type == :Var
        updatemodel!(m)
        m.last_op_type = :Con
    end
    add_qconstr!(m.inner, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
end

MPB.getsolvetime(m::GurobiMathProgModel) = get_runtime(m.inner)
MPB.getnodecount(m::GurobiMathProgModel) = get_node_count(m.inner)
MPB.getsimplexiter(m::GurobiMathProgModel) = get_dblattr(m.inner, "IterCount")


text:
struct Monomials{dim,order, T}
    indices::Union{Vector{NTuple{dim,Int}},Vector{Int}}
    centroid::Tensors.Vec{dim,T}
    diameter::T
end
function Monomials(dim,degree, centroid, diameter)
    if dim == 1
        return Monomials1D(degree, centroid, diameter)
    elseif dim == 2
        return Monomials2D(degree, centroid, diameter)
    else
        throw("Monomials of dim $dim are not implemented!")
    end
end
function Monomials1D(degree::Int, centroid::Tensors.Vec{1,T}, diameter::T) where {T}
    indices = 0:degree
    Monomials{1,degree, T}(collect(indices), centroid, diameter)
end
function Monomials2D(degree::Int, centroid::Tensors.Vec{2,T}, diameter::T) where {T}
    indices = get_monomial2DIndices(degree)
    Monomials{2,degree, T}(indices, centroid, diameter)
end

function get_monomial2DIndices(degree)
    if degree == 0
        return [(0,0)]
    else
        tmp = Vector{NTuple{2,Int}}()
        a = (degree,0)
        for i in 1:(binomial(degree+1,1))
            push!(tmp,a)
            a = (a[1]-1,a[2]+1)
        end
        return [get_monomial2DIndices(degree-1)...,tmp...]
    end
end

getnbasefunctions(::Monomials{2,order}) where {order} = Int((order+1)*(order+2)/2)
getnbasefunctions(::Monomials{1,order}) where {order} = order+1

# Evaluate monomials
function value(ip::Monomials{dim,order}, k::Int, ξ::Tensors.Vec{dim,T}) where {order, T, dim}
    if k > getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    prod((1/ip.diameter*(ξ-ip.centroid)).^ip.indices[k])
end

function gradient_value(ip::Monomials{dim,order}, k::Int, ξ::Tensors.Vec{dim,T}) where {dim,order,T}
    if k >getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    Tensors.gradient(ξ -> value(ip, k, ξ), ξ)
end

function laplace_value(ip::Monomials{dim,order}, k::Int, ξ::Tensors.Vec{dim,T}) where {dim,order,T}
    if k >getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    Tensors.laplace(ξ -> value(ip, k, ξ), ξ)
end

struct RotMonomials{dim,order, T}
    indices::Union{Vector{NTuple{dim,Int}},Vector{Int}}
    centroid::Tensors.Vec{dim,T}
    diameter::T
end
function RotMonomials(dim,order, centroid, diameter)
    if dim == 2
        return RotMonomials2D(order, centroid, diameter)
    else
        throw("Monomials of dim $dim are not implemented!")
    end
end
function RotMonomials2D(order::Int, centroid::Tensors.Vec{2,T}, diameter::T) where {T}
    indices = get_monomial2DIndices(order-1)
    RotMonomials{2,order, T}(indices, centroid, diameter)
end

getnbasefunctions(::RotMonomials{2,order}) where {order} = Int((order)*(order+1)/2)

# Evaluate monomials

rot(ξ::Tensors.Vec{2}) = Tensors.Vec((-ξ[2],ξ[1]))

function value(ip::RotMonomials{dim,order}, k::Int, ξ::Tensors.Vec{dim,T}) where {order, T, dim}
    if k > getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    rot(ξ)*prod((1/ip.diameter*(ξ-ip.centroid)).^ip.indices[k])
end

function gradient_value(ip::RotMonomials{dim,order}, k::Int, ξ::Tensors.Vec{dim,T}) where {dim,order,T}
    if k >getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    Tensors.gradient(ξ -> value(ip, k, ξ), ξ)
end

function laplace_value(ip::RotMonomials{dim,order}, k::Int, ξ::Tensors.Vec{dim,T}) where {dim,order,T}
    if k >getnbasefunctions(ip);throw(ArgumentError("no shape function $k for interpolation $ip"));end
    Tensors.laplace(ξ -> value(ip, k, ξ), ξ)
end

ext: jl

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
# Problem 1 - Project Euler
# http://projecteuler.net/index.php?section=problems&id=001
function summul(n, x)
    count = div(n, x)
    div(x * count * (count + 1), 2)
end

N = 1000 - 1
println(summul(N, 3) + summul(N, 5) - summul(N, 3 * 5))

max_stars_repo_name: emergent/ProjectEuler

max_issues_repo_name: emergent/ProjectEuler

max_forks_repo_name: emergent/ProjectEuler


text:
f(a, b) = a + b

@code_lowered f(1, 1)
@code_typed f(1, 1)
@code_llvm f(1, 1)
@code_native f(1, 1)

@code_native f(1.0, 3.0)
@code_native f(1.2, 4)

using BenchmarkTools, LinearAlgebra

arr = rand(Float32, 1_000, 10)
arr = rand(Float32, 100_000, 10)

function f(a)
    b = a * a'
    b[diagind(b)] .= 0
    sum(b, dims=2)
end

f(arr)
@benchmark f(arr)

function g(a::Array{Float32, 2})::Array{Float32}
    """
    # vectorized version of the following operations with M (n, m) => NM (n, n) => n
    out = a * a' => creates a a we cannot store in RAM
    out[diagind(out)] .= 0
    out = sum(out, dims=2)
    """
    N, M = size(a)

    out = Array{Float32, 1}(undef, N)

    Threads.@threads for i in 1:N
        total = zero(Float32)
        @inbounds for ii in 1:N
            if (i == ii) continue end
            @inbounds @simd for j in 1:M
                total += a[i, j] * a[ii, j]
            end
        end
        @inbounds out[i] = total
    end
    return out
end

g(arr)
@benchmark g(arr)

using CUDA
function g_kernel(matrix::CuDeviceMatrix{Float32}, out::CuDeviceVector{Float32}, len::Int, N::Int, M::Int)::Nothing
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x

    if (index) < len
        for i in 1:N
            if index == i continue end
            for j in 1:M
                @inbounds out[index] += matrix[index, j] * matrix[i, j]
            end
        end
    end
    return nothing
end

function g(matrix::CuArray{Float32, 2, CUDA.Mem.DeviceBuffer})::Array{Float32}
    len = length(matrix)
    N, M = size(matrix)
    threads_per_block = 256
    blocks = Int(ceil(N / threads_per_block))

    out = CUDA.zeros(Float32, N)

    @cuda threads=threads_per_block blocks=blocks g_kernel(matrix, out, len, N, M)
    return Array(out)
end

arr = CUDA.rand(Float32, 100_000, 10)
g(arr)

@CUDA.time g(arr)

max_stars_repo_name: monkeyusage/brownbaglunch

max_issues_repo_name: monkeyusage/brownbaglunch

max_forks_repo_name: monkeyusage/brownbaglunch


text:
module IUP_CD

export
	# From libcd.jl
	cdCanvas,
	cdCreateCanvas,
	cdKillCanvas,
	cdCanvasGetContext,
	cdCanvasActivate,
	cdCanvasDeactivate,
	cdUseContextPlus,
	cdInitContextPlus,
	cdFinishContextPlus,
	cdContextRegisterCallback,
	cdContextCaps,
	cdContextIsPlus,
	cdContextType,
	cdContext,
	cdContextNativeWindow,
	cdCanvasSimulate,
	cdCanvasFlush,
	cdCanvasClear,
	cdCanvasSaveState,
	cdCanvasRestoreState,
	cdReleaseState,
	cdCanvasSetAttribute,
	cdCanvasGetAttribute,
	cdCanvasPlay,
	cdCanvasGetSize,
	cdCanvasUpdateYAxis,
	cdfCanvasUpdateYAxis,
	cdCanvasInvertYAxis,
	cdfCanvasInvertYAxis,
	cdCanvasMM2Pixel,
	cdCanvasPixel2MM,
	cdfCanvasMM2Pixel,
	cdfCanvasPixel2MM,
	cdCanvasOrigin,
	cdfCanvasOrigin,
	cdCanvasGetOrigin,
	cdfCanvasGetOrigin,
	cdCanvasTransform,
	cdCanvasGetTransform,
	cdCanvasTransformMultiply,
	cdCanvasTransformRotate,
	cdCanvasTransformScale,
	cdCanvasTransformTranslate,
	cdCanvasTransformPoint,
	cdfCanvasTransformPoint,
	cdCanvasClip,
	cdCanvasClipArea,
	cdCanvasGetClipArea,
	cdfCanvasClipArea,
	cdfCanvasGetClipArea,
	cdCanvasIsPointInRegion,
	cdCanvasOffsetRegion,
	cdCanvasGetRegionBox,
	cdCanvasRegionCombineMode,
	cdCanvasPixel,
	cdCanvasMark,
	cdCanvasBegin,
	cdCanvasPathSet,
	cdCanvasEnd,
	cdCanvasLine,
	cdCanvasVertex,
	cdCanvasRect,
	cdCanvasBox,
	cdCanvasArc,
	cdCanvasSector,
	cdCanvasChord,
	cdCanvasText,
	cdfCanvasLine,
	cdfCanvasVertex,
	cdfCanvasRect,
	cdfCanvasBox,
	cdfCanvasArc,
	cdfCanvasSector,
	cdfCanvasChord,
	cdfCanvasText,
	cdCanvasSetBackground,
	cdCanvasSetForeground,
	cdCanvasBackground,
	cdCanvasForeground,
	cdCanvasBackOpacity,
	cdCanvasWriteMode,
	cdCanvasLineStyle,
	cdCanvasLineStyleDashes,
	cdCanvasLineWidth,
	cdCanvasLineJoin,
	cdCanvasLineCap,
	cdCanvasInteriorStyle,
	cdCanvasHatch,
	cdCanvasStipple,
	cdCanvasGetStipple,
	cdCanvasPattern,
	cdCanvasGetPattern,
	cdCanvasFillMode,
	cdCanvasFont,
	cdCanvasGetFont,
	cdCanvasNativeFont,
	cdCanvasTextAlignment,
	cdCanvasTextOrientation,
	cdCanvasMarkType,
	cdCanvasMarkSize,
	cdCanvasVectorText,
	cdCanvasMultiLineVectorText,
	cdCanvasVectorFont,
	cdCanvasVectorTextDirection,
	cdCanvasVectorTextTransform,
	cdCanvasVectorTextSize,
	cdCanvasVectorCharSize,
	cdCanvasVectorFontSize,
	cdCanvasGetVectorFontSize,
	cdCanvasGetVectorTextSize,
	cdCanvasGetVectorTextBounds,
	cdCanvasGetVectorTextBox,
	cdfCanvasVectorTextDirection,
	cdfCanvasVectorTextSize,
	cdfCanvasGetVectorTextSize,
	cdfCanvasVectorCharSize,
	cdfCanvasVectorText,
	cdfCanvasMultiLineVectorText,
	cdfCanvasGetVectorTextBounds,
	cdfCanvasGetVectorTextBox,
	cdCanvasGetFontDim,
	cdCanvasGetTextSize,
	cdCanvasGetTextBox,
	cdCanvasGetTextBounds,
	cdCanvasGetColorPlanes,
	cdCanvasPalette,
	cdCanvasGetImageRGB,
	cdCanvasPutImageRectRGB,
	cdCanvasPutImageRectRGBA,
	cdCanvasPutImageRectMap,
	cdCanvasCreateImage,
	cdKillImage,
	cdCanvasGetImage,
	cdCanvasPutImageRect,
	cdCanvasScrollArea,
	cdCreateBitmap,
	cdKillBitmap,
	cdBitmapGetData,
	cdBitmapSetRect,
	cdCanvasPutBitmap,
	cdCanvasGetBitmap,
	cdBitmapRGB2Map,
	cdEncodeColor,
	cdDecodeColor,
	cdDecodeAlpha,
	cdEncodeAlpha,
	cdRGB2Map,
	cdActivate,
	cdActiveCanvas,
	cdSimulate,
	cdFlush,
	cdClear,
	cdSaveState,
	cdRestoreState,
	cdSetAttribute,
	cdGetAttribute,
	cdGetContext,
	cdRegisterCallback,
	cdPlay,
	cdGetCanvasSize,
	cdUpdateYAxis,
	cdMM2Pixel,
	cdPixel2MM,
	cdOrigin,
	cdClip,
	cdGetClipPoly,
	cdClipArea,
	cdGetClipArea,
	cdPointInRegion,
	cdOffsetRegion,
	cdRegionBox,
	cdRegionCombineMode,
	cdPixel,
	cdMark,
	cdLine,
	cdBegin,
	cdVertex,
	cdEnd,
	cdRect,
	cdBox,
	cdArc,
	cdSector,
	cdChord,
	cdText,
	cdBackground,
	cdForeground,
	cdBackOpacity,
	cdWriteMode,
	cdLineStyle,
	cdLineStyleDashes,
	cdLineWidth,
	cdLineJoin,
	cdLineCap,
	cdInteriorStyle,
	cdHatch,
	cdStipple,
	cdGetStipple,
	cdPattern,
	cdGetPattern,
	cdFillMode,
	cdFont,
	cdGetFont,
	cdNativeFont,
	cdTextAlignment,
	cdTextOrientation,
	cdMarkType,
	cdMarkSize,
	cdMultiLineVectorText,
	cdGetVectorTextSize,
	cdGetVectorTextBounds,
	cdFontDim,
	cdImage,
	cdTextSize,
	cdTextBox,
	cdTextBounds,
	cdGetColorPlanes,
	cdPalette,
	cdGetImageRGB,
	cdPutImageRectRGB,
	cdPutImageRectRGBA,
	cdPutImageRectMap,
	cdCreateImage,
	cdGetImage,
	cdPutImageRect,
	cdScrollArea,
	cdPutBitmap,
	cdGetBitmap,
	cdVersion,
	cdVersionDate,
	cdVersionNumber,
	cdVectorText,
#	cdVectorFont,
	cdVectorTextDirection,
	cdVectorTextTransform,
	cdVectorTextSize,
	cdVectorCharSize,
	#
	# From wd.jl
	#
	wdCanvasWindow,
	wdCanvasGetWindow,
	wdCanvasViewport,
	wdCanvasGetViewport,
	wdCanvasWorld2Canvas,
	wdCanvasWorld2CanvasSize,
	wdCanvasCanvas2World,
	wdCanvasSetTransform,
	wdCanvasGetTransform,
	wdCanvasTranslate,
	wdCanvasScale,
	wdCanvasClipArea,
	wdCanvasGetClipArea,
	wdCanvasIsPointInRegion,
	wdCanvasOffsetRegion,
	wdCanvasGetRegionBox,
	wdCanvasHardcopy,
	wdCanvasPixel,
	wdCanvasMark,
	wdCanvasLine,
	wdCanvasVertex,
	wdCanvasRect,
	wdCanvasBox,
	wdCanvasArc,
	wdCanvasSector,
	wdCanvasChord,
	wdCanvasText,
	wdCanvasPutImageRect,
	wdCanvasPutImageRectRGB,
	wdCanvasPutImageRectRGBA,
	wdCanvasPutImageRectMap,
	wdCanvasPutBitmap,
	wdCanvasLineWidth,
	wdCanvasFont,
	wdCanvasGetFont,
	wdCanvasMarkSize,
	wdCanvasGetFontDim,
	wdCanvasGetTextSize,
	wdCanvasGetTextBox,
	wdCanvasGetTextBounds,
	wdCanvasStipple,
	wdCanvasPattern,
	wdCanvasVectorTextDirection,
	wdCanvasVectorTextSize,
	wdCanvasGetVectorTextSize,
	wdCanvasVectorCharSize,
	wdCanvasVectorText,
	wdCanvasMultiLineVectorText,
	wdCanvasGetVectorTextBounds,
	wdCanvasGetVectorTextBox,
	wdWindow,
	wdGetWindow,
	wdViewport,
	wdGetViewport,
	wdWorld2Canvas,
	wdWorld2CanvasSize,
	wdCanvas2World,
	wdClipArea,
	wdGetClipArea,
	wdGetClipPoly,
	wdPointInRegion,
	wdOffsetRegion,
	wdRegionBox,
	wdHardcopy,
	wdPixel,
	wdMark,
	wdLine,
	wdVertex,
	wdRect,
	wdBox,
	wdArc,
	wdSector,
	wdChord,
	wdText,
	wdPutImageRect,
	wdPutImageRectRGB,
	wdPutImageRectRGBA,
	wdPutImageRectMap,
	wdPutBitmap,
	wdLineWidth,
	wdFont,
	wdGetFont,
	wdMarkSize,
	wdFontDim,
	wdTextSize,
	wdTextBox,
	wdTextBounds,
	wdStipple,
	wdPattern,
	wdVectorTextDirection,
	wdVectorTextSize,
	wdGetVectorTextSize,
	wdVectorCharSize,
	wdVectorText,
	wdMultiLineVectorText,
	wdGetVectorTextBound,
	#
	# From cdiup.jl
	#
	cdContextIup,
	# From cdxdbuf
	cdContextDBuffer,
	#
	# From handlegraphics
	#
	guidata, setappdata, getappdata, WindowButtonDownFcn, WindowButtonMotionFcn, KeyPressFcn,
	set_current_point, get_current_point, set_current_character, get_current_character, isnada, Handles,
	gc_ref, gc_unref,
	# local
	imcdCanvasPutImage

include("libcd_h.jl")
include("libcd.jl")
include("cdiup.jl")
include("wd.jl")
include("im_image_h.jl")
include("handlegraphics.jl")

# ----------------------------------------------------------------------------------------
# Utility function to draw the image in a CD library canvas (from im_image.h)
# Works only for data_type IM_BYTE, and color spaces: IM_RGB, IM_MAP, IMGRAY and IM_BINARY.
function imcdCanvasPutImage(_canvas, _image, _x, _y, _w, _h, _xmin, _xmax, _ymin, _ymax)
	if (isa(_image, Ptr{imImage}))
		_image = unsafe_load(_image)
	end
	const IM_RGB = 0		# Redefining here to not have to import it from libim.h (Ghrrr)
	if (_image.color_space == IM_RGB)
		data = [convert(Ptr{UInt8}, unsafe_load(_image.data,1)),
			convert(Ptr{UInt8}, unsafe_load(_image.data,2)),
			convert(Ptr{UInt8}, unsafe_load(_image.data,3))]

		if (_image.has_alpha != 0)
			alfa = convert(Ptr{UInt8}, unsafe_load(_image.data,4))
			cdCanvasPutImageRectRGBA(_canvas, _image.width, _image.height, data[1], data[2],
			                         data[3], alfa, _x, _y, _w, _h, _xmin, _xmax, _ymin, _ymax)
		else
			cdCanvasPutImageRectRGB(_canvas, _image.width, _image.height, data[1], data[2],
			                        data[3], _x, _y, _w, _h, _xmin, _xmax, _ymin, _ymax)
		end
	else
		data = [convert(Ptr{UInt8}, unsafe_load(_image.data,1))]
		cdCanvasPutImageRectMap(_canvas, _image.width, _image.height, data[1], _image.palette,
		                        _x, _y, _w, _h, _xmin, _xmax, _ymin, _ymax)
	end
end

end  # module


text:
using Base.Iterators
using CVChannel
using LinearAlgebra
using QBase
using Random
using Test

"""
This script verifies that the see-saw method can reproduce the known super-multiplicativity
of the anti-symmetric Werner-Holevo channel.
We demonstrate the see-saw method's ability to scale to dimensions beyond what can be found
using the PPT relaxation on the Choi operator approach to the communication value (CV).
The verification procedure below takes about 10-20 minutes to complete.
"""

@time @testset "see-saw verification of anti-symmetric werner-holevo channel super-multiplicativiity" for d in 2:6
@time @testset "(d = $d)" begin
    Random.seed!(111) # seeding random numbers for reproducibility

    println("\nd = ", d)
    anti_sym_choi = wernerState(d, 0) * d
    nonzero_ids = filter(i -> anti_sym_choi[:,i] != zeros(d^2), 1:d^2)
    anti_sym_kraus_ops = map(i -> reshape(anti_sym_choi[:,i], (d,d)), nonzero_ids)*sqrt((d-1)/2)
    par_anti_sym_kraus_ops = collect(flatten(
        map(k1 -> map(k2 -> kron(k1,k2), anti_sym_kraus_ops), anti_sym_kraus_ops)
    ))

    @test sum(k -> k' * k, anti_sym_kraus_ops) ≈ I
    @test choi(anti_sym_kraus_ops) ≈ anti_sym_choi
    @test sum(k -> k' * k, par_anti_sym_kraus_ops) ≈ I

    # warning: this computation is slow for d=6 and larger
    @test choi(par_anti_sym_kraus_ops) ≈ parChoi(
        Choi(anti_sym_choi, d, d), Choi(anti_sym_choi, d, d)
    ).JN

    init_states = haarStates(d, d)
    max_cv_tuple, = seesawCV(init_states, anti_sym_kraus_ops, 3)

    @test max_cv_tuple[1] ≈ d/(d-1) atol=1e-6
    @test all(ρ -> is_density_matrix(ρ, atol=1e-6), max_cv_tuple[2])
    @test is_povm(max_cv_tuple[3], atol=1e-6)

    par_init_states = haarStates(d^2, d^2)
    par_max_cv_tuple, = seesawCV(
        par_init_states,
        par_anti_sym_kraus_ops,
        5,
    )

    @test par_max_cv_tuple[1] ≈ 2*d/(d-1) atol=5e-2
    @test all(ρ -> is_density_matrix(ρ, atol=1e-5), par_max_cv_tuple[2])
    @test is_povm(par_max_cv_tuple[3], atol=1e-5)

    println("max_cv : ", max_cv_tuple[1], ", match : ", d/(d-1))
    println("par_max_cv : ", par_max_cv_tuple[1], ", match : ", 2*d/(d-1), "\n")
end
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
### A Pluto.jl notebook ###
# v0.19.3

using Markdown
using InteractiveUtils

# ╔═╡ 19afaf4e-b19b-47a3-8c4c-31b8879f392d
using JSON, StanSample, Statistics, NamedTupleTools, Random

# ╔═╡ f19cee90-c255-4760-abdc-3c3da106ff9b
stan_chris = "
data {
    int n_rows;
    int<lower=1> n_cols;
    matrix<lower=0>[n_rows,n_cols] x;
}

parameters {
   real mu;
}

model {
    mu ~ normal(0, 1);
}";

# ╔═╡ 52d5a30c-37ac-4726-abb4-e4e77f9a1871
sm = SampleModel("chris", stan_chris);

# ╔═╡ 7a0e2a02-46de-403e-9c88-f0c39c86951d
begin
	seed = 65445
	Random.seed!(seed)
end;

# ╔═╡ 07e20b3a-f0e9-4842-b3eb-f8016f24331d
begin
	n_rows = 8
	n_cols = 2
	x = zeros(n_rows, n_cols)
end

# ╔═╡ 3a61ff36-66c3-4612-aef6-7c39fa8387b2
indata = (x = x, n_rows = n_rows, n_cols = n_cols)

# ╔═╡ 9fc648c2-5a17-4d4e-9f94-f03e96e39c34
size(indata.x)

# ╔═╡ a5c57d03-8bbc-41af-81b3-db4239d74037
function convert_matrices_in_nt_or_dict(d::Union{NamedTuple, Dict})
	dct = typeof(d) <: NamedTuple ? dct = convert(Dict, d) : d
	for key in keys(dct)
		if typeof(dct[key]) <: Matrix
			dct[key] = Matrix(dct[key]')
		end
	end
	dct
end

# ╔═╡ f4e2869a-686d-4b00-958d-7d5dc957de51
data = convert_matrices_in_nt_or_dict(indata)

# ╔═╡ c92e7631-2fa3-4867-b417-c7661bbcbce8
rc = stan_sample(sm; data, seed)

# ╔═╡ 43512b07-2ea2-4f23-b6da-279202cc61b5
if success(rc)
	df = read_samples(sm , :dataframe)
	[mean(df.mu), std(df.mu)]
end

# ╔═╡ e7bb012a-2f39-4ce3-818c-9976f5ada1d5
function convert_matrices(d::T) where {T <: Vector}
    dd = copy(d)
    for i in 1:length(dd)
        dct = typeof(dd[i]) <: NamedTuple ? dct = convert(Dict, dd[i]) : dd[i]
        for key in keys(dct)
            if typeof(dct[key]) <: Matrix
                dct[key] = Matrix(dct[key]')
            end
        end
        dd[i] = dct
    end
    dd
end

# ╔═╡ 80b8fde0-eeb6-4035-8ca6-042117ad9fcc
dd = [(x = x, n_rows = n_rows, n_cols = n_cols) for i in 1:4]

# ╔═╡ c26b973e-6cc9-45d6-982c-632136e049c8

# ╔═╡ 155c4b2d-b151-4212-b390-a90cfc52c579
begin
	typeof(dd[1])
	dct = convert(Dict, dd[1])
	typeof(dct)
end

# ╔═╡ 2525e4b5-abd6-4209-8fea-6d69b939d694
let
	rc = stan_sample(sm; data=dct, seed)
	if success(rc)
		df = read_samples(sm , :dataframe)
		[mean(df.mu), std(df.mu)]
	end
	[mean(df.mu), std(df.mu)]
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
StanSample = "c1514b29-d3a0-5178-b312-660c88baa699"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
JSON = "~0.21.3"
NamedTupleTools = "~0.14.0"
StanSample = "~6.4.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.9.0-DEV"
manifest_format = "2.0"
project_hash = "ff9bc0ea24be11f2223101e1a5610a25d2106397"

[[deps.ANSIColoredPrinters]]
git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
version = "0.0.1"

[[deps.ArgCheck]]
git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4"
uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197"
version = "2.3.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "c933ce606f6535a7c7b98e1d86d5d1014f730596"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "5.0.7"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.BitTwiddlingConvenienceFunctions]]
deps = ["Static"]
git-tree-sha1 = "28bbdbf0354959db89358d1d79d421ff31ef0b5e"
uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b"
version = "0.1.3"

[[deps.CPUSummary]]
deps = ["CpuId", "IfElse", "Static"]
git-tree-sha1 = "baaac45b4462b3b0be16726f38b789bf330fcb7a"
uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
version = "0.1.21"

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "873fb188a4b9d76549b81465b1f75c82aaf59238"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.4"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.43.0"

[[deps.CompatHelperLocal]]
deps = ["DocStringExtensions", "Pkg", "UUIDs"]
git-tree-sha1 = "908a50c2805d2571379fe1595585148c8352f63e"
uuid = "5224ae11-6099-4aaa-941d-3aab004bd678"
version = "0.1.20"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.5.2+0"

[[deps.CpuId]]
deps = ["Markdown"]
git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406"
uuid = "adafc99b-e345-5852-983c-f28acb93d879"
version = "0.3.1"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "fb5f5316dd3fd4c5e7c30a24d50643b73e37cd40"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.10.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "6c19003824cbebd804a51211fd3bbd81bf1ecad5"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.3"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "70f5bfdfbdc6c9d2b7a143d70ae88f4cb7b193b1"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.56"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Documenter]]
deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
git-tree-sha1 = "6edbf28671b4df4f692e54ae72f1e35851cfbf38"
uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
version = "0.27.16"

[[deps.Downloads]]
deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
version = "1.6.0"

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "129b104185df66e408edd6625d480b7f9e9823a0"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.18"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.HostCPUFeatures]]
deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
git-tree-sha1 = "18be5268cf415b5e27f34980ed25a7d34261aa83"
uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
version = "0.1.7"

[[deps.Hwloc]]
deps = ["Hwloc_jll"]
git-tree-sha1 = "92d99146066c5c6888d5a3abc871e6a214388b91"
uuid = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
version = "2.0.0"

[[deps.Hwloc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "303d70c961317c4c20fafaf5dbe0e6d610c38542"
uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
version = "2.7.1+0"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "61feba885fac3a407465726d0c330b3055df897f"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LayoutPointers]]
deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static"]
git-tree-sha1 = "b651f573812d6c36c22c944dd66ef3ab2283dfa1"
uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
version = "0.1.6"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
version = "0.6.3"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
version = "7.81.0+0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
version = "1.10.2+0"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "76c987446e8d555677f064aaac1145c4c17662f8"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.14"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.ManualMemory]]
git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667"
version = "0.1.8"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
version = "2.28.0+0"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MonteCarloMeasurements]]
deps = ["Distributed", "Distributions", "LinearAlgebra", "MacroTools", "Random", "RecipesBase", "Requires", "SLEEFPirates", "StaticArrays", "Statistics", "StatsBase", "Test"]
git-tree-sha1 = "03619e255664666b352a5e5f6b45e8b00d439870"
uuid = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
version = "1.0.8"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
version = "2022.2.1"

[[deps.NamedTupleTools]]
git-tree-sha1 = "befc30261949849408ac945a1ebb9fa5ec5e1fd5"
uuid = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
version = "0.14.0"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
version = "1.2.0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
version = "0.3.20+0"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
version = "0.8.1+0"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "3114946c67ef9925204cc024a73c9e679cebe0d7"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.8"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "1285416549ccfcdf0c50d4997a94331e88d68413"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.3.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
version = "1.8.0"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.2"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.3.0"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
version = "0.7.0"

[[deps.SIMDTypes]]
git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c"
uuid = "94e857df-77ce-4151-89e5-788b33177be4"
version = "0.1.0"

[[deps.SLEEFPirates]]
deps = ["IfElse", "Static", "VectorizationBase"]
git-tree-sha1 = "ac399b5b163b9140f9c310dfe9e9aaa225617ff6"
uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
version = "0.6.32"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "6a2f7d70512d205ca8c7ee31bfa9f142fe74310c"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.12"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.StanBase]]
deps = ["CSV", "DataFrames", "DelimitedFiles", "Distributed", "DocStringExtensions", "Documenter", "JSON", "Parameters", "Random", "StanDump", "Unicode"]
git-tree-sha1 = "fe0c4b9057c1e23e1f8e27666e63f71eedf7b023"
uuid = "d0ee94f6-a23d-54aa-bbe9-7f572d6da7f5"
version = "4.4.0"

[[deps.StanDump]]
deps = ["ArgCheck", "DocStringExtensions"]
git-tree-sha1 = "bfaebe19ada44a52a6c797d48473f1bb22fd0853"
uuid = "9713c8f3-0168-54b5-986e-22c526958f39"
version = "0.2.0"

[[deps.StanSample]]
deps = ["CSV", "CompatHelperLocal", "DataFrames", "DelimitedFiles", "Distributed", "DocStringExtensions", "JSON", "MonteCarloMeasurements", "NamedTupleTools", "OrderedCollections", "Parameters", "Random", "Reexport", "Requires", "StanBase", "StanDump", "TableOperations", "Tables", "Unicode"]
git-tree-sha1 = "7301a753d5e78694b0a9bd4fe7df6c366ef500ec"
uuid = "c1514b29-d3a0-5178-b312-660c88baa699"
version = "6.4.0"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "91181e5820a400d1171db4382aa36e7fd19bee27"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.3"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "cd56bf18ed715e8b09f06ef8c6b781e6cdc49911"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.4"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c82aaa13b44ea00134f8c9c89819477bd3986ecd"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.3.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "5950925ff997ed6fb3e985dcce8eb1ba42a0bbe7"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.18"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
version = "1.0.0"

[[deps.TableOperations]]
deps = ["SentinelArrays", "Tables", "Test"]
git-tree-sha1 = "e383c87cf2a1dc41fa30c093b2a19877c83e1bc1"
uuid = "ab02a1b2-a7df-11e8-156e-fb1833f50b87"
version = "1.2.0"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
version = "1.10.0"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.VectorizationBase]]
deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "Hwloc", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static"]
git-tree-sha1 = "858e541ffc21873e45aeaf744e0d015966e0328e"
uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
version = "0.21.30"

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.2"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
version = "1.2.12+3"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
version = "5.1.0+0"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
version = "1.41.0+1"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
version = "16.2.1+1"
"""

# ╔═╡ Cell order:
# ╠═19afaf4e-b19b-47a3-8c4c-31b8879f392d
# ╠═f19cee90-c255-4760-abdc-3c3da106ff9b
# ╠═52d5a30c-37ac-4726-abb4-e4e77f9a1871
# ╠═7a0e2a02-46de-403e-9c88-f0c39c86951d
# ╠═07e20b3a-f0e9-4842-b3eb-f8016f24331d
# ╠═3a61ff36-66c3-4612-aef6-7c39fa8387b2
# ╠═9fc648c2-5a17-4d4e-9f94-f03e96e39c34
# ╠═a5c57d03-8bbc-41af-81b3-db4239d74037
# ╠═f4e2869a-686d-4b00-958d-7d5dc957de51
# ╠═c92e7631-2fa3-4867-b417-c7661bbcbce8
# ╠═43512b07-2ea2-4f23-b6da-279202cc61b5
# ╠═e7bb012a-2f39-4ce3-818c-9976f5ada1d5
# ╠═80b8fde0-eeb6-4035-8ca6-042117ad9fcc
# ╠═c26b973e-6cc9-45d6-982c-632136e049c8
# ╠═155c4b2d-b151-4212-b390-a90cfc52c579
# ╠═2525e4b5-abd6-4209-8fea-6d69b939d694
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using CartesianGP
using Base.Test

# Utility functions

function default_parameters()
    numperlevel = 2
    numlevels = 1
    numlevelsback = 1
    numinputs = 2
    numoutputs = 2

    p = Parameters(numinputs, numoutputs, numperlevel, numlevels, numlevelsback)

    return p
end

function chromosome_half_adder(numperlevel=2, numlevels=1, numlevelsback=1)
    numinputs = 2
    numoutputs = 2

    p = Parameters(numinputs, numoutputs, numperlevel, numlevels, numlevelsback)
    c = Chromosome(p)

    c.inputs = [InputNode(1), InputNode(2)]
    c.interiors = [InteriorNode(XOR, [(0, 1), (0, 2)]) InteriorNode(AND, [(0, 1), (0, 2)])]
    c.outputs = [OutputNode((1, 1)), OutputNode((1, 2))]

    return c
end

# 1100
# 1010

function goal_half_adder()
    g = Goal(2, (0b0110, 0b1000))

    return g
end

include("Func.jl")
include("Chromosome.jl")
include("Goal.jl")
include("Fitness.jl")
include("Evolution.jl")

numinputs = 2
numoutputs = 1
numperlevel = 1
numlevels = 10
numlevelsback = 10

funcs = default_funcs()
p = Parameters(numinputs, numoutputs, numperlevel, numlevels, numlevelsback)

for _ = 1:100
    c0 = random_chromosome(p)
    # Executing on these inputs tests all possible bit combinations for inputs
    execute_chromosome(c0, [convert(BitString, 0xC), convert(BitString, 0xA)])

    # Try a mutation, right now we just want to make sure it doesn't crash
    c1 = mutate(c0, funcs)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Yosys_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Yosys")
JLLWrappers.@generate_main_file("Yosys", UUID("2eded7df-79ee-58af-89be-7676cda93c14"))
end  # module Yosys_jll


text:
# Hessenberg factorization

if VERSION < v"1.3"
    struct HessenbergFactorization{T, S<:StridedMatrix, U} <: Factorization{T}
        data::S
        τ::Vector{U}
    end
    Base.size(H::HessenbergFactorization, args...) = size(H.data, args...)
end

if VERSION < v"1.3"

    function _hessenberg!(A::StridedMatrix{T}) where T
        n = LinearAlgebra.checksquare(A)
        τ = Vector{Householder{T}}(undef, n - 1)
        for i = 1:n - 1
            xi = view(A, i + 1:n, i)
            t  = LinearAlgebra.reflector!(xi)
            H  = Householder{T,typeof(xi)}(view(xi, 2:n - i), t)
            τ[i] = H
            lmul!(H', view(A, i + 1:n, i + 1:n))
            rmul!(view(A, :, i + 1:n), H)
        end
        return HessenbergFactorization{T, typeof(A), eltype(τ)}(A, τ)
    end

    function _materializeQ(H::HessenbergFactorization{T}) where T
        n = size(H.data, 1)
        Z = Matrix{T}(I, n, n)
        for j=n-1:-1:1
            lmul!(H.τ[j], view(Z, j+1:n, j:n))
            Z[1:j-1,j] .= 0
        end
        Z
    end
else # version v1.3+

function _hessenberg!(A::StridedMatrix{T}) where T
    n = LinearAlgebra.checksquare(A)
    τ = Vector{T}(undef, n - 1)
    for i = 1:n - 1
        xi = view(A, i + 1:n, i)
        t  = LinearAlgebra.reflector!(xi)
        H  = Householder{T,typeof(xi)}(view(xi, 2:n - i), t)
        τ[i] = H.τ
        lmul!(H', view(A, i + 1:n, i + 1:n))
        rmul!(view(A, :, i + 1:n), H)
    end
    return Hessenberg(A, τ)
end
using LinearAlgebra: QRPackedQ
function _materializeQ(H::Hessenberg{T}) where {T}
    A = copy(H.Q.factors)
    n = checksquare(A)
    # shift reflectors one column rightwards
    @inbounds for j in n:-1:2
        A[1:j-1,j] .= zero(T)
        for i in j+1:n
            A[i,j] = A[i,j-1]
        end
    end
    A[2:n,1] .= zero(T)
    A[1,1] = one(T)
    Q1 = QRPackedQ(view(A,2:n,2:n), H.Q.τ)
    A[2:n,2:n] .= Matrix{T}(Q1)
    A
end

end # v1.3+ branch

LinearAlgebra.hessenberg!(A::StridedMatrix) = _hessenberg!(A)


text:
using DataStructures
import Base.convert, Base.median

export Slider, slider, ToggleButton, togglebutton, Button, button,
       Options, Checkbox, checkbox, Textbox, textbox, Textarea, textarea,
       RadioButtons, radiobuttons, Dropdown, dropdown, Select, select,
       ToggleButtons, togglebuttons, HTML, html, Latex, latex,
       Progress, progress, widget

median(r::Range) = r[(1+length(r))>>1]

### Input widgets

########################## Slider ############################

type Slider{T<:Number} <: InputWidget{T}
    signal::Input{T}
    label::String
    value::T
    range::Range{T}
end

slider(args...) = Sldier(args...)
slider{T}(range::Range{T};
          value=median(range),
          signal::Signal{T}=Input(value),
          label="") =
              Slider(signal, label, value, range)

######################### Checkbox ###########################

type Checkbox <: InputWidget{Bool}
    signal::Input{Bool}
    label::String
    value::Bool
end

checkbox(args...) = Checkbox(args...)
checkbox(value::Bool; signal=Input(value), label="") =
    Checkbox(signal, label, value)
checkbox(; label="", value=false, signal=Input(value)) =
    Checkbox(signal, label, value)

###################### ToggleButton ########################

type ToggleButton <: InputWidget{Bool}
    signal::Input{Bool}
    label::String
    value::Bool
end

togglebutton(args...) = ToggleButton(args...)

togglebutton(; label="", value=false, signal=Input(value)) =
    ToggleButton(signal, label, value)

togglebutton(label; kwargs...) =
    togglebutton(label=label; kwargs...)

######################### Button ###########################

type Button{T} <: InputWidget{T}
    signal::Input{T}
    label::String
    value::T
end

button(; value=nothing, label="", signal=Input(value)) =
    Button(signal, label, value)

button(label; kwargs...) =
    button(label=label; kwargs...)

######################## Textbox ###########################

type Textbox{T <: Union(Number, String)} <: InputWidget{T}
    signal::Input{T}
    label::String
    range::Union(Nothing, Range)
    value::T
end

function empty(t::Type)
    if is(t, Number) zero(t)
    elseif is(t, String) ""
    end
end

function Textbox(; typ=String, label="",
                 value=empty(typ),
                 range=nothing,
                 signal=Input(value))
    if isa(value, String) && !isa(range, Nothing)
        throw(ArgumentError(
               "You cannot set a range on a string textbox"
             ))
    end
    Textbox(signal, label, range, value)
end

textbox(;kwargs...) = Textbox(;kwargs...)
textbox(val; kwargs...) =
    Textbox(value=val; kwargs...)

function parse{T<:Number}(val, w::Textbox{T})
    v = convert(T, val)
    if isa(w.range, Range)
        # force value to stay in range
        v = max(first(w.range),
                min(last(w.range), v))
    end
    v
end

######################### Textarea ###########################

type Textarea{String} <: InputWidget{String}
    signal::Input{String}
    label::String
    value::String
end

textarea(args...) = Textarea(args...)

textarea(; label="",
         value="",
         signal=Input(value)) =
    Textarea(signal, label, value)

textarea(val; kwargs...) =
    textarea(value=val; kwargs...)

##################### SelectionWidgets ######################

type Options{view, T} <: InputWidget{T}
    signal::Input{T}
    label::String
    value::T
    value_label::String
    options::OrderedDict{String, T}
    # TODO: existential checks
end

Options{T}(view::Symbol, options::OrderedDict{String, T};
        label = "",
        value_label=first(options)[1],
        value=options[value_label],
        signal=Input(value)) =
            Options{view, T}(signal, label, value, value_label, options)

function Options{T}(view::Symbol,
                    options::AbstractArray{T};
                    kwargs...)
    opts = OrderedDict{String, T}()
    map(v -> opts[string(v)] = v, options)
    Options(view, opts; kwargs...)
end

function Options{K, V}(view::Symbol,
                    options::Associative{K, V};
                    kwargs...)
    opts = OrderedDict{String, V}()
    map(v->opts[string(v[1])] = v[2], options)
    Options(view, opts; kwargs...)
end

dropdown(opts; kwargs...) =
    Options(:Dropdown, opts; kwargs...)

radiobuttons(opts; kwargs...) =
    Options(:RadioButtons, opts; kwargs...)

select(opts; kwargs...) =
    Options(:Select, opts; kwargs...)

togglebuttons(opts; kwargs...) =
    Options(:ToggleButtons, opts; kwargs...)

### Output Widgets

export HTML, Latex, Progress

type HTML <: Widget
    label::String
    value::String
end
html(label, value) = HTML(label, value)
html(value; label="") = HTML(label, value)

# assume we already have HTML
## writemime(io::IO, m::MIME{symbol("text/html")}, h::HTML) =
##     write(io, h.value)

type Latex <: Widget
    label::String
    value::String
end
latex(label, value) = Latex(label, value)
latex(value; label="") = Latex(label, value)

## # assume we already have Latex
## writemime(io::IO, m::MIME{symbol("application/x-latex")}, l::Latex) =
##     write(io, l.value)

type Progress <: Widget
    label::String
    value::Int
    range::Range
end

progress(args...) = Progress(args...)
progress(;label="", value=0, range=0:100) =
    Progress(label, value, range)

# Make a widget out of a domain
widget(x::Signal, label="") = x
widget(x::Widget, label="") = x
widget(x::Range, label="") = slider(x, label=label)
widget(x::AbstractVector, label="") = togglebuttons(x, label=label)
widget(x::Associative, label="") = togglebuttons(x, label=label)
widget(x::Bool, label="") = checkbox(x, label=label)
widget(x::String, label="") = textbox(x, label=label)
widget{T <: Number}(x::T, label="") = textbox(typ=T, value=x, label=label)


text:
struct MInfo{T}
    idx_q1::Vector{Int}
    idx_q2::Vector{Int}
    idx_u1::Vector{Int}
    v1::Vector{T}
end

function _step(sim::Simulator, info::MInfo, q1, q2, u1)
    info.v1 .= q2
    info.v1 .-= q1
    info.v1 ./= sim.h
    RoboDojo.step!(sim, q2, info.v1, u1, 1)
end

struct GradientBundle{T}
    ls::LeastSquares{T}
    q1η::Vector{T}
    q2η::Vector{T}
    u1η::Vector{T}
    dz::Matrix{T}
    ny::Int
    nz::Int
    info::MInfo{T}
end

function GradientBundle(model; N=100, ϵ=1.0e-4)
    nx = 2 * model.nq
    ny = model.nq
    nu = model.nu
    nz = nx + nu
    nθ = ny * (nz)

    @variables fz[1:ny] fη[1:ny] η[1:nz] θ[1:nθ]

    function cost(fz, fη, η, θ)
        M = reshape(θ, ny, nz)
        r = fη - fz - M * η
        return [transpose(r) * r]
    end

    c = cost(fz, fη, η, θ)
    cθ = Symbolics.gradient(c[1], θ)
    cθθ = Symbolics.hessian(c[1], θ)

    c_func_ls = eval(Symbolics.build_function(c, fz, fη, η, θ)[2])
    cθ_func_ls = eval(Symbolics.build_function(cθ, fz, fη, η, θ)[2])
    cθθ_func_ls = eval(Symbolics.build_function(cθθ, fz, fη, η, θ)[2])

    η = [zeros(nz) for i = 1:N]

    for i = 1:N
        w = ϵ * randn(1)[1]
        η[i][rand([j for j = 1:nz], 1)[1]] = w
    end

    # for i = 1:nz
    #     η[i][i] = ϵ
    #     η[i+nz][i] = -ϵ
    # end

    info = MInfo(
        collect(1:model.nq),
        collect(model.nq .+ (1:model.nq)),
        collect(2 * model.nq .+ (1:model.nu)),
        zeros(model.nq))

    fz = zeros(model.nq)
    fη = [zeros(model.nq) for i = 1:N]

    ls = LeastSquares(N,
        fz, fη, η,
        c_func_ls, cθ_func_ls, cθθ_func_ls,
        [0.0], zeros(nθ), zeros(nθ, nθ),
        [0.0], zeros(nθ), zeros(nθ, nθ),
        zeros(nθ), zeros(nθ),
        lu_solver(zeros(nθ, nθ)))

    GradientBundle(ls,
        zeros(nq),
        zeros(nq),
        zeros(nu),
        zeros(ny, nz),
        ny, nz,
        info)
end

function gradient!(sim::Simulator, gb::GradientBundle, q1, q2, u1; verbose=false)
    gb.ls.fz .= _step(sim, gb.info, q1, q2, u1)
    for i = 1:gb.ls.N
        ηq1 = @views gb.ls.η[i][gb.info.idx_q1]
        ηq2 = @views gb.ls.η[i][gb.info.idx_q2]
        ηu1 = @views gb.ls.η[i][gb.info.idx_u1]
        gb.q1η .= q1
        gb.q1η .+= ηq1
        gb.q2η .= q2
        gb.q2η .+= ηq2
        gb.u1η .= u1
        gb.u1η .+= ηu1
        gb.ls.fη[i] .= _step(sim, gb.info, gb.q1η, gb.q2η, gb.u1η)
    end
    update!(gb.ls, verbose=verbose)
    gb.dz .= Base.ReshapedArray(gb.ls.θ, (gb.ny, gb.nz), ())
    return gb.dz
end

# gradient!(im_dyn.eval_sim, gb, q1, q2, u1)
# @benchmark gradient!($im_dyn.eval_sim, $gb, $q1, $q2, $u1)

function fx_gb(dx, model::ImplicitDynamics, x, u, w)
	q1 = @views x[model.idx_q1]
	q2 = @views x[model.idx_q2]
	model.v1 .= q2
	model.v1 .-= q1
	model.v1 ./= model.eval_sim.h
	RoboDojo.step!(model.eval_sim, q2, model.v1, u, 1)
	nq = model.eval_sim.model.nq
	for i = 1:nq
		dx[model.idx_q1[i], model.idx_q2[i]] = 1.0
	end
    gradient!(model.eval_sim, model.info, q1, q2, u)
    ∂q3∂q1 = @views model.info.dz[model.idx_q1, model.idx_q1]
    ∂q3∂q2 = @views model.info.dz[model.idx_q1, model.idx_q2]
	dx[model.idx_q2, model.idx_q1] = ∂q3∂q1
	dx[model.idx_q2, model.idx_q2] = ∂q3∂q2
	return dx
end

# dx = zeros(2 * planarpush.nq, 2 * planarpush.nq)
# x = rand(2 * planarpush.nq)
# u = rand(planarpush.nu)
# w = zeros(0)
# fx_gb(dx, im_dyn, x, u, w)
# @benchmark fx_gb($dx, $im_dyn, $x, $u, $w)
# @code_warntype fx_gb(dx, im_dyn, x, u, w)

function fu_gb(du, model::ImplicitDynamics, x, u, w)
	q1 = @views x[model.idx_q1]
	q2 = @views x[model.idx_q2]
	model.v1 .= q2
	model.v1 .-= q1
	model.v1 ./= model.eval_sim.h
	RoboDojo.step!(model.eval_sim, q2, model.v1, u, 1)
    gradient!(model.eval_sim, model.info, q1, q2, u)
    ∂q3∂u1 = @views model.info.dz[model.idx_q1, model.info.info.idx_u1]
	du[model.idx_q2, :] = ∂q3∂u1
	return du
end

# du = zeros(2 * planarpush.nq, planarpush.nu)
# x = rand(2 * planarpush.nq)
# u = rand(planarpush.nu)
# w = zeros(0)
# fu_gb(du, im_dyn, x, u, w)
# @benchmark fu_gb($du, $im_dyn, $x, $u, $w)
# @code_warntype fu_gb(du, im_dyn, x, u, w)


text:
using LightGraphs: SimpleGraph

@testset "Chordal Graph tests" begin
    # From Garstka paper
    n1 = 9
    nonzero_inds = vcat([
            (3,1), (3,2), (5,4), (6,1), (6,3), (7,3), (7,6), (8,3), (8,4),
            (8,5), (8,6), (8,7), (9,6), (9,7), (9,8)
        ],
        [(i,i) for i in 1:n1]
    )
    sp = sparse(CD.unzip(nonzero_inds)..., [1.0i+j for (i, j) in nonzero_inds])
    A1 = sp + tril(sp, -1)'

    # From VA Figure 4.2
    n2 = 17
    ijs = [(3,1), (3,2), (4,1), (4,2), (4,3), (5,1), (5,3), (5,4), (8,7),
           (9,5), (9,6), (9,7), (9,8), (11,10), (13,10), (13,11), (13,12),
           (14,10), (14,11), (14,12), (14,13),
           (15, 1), (15,3), (15,4), (15,5), (15,7), (15,8), (15,9),
           (16,5), (16,6), (16,9), (16,12), (16,13), (16,14), (16,15),
           (17,10), (17,11), (17,12), (17,13), (17,14), (17,15), (17,16)]
    append!(ijs, [(i,i) for i in 1:n2])
    II, JJ = Chordal.unzip(ijs)
    A2 = sparse(II, JJ, ones(length(II)))
    A2 = (A2 + A2') / 2

    # TODO: Figure out error from LightGraphs
    peo = CD.maximum_cardinality_search(A1)
    # @test peo == CD.maximum_cardinality_search(SimpleGraph(A1))
    @test CD.is_chordal(A1; peo=peo)
    # @test CD.is_chordal(SimpleGraph(A1); peo=1:n1)

    peo = CD.maximum_cardinality_search(A2)
    # @test peo == CD.maximum_cardinality_search(SimpleGraph(A2))
    @test CD.is_chordal(A2; peo=peo)
    # @test CD.is_chordal(SimpleGraph(A2); peo=1:n2)
end


text:
@testset "CUDA: Galileo E1B" begin
    galileo_e1b = GalileoE1B(use_gpu = Val(true))
    @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 29)) == true

    @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 6)) == false

    @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xc, 40)) == false

    @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf, 8)) == true

    @test @inferred(Tracking.is_upcoming_integration_new_bit(galileo_e1b, 0xf0, 8)) == true

    @test @inferred(Tracking.get_default_correlator(galileo_e1b, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1))
    @test @inferred(Tracking.get_default_correlator(galileo_e1b, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3))
end


text:
# ------------ Constraints to direct addresses ------------ #

struct ConstrainedByAddress <: ConstrainedSelectQuery
    query::Dict{Any, Any}
    ConstrainedByAddress() = new(Dict())
    ConstrainedByAddress(d::Dict) = new(d)
end

has_top(csa::ConstrainedByAddress, addr) = haskey(csa.query, addr)

function dump_queries(csa::ConstrainedByAddress)
    v = Vector{Pair}()
    for (k, l) in csa.query
        push!(v, (k, ) => l)
    end
    return v
end

get_top(csa::ConstrainedByAddress, addr) = getindex(csa.query, addr)

isempty(csa::ConstrainedByAddress) = isempty(csa.query)

push!(sel::ConstrainedByAddress, addr, val) = sel.query[addr] = val

merge!(sel1::ConstrainedByAddress, sel2::ConstrainedByAddress) = Base.merge!(sel1.query, sel2.query)

addresses(csa::ConstrainedByAddress) = keys(csa.query)

function compare(chs::ConstrainedByAddress, v::Visitor)
    addrs = []
    for addr in addresses(chs)
        addr in v.addrs && continue
        push!(addrs, addr)
    end
    return isempty(addrs), addrs
end

function ==(cba1::ConstrainedByAddress, cba2::ConstrainedByAddress)
    for (k, v) in cba1.query
        k in keys(cba2.query) || return false
        cba2.query[k] == v || return false
    end
    return true
end

# Functional filter.
function filter(k_fn::Function, v_fn::Function, query::ConstrainedByAddress) where T <: Address
    top = ConstrainedByAddress()
    for (k, v) in query.query
        k_fn(k) && v_fn(v) && push!(top, k, v)
    end
    return top
end

# ------------ Utility for pretty printing ------------ #

function collect!(par, addrs, chd, query::ConstrainedByAddress)
    for (k, v) in query.query
        push!(addrs, (par..., k))
        chd[(par..., k)] = v
    end
end

function collect!(addrs, chd, query::ConstrainedByAddress)
    for (k, v) in query.query
        push!(addrs, (k, ))
        chd[(k, )] = v
    end
end

# ----------- Selection to direct addresses ------------ #

struct UnconstrainedByAddress <: UnconstrainedSelectQuery
    query::Vector{Any}
    UnconstrainedByAddress() = new(Any[])
end

has_top(csa::UnconstrainedByAddress, addr) = addr in csa.query

dump_queries(csa::UnconstrainedByAddress) = csa.query

isempty(csa::UnconstrainedByAddress) = isempty(csa.query)

push!(sel::UnconstrainedByAddress, addr::T) where T <: Address = push!(sel.query, addr)

merge!(sel1::UnconstrainedByAddress, sel2::UnconstrainedByAddress) = union(sel1.query, sel2.query)

addresses(usa::UnconstrainedByAddress) = usa.query

function compare(uba::UnconstrainedByAddress, v::Visitor)
    addrs = []
    for addr in addresses(uba)
        addr in v.addrs && continue
        push!(addrs, addr)
    end
    return isempty(addrs), addrs
end

function ==(cba1::UnconstrainedByAddress, cba2::UnconstrainedByAddress)
    cba1.query == cba2.query
end

# Functional filter.
function filter(k_fn::Function, v_fn::Function, query::UnconstrainedByAddress) where T <: Address
    top = UnconstrainedByAddress()
    for k in query.query
        k_fn(k) && push!(top, k)
    end
    return top
end

# ------------ Utility for pretty printing ------------ #

function collect!(par, addrs, query::UnconstrainedByAddress)
    for k in query.query
        push!(addrs, (par..., k))
    end
end

function collect!(addrs, query::UnconstrainedByAddress)
    for k in query.query
        push!(addrs, (k, ))
    end
end


text:
using RobotDynamics: state, control, states, controls

@testset "Trajectories" begin
n,m,N = 5,3,11
t,dt = 0.0,0.1
x = @SVector rand(n)
u = @SVector rand(m)
z = RD.KnotPoint(x,u,t,dt)

#--- Empty contructor
Z = RD.Traj(n,m,dt,N)
@test all(isnan, state(Z[rand(1:N)]))
@test control(Z[rand(1:N-1)]) == zeros(m)
@test length(controls(Z)) == N-1
@test length(states(Z)) == N
@test eltype(states(Z)) <: SVector{n}
@test eltype(controls(Z)) <: SVector{m}
@test Z[1].dt ≈ dt
@test Z[end].dt ≈ 0
@test Z[1].t ≈ 0
@test Z[N].t ≈ (N-1)*dt
@test RobotDynamics.is_terminal(Z[end]) == true
@test RobotDynamics.traj_size(Z) == (n,m,N)
@test RobotDynamics.num_vars(Z) == N*(n+m)-m

Z = RD.Traj(n,m,dt,N, equal=true)
@test all(isnan, state(Z[rand(1:N)]))
@test control(Z[rand(1:N)]) == zeros(m)
@test length(controls(Z)) == N
@test length(states(Z)) == N
@test eltype(states(Z)) <: SVector{n}
@test eltype(controls(Z)) <: SVector{m}
@test Z[1].dt ≈ dt
@test Z[end].dt ≈ dt
@test Z[1].t ≈ 0
@test Z[N].t ≈ (N-1)*dt
@test RobotDynamics.is_terminal(Z[end]) == false
@test RobotDynamics.num_vars(Z) == N*(n+m)

#--- Copy single constructor
Z = RD.Traj(x, u, dt, N)
@test length(Z) == N
@test eltype(Z) <: RD.KnotPoint{n,m,<:Any,Float64}
@test state(Z[1]) == x
Z[1].z = 2*Z[1].z
@test state(Z[1]) ≈ 2x
@test state(Z[2]) ≈ x
@test control(Z[1]) ≈ 2u
@test control(Z[2]) ≈ u
@test Z[1].dt ≈ dt
@test Z[end].dt ≈ 0
@test Z[1].t ≈ 0
@test Z[N].t ≈ (N-1)*dt

#--- Vector constructor
X = [@SVector rand(n) for k = 1:N]
U = [@SVector rand(m) for k = 1:N-1]
Z = RD.Traj(X,U,fill(dt,N))
@test states(Z) ≈ X
@test controls(Z) ≈ U
@test RobotDynamics.gettimes(Z) ≈ range(0, length=N, step=dt)

Z2 = copy(Z)
RobotDynamics.setstate!(Z[1], x)
@test !(state(Z[1]) ≈ state(Z2[1]))
@test state(Z[2]) ≈ state(Z2[2])
X = 2 .* X
U = 3 .* U
X[1] = x
RobotDynamics.setstates!(Z2, X)
@test state(Z[1]) ≈ state(Z2[1])
@test !(control(Z2[2]) ≈ U[2])
RobotDynamics.setcontrols!(Z2, U)
@test control(Z2[2]) ≈ U[2]

X2 = hcat(3X...)
RobotDynamics.setstates!(Z2, X2)
state(Z2[4]) ≈ 3X[4]
U2 = hcat(2U...)
RobotDynamics.setcontrols!(Z2, U2)
@test control(Z2[2]) ≈ 2U[2]

RobotDynamics.setcontrols!(Z2, u)
@test control(Z2[3]) ≈ u

times = range(0,length=N,step=dt*2)
@test RobotDynamics.gettimes(Z) ≈ times ./2
RobotDynamics.settimes!(Z, times)
@test Z[end].t ≈ 2
@test Z[1].dt ≈ 2dt

# Test approx
RobotDynamics.settimes!(Z, times ./ 2)
@test !(Z ≈ Z2)
copyto!(Z, Z2)
@test Z[1] ≈ Z2[1]
@test Z ≈ Z2

@test RobotDynamics.is_terminal(Z[end]) == true
push!(U, 2*U[end])
Z = RD.Traj(X,U,fill(dt,N))
@test length(controls(Z)) == N
@test RobotDynamics.is_terminal(Z[end]) == false

#--- Test copyto!
Z0 = [RD.KnotPoint(3*X[k], 2*U[k], dt, dt*(k-1)) for k = 1:N]
@test state.(Z0) ≈ 3 .* states(Z)
@test control.(Z0) ≈ 2 .* controls(Z)
copyto!(Z0, Z)
@test state.(Z0) ≈ states(Z)
@test control.(Z0) ≈ controls(Z)

Z2 = RD.Traj(rand() .* X, rand() .* U,fill(dt,N))
@test !(states(Z) ≈ states(Z2))
@test !(controls(Z) ≈ controls(Z2))
copyto!(Z2, Z)
@test states(Z) ≈ states(Z2)
@test controls(Z) ≈ controls(Z2)

# Test iteration
Z = RD.Traj(X,U,fill(dt,N))
@test Z[1] ≈ RD.KnotPoint(X[1],U[1],0.0,dt)
@test Z[end] ≈ RD.KnotPoint(X[end], U[end], dt, dt*(N-1))
@test Base.IndexStyle(Z) == IndexLinear()
Z_ = [z for z in Z]
@test Z_[1] === Z[1]
@test Base.IteratorSize(Z) == Base.HasLength()
@test Base.IteratorEltype(Z) == Base.HasEltype()

@test states(Z, 2) == [x[2] for x in X]
@test states(Z, 1:3) == [[x[i] for x in X] for i = 1:3]

#--- Test functions on trajectories
model = Cartpole()
dmodel = RD.DiscretizedDynamics{RD.RK4}(model)
n,m = size(model)
fVal = [@SVector zeros(n) for k = 1:N]
X = [@SVector rand(n) for k = 1:N]
U = [@SVector rand(m) for k = 1:N-1]
Z = RD.Traj(X,U,fill(dt,N))

# Test shift fill
Z = RD.Traj(X,U,fill(dt,N))
RobotDynamics.shift_fill!(Z)
@test state(Z[1]) ≈ X[2]
@test control(Z[2]) ≈ U[3]
@test Z[1].t ≈ dt
@test Z[end].t ≈ N*dt
@test Z[N].dt == 0
@test Z[N-1].dt ≈ dt

# # Test dynamics evaluation
# discrete_dynamics!(RK3, fVal, model, Z)
# @test fVal[1] ≈ discrete_dynamics(RK3, model, X[1], U[1], 0.0, dt)
# dyn = TO.DynamicsConstraint{RK3}(model, N)
# conval = TO.ConVal(n,m, dyn, 1:N-1)
# TO.evaluate!(conval, Z)
# @test conval.vals ≈ fVal[1:N-1] .- X[2:N]
# @test !(conval.vals ≈ [zeros(n) for k = 1:N-1])

# Test rollout
RD.rollout!(RD.StaticReturn(), dmodel, Z, X[1])
Zmut = RD.Traj([RD.KnotPoint{n,m}(MVector(z.z), z.t, z.dt) for z in Z])
RD.rollout!(RD.InPlace(), dmodel, Zmut, X[1])
@test Z ≈ Zmut
end


text:
const BitSSigned32{V}      = Union{SInt8{V},SInt16{V},SInt32{V}}
const BitSUnsigned32{V}    = Union{SUInt8{V},SUInt16{V},SUInt32{V}}
const BitSInteger32{V}     = Union{BitSSigned32{V},BitSUnsigned32{V}}

const BitSSigned64{V}      = Union{BitSSigned32{V},SInt{V}}
const BitSUnsigned64{V}    = Union{BitSUnsigned32{V},SUInt{V}}
const BitSInteger64{V}     = Union{BitSSigned64{V},BitSUnsigned64{V}}

const BitSSigned{V}        = Union{BitSSigned64{V},SInt128{V}}
const BitSUnsigned{V}      = Union{BitSUnsigned64{V},SUInt128{V}}
const BitSInteger{V}       = Union{BitSSigned{V},BitSUnsigned{V}}

if Int === Int64
    const BitSSignedSmall{V} = Union{SInt8{V},SInt16{V},SInt32{V}}
else
    const BitSSignedSmall{V} = Union{SInt8{V},SInt16{V}}
end

if UInt === UInt64
    const BitSUnsignedSmall{V} = Union{SUInt8{V},SUInt16{V},SUInt32{V}}
else
    const BitSUnsignedSmall{V} = Union{SUInt8{V},SUInt16{V}}
end
const BitSIntegerSmall{V}  = Union{BitSSignedSmall{V},BitSUnsignedSmall{V}}

const BitSSigned64T{V} = Union{Type{SInt8{V}},Type{SInt16{V}},Type{SInt32{V}},Type{SInt64{V}}}
const BitSUnsigned64T{V} = Union{Type{SUInt8{V}},Type{SUInt16{V}},Type{SUInt32{V}},Type{SUInt64{V}}}
const BitIntegerT{V} = Union{Type{SInt8{V}},Type{SInt16{V}},Type{SInt32{V}},Type{SInt64{V}},
                             Type{SUInt8{V}},Type{SUInt16{V}},Type{SUInt32{V}},Type{SUInt64{V}}}

    >>(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{>>(V1::$BT, V2::$BT)}()
    <<(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{<<(V1::$BT, V2::$BT)}()
    >>>(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{>>>(V1::$BT, V2::$BT)}()

for ST in static_integer
    @eval begin

        (*)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{(*)(V1::$BT, V2::$BT)}()
        (+)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{(+)(V1::$BT, V2::$BT)}()
        (-)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{(-)(V1::$BT, V2::$BT)}()
        (-)(::$ST{V}) where V = $ST{-V::$BT}()
        (|)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{(|)(V1::$BT, V2::$BT)}()
        (&)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{(&)(V1::$BT, V2::$BT)}()
        Base.xor(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{(|)(V1::$BT, V2::$BT)}()
    end
end
function Base.fld(x::SInteger, y::SInteger)
    d = div(x, y)
    return d - (signbit(x ⊻ y) & (d * y != x))
end

function cld(x::SInteger, y::SInteger)
    d = div(x, y)
    return d + (((x > 0) == (y > 0)) & (d * y != x))
end

Base.fld(x::SSigned, y::SUnsigned) = div(x, y) - (signbit(x) & (rem(x, y) != 0))
Base.fld(x::SUnsigned, y::SSigned) = div(x, y) - (signbit(y) & (rem(x, y) != 0))

for ST in static_unsigned
    @eval begin
        Base.fld(x::$ST, y::$ST) = div(x,y)
        function cld(x::$ST, y::$ST) where T
            d = div(x, y)
            return d + (d * y != x)
        end
    end
end

Base.mod(x::BitSSigned, y::SUnsigned) = rem(y + unsigned(rem(x, y)), y)
Base.mod(x::SUnsigned, y::SSigned) = rem(y + signed(rem(x, y)), y)
Base.mod(x::SUnsigned, y::SUnsigned) = rem(x, y)

Base.cld(x::SSigned, y::SUnsigned) = div(x, y) + (!signbit(x) & (rem(x, y) != 0))
Base.cld(x::SUnsigned, y::SSigned) = div(x, y) + (!signbit(y) & (rem(x, y) != 0))


text:
"""
	normcdf(x)

Probablity that Z ≤ x, where Z is a standard normal random variable.
"""
normcdf(x::Real) = erfc(-x / √two(x)) / 2

"""
	normcdf(a, b)

Probablity that a ≤ Z ≤ b, where Z is a standard normal random variable.
WARNING: Silently returns a negative value if a > b.
"""
normcdf(a::Real, b::Real) = erf(a / √two(a), b / √two(b)) / 2

"""
	normcdfinv(x)

Inverse of normcdf.
"""
normcdfinv(x::Real) = -√two(x) * erfcinv(2x)

"""
	lognormcdf(a, b)

Computes log(normcdf(a, b)), but retaining accuracy.
"""
lognormcdf(a::Real, b::Real) = logerf(a / √two(a), b / √two(b)) - log(two(a + b))

"""
	tnmean(a)

Mean of the standard normal distribution,
truncated to the interval (a, +∞).
"""
tnmean(a::Real) = sqrt(two(a)/π) / erfcx(a/√two(a))

"""
	tnvar(a)

Variance of the standard normal distribution,
truncated to the interval (a, +∞).
WARNING: Fails for very very large values of `a`.
"""
function tnvar(a::Real)
	μ = tnmean(a)
	return one(μ) - (μ - a) * μ
end

"""
	tnstd(a)

Standard deviation of the standard normal distribution,
truncated to the interval (a, +∞).
"""
tnstd(a::Real) = √tnvar(a)

"""
	mills(x::Real)

Mills ratio of the standard normal distribution.
Defined as (1 - cdf(x)) / pdf(x).
"""
mills(x::Real) = erfcx(x / √two(x)) * √(π/two(x))


text:
__precompile__()

module FinancialMonteCarlo

using Requires # for conditional dependencies
using Random # for randn! and related
using LinearAlgebra # for Longstaff Schwartz
function __init__()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("deps/cuda_dependencies.cujl")
    @require ArrayFire = "b19378d9-d87a-599a-927f-45f220a2c452" include("deps/af_dependencies.cujl")
    @require DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" include("deps/dual_dependencies.jl")
    @require DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" include("deps/diffeq_dependencies.jl")
    @require Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" include("deps/distributed_dependencies.jl")
end
include("utils.jl")
include("models.jl")
include("options.jl")
include("io_utils.jl")
include("operations.jl")
include("metrics.jl")
include("output_type.jl")
include("multi_threading.jl")
export pricer, variance, confinter, delta, simulate, payoff, →

end#End Module


text:
using PWS
using Catalyst
using DrWatson
using LinearAlgebra

import Distributions: MvNormal, Poisson

sn = @reaction_network begin
    κ, ∅ --> S
    λ, S --> ∅
end κ λ

rn = @reaction_network begin
    ρ, S --> X + S
    μ, X --> ∅
end ρ μ

function get_system(mean_s, mean_x, tau_s, tau_x, duration=tau_s)
    κ = mean_s / tau_s
    λ = 1 / tau_s
    ρ = (mean_x / mean_s) / tau_x
    μ = 1 / tau_x
    PWS.JumpSystem(sn, rn, [κ, λ], [ρ, μ], mean_s, mean_x, duration)
end

nothing


text: # programas isso daqui---