text:
struct XpressError <: Exception
    errorcode::Int
    msg::String
    # TODO try to append string from Xpress.getlasterror(model) here
    # only possible on call that had the model as parameter
end

function Base.showerror(io::IO, err::XpressError)
    print(io, "XpressError($(err.errorcode)): ")
    if err.errorcode == 1
        print(io, "Bad input encountered.")
    elseif err.errorcode == 2
        print(io, "Bad or corrupt file - unrecoverable.")
    elseif err.errorcode == 4
        print(io, "Memory error.")
    elseif err.errorcode == 8
        print(io, "Corrupt use.")
    elseif err.errorcode == 16
        print(io, "Program error.")
    elseif err.errorcode == 32
        print(io, "Subroutine not completed successfully, possibly due to invalid argument.")
    elseif err.errorcode == 128
        print(io, "Too many users.")
    else
        print(io, "Unrecoverable error.")
    end
    print(io, " $(err.msg)")
end

function fixinfinity(val::Float64)
    if val == Inf
        return Xpress.Lib.XPRS_PLUSINFINITY
    elseif val == -Inf
        return Xpress.Lib.XPRS_MINUSINFINITY
    else
        return val
    end
end

function fixinfinity!(vals::Vector{Float64})
    map!(fixinfinity, vals, vals)
end

"""
    Xpress.CWrapper

abstract type Xpress.CWrapper
"""
abstract type CWrapper end

Base.unsafe_convert(::Type{Ptr{Cvoid}}, t::CWrapper) = (t.ptr == C_NULL) ? throw(XpressError(255, "Received null pointer in CWrapper. Something must be wrong.")) : t.ptr

mutable struct XpressProblem <: CWrapper
    ptr::Lib.XPRSprob
    logfile::String
    function XpressProblem(ptr::Lib.XPRSprob; finalize_env::Bool = true, logfile = "")
        model = new(ptr, logfile)
        if finalize_env
            finalizer(destroyprob, model)
        end
        return model
    end
end

function get_xpress_error_message(prob::XpressProblem)
    lstrip(Xpress.getlasterror(prob), ['?'])
end

function XpressProblem(; logfile = "")
    ref = Ref{Lib.XPRSprob}()
    createprob(ref)
    @assert ref[] != C_NULL "Failed to create XpressProblem. Received null pointer from Xpress C interface."
    p = XpressProblem(ref[], logfile = logfile)
    if logfile != ""
        setlogfile(p, logfile)
    end
    return p
end

addcolnames(prob::XpressProblem, names::Vector{String}) = addnames(prob, names, 2)
addrownames(prob::XpressProblem, names::Vector{String}) = addnames(prob, names, 1)

"""
    getcontrol(prob::XpressProblem, control::Integer)
    getcontrol(prob::XpressProblem, control::String)
    getcontrol(prob::XpressProblem, control::Symbol)

Get parameter of any type
"""
function getcontrol(prob::XpressProblem, control::Integer)
    # TODO: this function is not type stable
    if convert(Int, control) in XPRS_INT_CONTROLS
        return getintcontrol(prob, control)
    elseif convert(Int, control) in XPRS_DBL_CONTROLS
        return getdblcontrol(prob, control)
    elseif convert(Int, control) in XPRS_STR_CONTROLS
        return getstrcontrol(prob, control)
    else
        error("Unrecognized control parameter: $(control).")
    end
end
getcontrol(prob::XpressProblem, control::Symbol) = getcontrol(prob, getproperty(Lib, Symbol("XPRS_$(String(control))")))
getcontrol(prob::XpressProblem, control::String) = getcontrol(prob, XPRS_ATTRIBUTES[control])

"""
    setcontrol!(prob::XpressProblem, control::Symbol, val::Any)
    setcontrol!(prob::XpressProblem, control::String, val::Any)
    setcontrol!(prob::XpressProblem, control::Integer, val::Any)

Set parameter of any type
"""
setcontrol!(prob::XpressProblem, control::Symbol, val::Any) = setcontrol!(prob, getproperty(Lib, Symbol("XPRS_$(String(control))")), val::Any)
setcontrol!(prob::XpressProblem, control::String, val::Any) = setcontrol!(prob, XPRS_ATTRIBUTES[control], val)
setcontrol!(prob::XpressProblem, control::Integer, val::Any) = setcontrol!(prob, Cint(control), val)
function setcontrol!(prob::XpressProblem, control::Cint, val)
    # TODO: dispatch on Val(control) instead?
    if convert(Int, control) in XPRS_INT_CONTROLS
        if isinteger(val)
            setintcontrol(prob, control, convert(Int, val))
        else
            error("Expected and integer and got $val")
        end
    elseif convert(Int, control) in XPRS_DBL_CONTROLS
        setdblcontrol(prob, control, val)
    elseif convert(Int, control) in XPRS_STR_CONTROLS
        setstrcontrol(prob, control, val)
    else
        error("Unrecognized control parameter: $(control).")
    end
end

"""
    setparams!(prob::XpressProblem;args...)

Set multiple parameters of any type
"""
function setcontrols!(prob::XpressProblem; args...)
    for (control, val) in args
        setcontrols!(prob, getproperty(Lib, control), val)
    end
end

XPRS_STR_CONTROLS = [
                        Lib.XPRS_MPSRHSNAME,
                        Lib.XPRS_MPSOBJNAME,
                        Lib.XPRS_MPSRANGENAME,
                        Lib.XPRS_MPSBOUNDNAME,
                        Lib.XPRS_OUTPUTMASK,
                    ]

XPRS_DBL_CONTROLS = [
                        Lib.XPRS_MATRIXTOL
                        Lib.XPRS_PIVOTTOL
                        Lib.XPRS_FEASTOL
                        Lib.XPRS_OUTPUTTOL
                        Lib.XPRS_SOSREFTOL
                        Lib.XPRS_OPTIMALITYTOL
                        Lib.XPRS_ETATOL
                        Lib.XPRS_RELPIVOTTOL
                        Lib.XPRS_MIPTOL
                        Lib.XPRS_MIPTOLTARGET
                        Lib.XPRS_MIPADDCUTOFF
                        Lib.XPRS_MIPABSCUTOFF
                        Lib.XPRS_MIPRELCUTOFF
                        Lib.XPRS_PSEUDOCOST
                        Lib.XPRS_PENALTY
                        Lib.XPRS_BIGM
                        Lib.XPRS_MIPABSSTOP
                        Lib.XPRS_MIPRELSTOP
                        Lib.XPRS_CROSSOVERACCURACYTOL
                        Lib.XPRS_BAROBJSCALE
                        Lib.XPRS_BARRHSSCALE
                        Lib.XPRS_CHOLESKYTOL
                        Lib.XPRS_BARGAPSTOP
                        Lib.XPRS_BARDUALSTOP
                        Lib.XPRS_BARPRIMALSTOP
                        Lib.XPRS_BARSTEPSTOP
                        Lib.XPRS_ELIMTOL
                        Lib.XPRS_PERTURB
                        Lib.XPRS_MARKOWITZTOL
                        Lib.XPRS_MIPABSGAPNOTIFY
                        Lib.XPRS_MIPRELGAPNOTIFY
                        Lib.XPRS_PPFACTOR
                        Lib.XPRS_REPAIRINDEFINITEQMAX
                        Lib.XPRS_BARGAPTARGET
                        Lib.XPRS_SBEFFORT
                        Lib.XPRS_HEURDIVERANDOMIZE
                        Lib.XPRS_HEURSEARCHEFFORT
                        Lib.XPRS_CUTFACTOR
                        Lib.XPRS_EIGENVALUETOL
                        Lib.XPRS_INDLINBIGM
                        Lib.XPRS_TREEMEMORYSAVINGTARGET
                        Lib.XPRS_GLOBALFILEBIAS
                        Lib.XPRS_INDPRELINBIGM
                        Lib.XPRS_RELAXTREEMEMORYLIMIT
                        Lib.XPRS_MIPABSGAPNOTIFYOBJ
                        Lib.XPRS_MIPABSGAPNOTIFYBOUND
                        Lib.XPRS_PRESOLVEMAXGROW
                        Lib.XPRS_HEURSEARCHTARGETSIZE
                        Lib.XPRS_CROSSOVERRELPIVOTTOL
                        Lib.XPRS_CROSSOVERRELPIVOTTOLSAFE
                        Lib.XPRS_DETLOGFREQ
                        Lib.XPRS_MAXIMPLIEDBOUND
                        Lib.XPRS_FEASTOLTARGET
                        Lib.XPRS_OPTIMALITYTOLTARGET
                        Lib.XPRS_PRECOMPONENTSEFFORT
                    ]

XPRS_INT_CONTROLS = [
                        Lib.XPRS_EXTRAROWS
                        Lib.XPRS_EXTRACOLS
                        Lib.XPRS_LPITERLIMIT
                        Lib.XPRS_LPLOG
                        Lib.XPRS_SCALING
                        Lib.XPRS_PRESOLVE
                        Lib.XPRS_CRASH
                        Lib.XPRS_PRICINGALG
                        Lib.XPRS_INVERTFREQ
                        Lib.XPRS_INVERTMIN
                        Lib.XPRS_MAXNODE
                        Lib.XPRS_MAXTIME
                        Lib.XPRS_MAXMIPSOL
                        Lib.XPRS_DEFAULTALG
                        Lib.XPRS_VARSELECTION
                        Lib.XPRS_NODESELECTION
                        Lib.XPRS_BACKTRACK
                        Lib.XPRS_MIPLOG
                        Lib.XPRS_KEEPNROWS
                        Lib.XPRS_MPSECHO
                        Lib.XPRS_MAXPAGELINES
                        Lib.XPRS_OUTPUTLOG
                        Lib.XPRS_BARSOLUTION
                        Lib.XPRS_CACHESIZE
                        Lib.XPRS_CROSSOVER
                        Lib.XPRS_BARITERLIMIT
                        Lib.XPRS_CHOLESKYALG
                        Lib.XPRS_BAROUTPUT
                        Lib.XPRS_CSTYLE
                        Lib.XPRS_EXTRAMIPENTS
                        Lib.XPRS_REFACTOR
                        Lib.XPRS_BARTHREADS
                        Lib.XPRS_KEEPBASIS
                        Lib.XPRS_VERSION
                        Lib.XPRS_BIGMMETHOD
                        Lib.XPRS_MPSNAMELENGTH
                        Lib.XPRS_PRESOLVEOPS
                        Lib.XPRS_MIPPRESOLVE
                        Lib.XPRS_MIPTHREADS
                        Lib.XPRS_BARORDER
                        Lib.XPRS_BREADTHFIRST
                        Lib.XPRS_AUTOPERTURB
                        Lib.XPRS_DENSECOLLIMIT
                        Lib.XPRS_CALLBACKFROMMASTERTHREAD
                        Lib.XPRS_MAXMCOEFFBUFFERELEMS
                        Lib.XPRS_REFINEOPS
                        Lib.XPRS_LPREFINEITERLIMIT
                        Lib.XPRS_MIPREFINEITERLIMIT
                        Lib.XPRS_DUALIZEOPS
                        Lib.XPRS_PRESORT
                        Lib.XPRS_PREPERMUTE
                        Lib.XPRS_PREPERMUTESEED
                        Lib.XPRS_MAXMEMORY
                        Lib.XPRS_MAXMEMORYSOFT
                        Lib.XPRS_MAXMEMORYHARD
                        Lib.XPRS_CUTFREQ
                        Lib.XPRS_SYMSELECT
                        Lib.XPRS_SYMMETRY
                        Lib.XPRS_LPTHREADS
                        Lib.XPRS_MIQCPALG
                        Lib.XPRS_QCCUTS
                        Lib.XPRS_QCROOTALG
                        Lib.XPRS_ALGAFTERNETWORK
                        Lib.XPRS_TRACE
                        Lib.XPRS_MAXIIS
                        Lib.XPRS_CPUTIME
                        Lib.XPRS_COVERCUTS
                        Lib.XPRS_GOMCUTS
                        Lib.XPRS_MPSFORMAT
                        Lib.XPRS_CUTSTRATEGY
                        Lib.XPRS_CUTDEPTH
                        Lib.XPRS_TREECOVERCUTS
                        Lib.XPRS_TREEGOMCUTS
                        Lib.XPRS_CUTSELECT
                        Lib.XPRS_TREECUTSELECT
                        Lib.XPRS_DUALIZE
                        Lib.XPRS_DUALGRADIENT
                        Lib.XPRS_SBITERLIMIT
                        Lib.XPRS_SBBEST
                        Lib.XPRS_MAXCUTTIME
                        Lib.XPRS_ACTIVESET
                        Lib.XPRS_BARINDEFLIMIT
                        Lib.XPRS_HEURSTRATEGY
                        Lib.XPRS_HEURFREQ
                        Lib.XPRS_HEURDEPTH
                        Lib.XPRS_HEURMAXSOL
                        Lib.XPRS_HEURNODES
                        Lib.XPRS_LNPBEST
                        Lib.XPRS_LNPITERLIMIT
                        Lib.XPRS_BRANCHCHOICE
                        Lib.XPRS_BARREGULARIZE
                        Lib.XPRS_SBSELECT
                        Lib.XPRS_LOCALCHOICE
                        Lib.XPRS_LOCALBACKTRACK
                        Lib.XPRS_DUALSTRATEGY
                        Lib.XPRS_L1CACHE
                        Lib.XPRS_HEURDIVESTRATEGY
                        Lib.XPRS_HEURSELECT
                        Lib.XPRS_BARSTART
                        Lib.XPRS_BARNUMSTABILITY
                        Lib.XPRS_BARORDERTHREADS
                        Lib.XPRS_EXTRASETS
                        Lib.XPRS_FEASIBILITYPUMP
                        Lib.XPRS_PRECOEFELIM
                        Lib.XPRS_PREDOMCOL
                        Lib.XPRS_HEURSEARCHFREQ
                        Lib.XPRS_HEURDIVESPEEDUP
                        Lib.XPRS_SBESTIMATE
                        Lib.XPRS_BARCORES
                        Lib.XPRS_MAXCHECKSONMAXTIME
                        Lib.XPRS_MAXCHECKSONMAXCUTTIME
                        Lib.XPRS_HISTORYCOSTS
                        Lib.XPRS_ALGAFTERCROSSOVER
                        Lib.XPRS_LINELENGTH
                        Lib.XPRS_MUTEXCALLBACKS
                        Lib.XPRS_BARCRASH
                        Lib.XPRS_HEURDIVESOFTROUNDING
                        Lib.XPRS_HEURSEARCHROOTSELECT
                        Lib.XPRS_HEURSEARCHTREESELECT
                        Lib.XPRS_MPS18COMPATIBLE
                        Lib.XPRS_ROOTPRESOLVE
                        Lib.XPRS_CROSSOVERDRP
                        Lib.XPRS_FORCEOUTPUT
                        Lib.XPRS_DETERMINISTIC
                        Lib.XPRS_PREPROBING
                        Lib.XPRS_EXTRAQCELEMENTS
                        Lib.XPRS_EXTRAQCROWS
                        Lib.XPRS_TREEMEMORYLIMIT
                        Lib.XPRS_TREECOMPRESSION
                        Lib.XPRS_TREEDIAGNOSTICS
                        Lib.XPRS_MAXGLOBALFILESIZE
                        Lib.XPRS_REPAIRINFEASMAXTIME
                        Lib.XPRS_IFCHECKCONVEXITY
                        Lib.XPRS_PRIMALUNSHIFT
                        Lib.XPRS_REPAIRINDEFINITEQ
                        Lib.XPRS_MAXLOCALBACKTRACK
                        Lib.XPRS_USERSOLHEURISTIC
                        Lib.XPRS_FORCEPARALLELDUAL
                        Lib.XPRS_BACKTRACKTIE
                        Lib.XPRS_BRANCHDISJ
                        Lib.XPRS_MIPFRACREDUCE
                        Lib.XPRS_CONCURRENTTHREADS
                        Lib.XPRS_MAXSCALEFACTOR
                        Lib.XPRS_HEURTHREADS
                        Lib.XPRS_THREADS
                        Lib.XPRS_HEURBEFORELP
                        Lib.XPRS_PREDOMROW
                        Lib.XPRS_BRANCHSTRUCTURAL
                        Lib.XPRS_QUADRATICUNSHIFT
                        Lib.XPRS_BARPRESOLVEOPS
                        Lib.XPRS_QSIMPLEXOPS
                        Lib.XPRS_CONFLICTCUTS
                        Lib.XPRS_PREPROTECTDUAL
                        Lib.XPRS_CORESPERCPU
                        Lib.XPRS_RESOURCESTRATEGY
                        Lib.XPRS_SLEEPONTHREADWAIT
                        Lib.XPRS_PREDUPROW
                        Lib.XPRS_CPUPLATFORM
                        Lib.XPRS_BARALG
                        Lib.XPRS_TREEPRESOLVE
                        Lib.XPRS_TREEPRESOLVE_KEEPBASIS
                        Lib.XPRS_TREEPRESOLVEOPS
                        Lib.XPRS_LPLOGSTYLE
                        Lib.XPRS_RANDOMSEED
                        Lib.XPRS_TREEQCCUTS
                        Lib.XPRS_PRELINDEP
                        Lib.XPRS_DUALTHREADS
                        Lib.XPRS_PREOBJCUTDETECT
                        Lib.XPRS_PREBNDREDQUAD
                        Lib.XPRS_PREBNDREDCONE
                        Lib.XPRS_PRECOMPONENTS
                        Lib.XPRS_MAXMIPTASKS
                        Lib.XPRS_MIPTERMINATIONMETHOD
                        Lib.XPRS_PRECONEDECOMP
                        Lib.XPRS_HEURSEARCHROOTCUTFREQ
                        Lib.XPRS_EXTRAELEMS
                        Lib.XPRS_EXTRAPRESOLVE
                        Lib.XPRS_EXTRASETELEMS
                    ]

# TODO list attributes by type

# originals are more important to be used everywhere, presolved are actually
# secondary
n_variables(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALCOLS)
n_constraints(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALROWS)
n_special_ordered_sets(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALSETS)
n_quadratic_constraints(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALQCONSTRAINTS)
n_non_zero_elements(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ELEMS)
n_quadratic_elements(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALQELEMS)
n_quadratic_row_coefficients(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALQCELEMS)
n_entities(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALMIPENTS)
n_setmembers(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALSETMEMBERS)

n_original_variables(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALCOLS)
n_original_constraints(prob::XpressProblem) = getintattrib(prob, Lib.XPRS_ORIGINALROWS)

objective_sense(prob::XpressProblem) = getdblattrib(prob, Lib.XPRS_OBJSENSE) == Lib.XPRS_OBJ_MINIMIZE ? :minimize : :maximize

# derived attribute functions

"""
    n_linear_constraints(prob::XpressProblem)
Return the number of purely linear contraints in the XpressProblem
"""
n_linear_constraints(prob::XpressProblem) = n_constraints(prob) - n_quadratic_constraints(prob)

"""
    is_qcp(prob::XpressProblem)
Return `true` if there are quadratic constraints in the XpressProblem
"""
is_quadratic_constraints(prob::XpressProblem) = n_quadratic_constraints(prob) > 0

"""
    is_mip(prob::XpressProblem)
Return `true` if there are integer entities in the XpressProblem
"""
is_mixedinteger(prob::XpressProblem) = (n_entities(prob) + n_special_ordered_sets(prob)) > 0

"""
    is_quadratic_objective(prob::XpressProblem)
Return `true` if there are quadratic terms in the objective in the XpressProblem
"""
is_quadratic_objective(prob::XpressProblem) = n_quadratic_elements(prob) > 0

"""
    problem_type(prob::XpressProblem)
Return a symbol enconding the type of the problem.]
Options are: `:LP`, `:QP` and `:QCP`
"""
function problem_type(prob::XpressProblem)
    is_quadratic_constraints(prob) ? (:QCP) :
    is_quadratic_objective(prob)  ? (:QP)  : (:LP)
end

"""
    show(io::IO, prob::XpressProblem)
Prints a simplified problem description
"""
function Base.show(io::IO, prob::XpressProblem)
    println(io, "Xpress Problem:"     )
    if is_mixedinteger(prob)
    println(io, "    type   : $(problem_type(prob)) (MIP)")
    else
    println(io, "    type   : $(problem_type(prob))")
    end
    println(io, "    sense  : $(objective_sense(prob))")
    println(io, "    number of variables                    = $(n_variables(prob))")
    println(io, "    number of linear constraints           = $(n_linear_constraints(prob))")
    println(io, "    number of quadratic constraints        = $(n_quadratic_constraints(prob))")
    println(io, "    number of sos constraints              = $(n_special_ordered_sets(prob))")
    println(io, "    number of non-zero coeffs              = $(n_non_zero_elements(prob))")
    println(io, "    number of non-zero qp objective terms  = $(n_quadratic_elements(prob))")
    println(io, "    number of non-zero qp constraint terms = $(n_quadratic_row_coefficients(prob))")
    println(io, "    number of integer entities             = $(n_entities(prob))")
end

const MIPSTATUS_STRING = Dict{Int,String}(
    Xpress.Lib.XPRS_MIP_NOT_LOADED => "0 Problem has not been loaded ( XPRS_MIP_NOT_LOADED).",
    Xpress.Lib.XPRS_MIP_LP_NOT_OPTIMAL => "1 Global search incomplete - the initial continuous relaxation has not been solved and no integer solution has been found ( XPRS_MIP_LP_NOT_OPTIMAL).",
    Xpress.Lib.XPRS_MIP_LP_OPTIMAL => "2 Global search incomplete - the initial continuous relaxation has been solved and no integer solution has been found ( XPRS_MIP_LP_OPTIMAL).",
    Xpress.Lib.XPRS_MIP_NO_SOL_FOUND => "3 Global search incomplete - no integer solution found ( XPRS_MIP_NO_SOL_FOUND).",
    Xpress.Lib.XPRS_MIP_SOLUTION => "4 Global search incomplete - an integer solution has been found ( XPRS_MIP_SOLUTION).",
    Xpress.Lib.XPRS_MIP_INFEAS => "5 Global search complete - no integer solution found ( XPRS_MIP_INFEAS).",
    Xpress.Lib.XPRS_MIP_OPTIMAL => "6 Global search complete - integer solution found ( XPRS_MIP_OPTIMAL).",
    Xpress.Lib.XPRS_MIP_UNBOUNDED => "7 Global search incomplete - the initial continuous relaxation was found to be unbounded. A solution may have been found ( XPRS_MIP_UNBOUNDED).",
)

function mip_solve_complete(stat)
    stat in [Xpress.Lib.XPRS_MIP_INFEAS, Xpress.Lib.XPRS_MIP_OPTIMAL]
end
function mip_solve_stopped(stat)
    stat in [Xpress.Lib.XPRS_MIP_INFEAS, Xpress.Lib.XPRS_MIP_OPTIMAL]
end

const LPSTATUS_STRING = Dict{Int,String}(
    Xpress.Lib.XPRS_LP_UNSTARTED => "0 Unstarted ( XPRS_LP_UNSTARTED).",
    Xpress.Lib.XPRS_LP_OPTIMAL => "1 Optimal ( XPRS_LP_OPTIMAL).",
    Xpress.Lib.XPRS_LP_INFEAS => "2 Infeasible ( XPRS_LP_INFEAS).",
    Xpress.Lib.XPRS_LP_CUTOFF => "3 Objective worse than cutoff ( XPRS_LP_CUTOFF).",
    Xpress.Lib.XPRS_LP_UNFINISHED => "4 Unfinished ( XPRS_LP_UNFINISHED).",
    Xpress.Lib.XPRS_LP_UNBOUNDED => "5 Unbounded ( XPRS_LP_UNBOUNDED).",
    Xpress.Lib.XPRS_LP_CUTOFF_IN_DUAL => "6 Cutoff in dual ( XPRS_LP_CUTOFF_IN_DUAL).",
    Xpress.Lib.XPRS_LP_UNSOLVED => "7 Problem could not be solved due to numerical issues. ( XPRS_LP_UNSOLVED).",
    Xpress.Lib.XPRS_LP_NONCONVEX => "8 Problem contains quadratic data, which is not convex ( XPRS_LP_NONCONVEX).",
)

const STOPSTATUS_STRING = Dict{Int, String}(
    Xpress.Lib.XPRS_STOP_NONE => "no interruption - the solve completed normally",
    Xpress.Lib.XPRS_STOP_TIMELIMIT => "time limit hit",
    Xpress.Lib.XPRS_STOP_CTRLC => "control C hit",
    Xpress.Lib.XPRS_STOP_NODELIMIT => "node limit hit",
    Xpress.Lib.XPRS_STOP_ITERLIMIT => "iteration limit hit",
    Xpress.Lib.XPRS_STOP_MIPGAP => "MIP gap is sufficiently small",
    Xpress.Lib.XPRS_STOP_SOLLIMIT => "solution limit hit",
    Xpress.Lib.XPRS_STOP_USER => "user interrupt.",
)


text:
const FT = Float64
import SciMLBase: step!
using OrdinaryDiffEq: ODEProblem, solve, SSPRK33, savevalues!
import ClimaCore.Utilities: PlusHalf
using DiffEqCallbacks

push!(LOAD_PATH, joinpath(@__DIR__, "..", "..", ".."))
using ClimaCoupler

include("atmos_rhs.jl")
include("ocean_rhs.jl")
include("land_rhs.jl")
include("coupledbc.jl")

# model parameters that the coupler overwrites
const atm_T_ini = FT(270.0)
const MSLP = FT(1e5)
const grav = FT(9.8)
const R_d = FT(287.058)
const γ = FT(1.4)
const C_p = FT(R_d * γ / (γ - 1))
const C_v = FT(R_d / (γ - 1))
const R_m = R_d
cpl_parameters = (
    # atmos parameters
    atm_μ = FT(0.0001), # diffusion coefficient
    atm_T_top = FT(280.0), # fixed temperature at the top of the domain_atm
    atm_T_ini = atm_T_ini, # initial condition of at temperature (isothermal) [K]
    MSLP = MSLP, # mean sea level pressure
    grav = grav, # gravitational constant
    R_d = R_d, # R dry (gas constant / mol mass dry air)
    γ = γ, # heat capacity ratio
    C_p = C_p, # heat capacity at constant pressure
    C_v = C_v, # heat capacity at constant volume
    R_m = R_m, # moist R, assumed to be dry
    # land slab parameters
    lnd_h = FT(0.5), # depth of slab layer [m]
    lnd_ρ = FT(1500), # density [kg m^-3]
    lnd_c = FT(800), # specific heat [J K^-1 kg^-1]
    lnd_T_ini = FT(260.0), # initial condition of at temperature (isothermal) [K]
    # ocean slab parameters
    ocn_h = FT(0.5), # depth of slab layer [m]
    ocn_ρ = FT(1025), # density [kg m^-3]
    ocn_c = FT(3850), # specific heat [J K^-1 kg^-1]
    ocn_T_ini = FT(260.0), # initial condition of at temperature (isothermal) [K]
    # coupling parameters
    C_H = FT(0.0015),
)

function make_dss_func()
    function _dss!(x::Fields.Field)
        Spaces.weighted_dss!(x)
    end
    function _dss!(::Any)
        nothing
    end
    dss_func(Y, t, integrator) = foreach(_dss!, Fields._values(Y))
    return dss_func
end
dss_func = make_dss_func()
dss_callback = FunctionCallingCallback(dss_func, func_start = true)

@info "Init Models and Maps"

# timestepping
t_start, t_end = (0.0, 1.0)
cpl_Δt = 0.1
saveat = 1e2
atm_nsteps, ocn_nsteps, lnd_nsteps = (5, 1, 1)

atm_Y_default, atm_bc, atm_domain = atm_init(
    xmin = -500,
    xmax = 500,
    zmin = 0,
    zmax = 1000,
    npoly = 4,
    helem = 20,
    velem = 20,
    bc = (ρθ = (bottom = CoupledFlux(), top = ZeroFlux()),),
)

ocn_Y_default, ocn_domain = ocn_init(xmin = -500, xmax = 0, helem = 10, npoly = 0)

lnd_Y_default, lnd_domain = lnd_init(xmin = 0, xmax = 500, helem = 10, npoly = 0)

# Build remapping operators
atm_boundary = Spaces.level(atm_domain.hv_face_space, PlusHalf(0))

atm_to_ocn = Operators.LinearRemap(ocn_domain, atm_boundary)
atm_to_lnd = Operators.LinearRemap(lnd_domain, atm_boundary)
ocn_to_atm = Operators.LinearRemap(atm_boundary, ocn_domain)
lnd_to_atm = Operators.LinearRemap(atm_boundary, lnd_domain)

# initialize coupling fields
atm_T_sfc = Operators.remap(ocn_to_atm, ocn_Y_default.T_sfc) .+ Operators.remap(lnd_to_atm, lnd_Y_default.T_sfc) # masked arrays; regrid to atm grid
atm_F_sfc = Fields.zeros(atm_boundary)
ocn_F_sfc = Fields.zeros(ocn_domain)
lnd_F_sfc = Fields.zeros(lnd_domain)

# init models
atm_Y = Fields.FieldVector(Yc = atm_Y_default.Yc, ρw = atm_Y_default.ρw, F_sfc = atm_F_sfc)
atm_p = (cpl_p = cpl_parameters, T_sfc = atm_T_sfc, bc = atm_bc)
atmos = AtmosSimulation(atm_Y, t_start, cpl_Δt / atm_nsteps, t_end, SSPRK33(), atm_p, saveat, dss_callback)

ocn_Y = Fields.FieldVector(T_sfc = ocn_Y_default.T_sfc)
ocn_p = (cpl_parameters, F_sfc = ocn_F_sfc)
ocean = OceanSimulation(ocn_Y, t_start, cpl_Δt / ocn_nsteps, t_end, SSPRK33(), ocn_p, saveat)

lnd_Y = Fields.FieldVector(T_sfc = lnd_Y_default.T_sfc)
lnd_p = (cpl_parameters, F_sfc = lnd_F_sfc)
land = LandSimulation(lnd_Y, t_start, cpl_Δt / lnd_nsteps, t_end, SSPRK33(), lnd_p, saveat)

# coupled simulation
struct AOLCoupledSimulation{FT, A <: AtmosSimulation, O <: OceanSimulation, L <: LandSimulation} <:
       ClimaCoupler.AbstractCoupledSimulation
    # Atmosphere Simulation
    atmos::A
    # Ocean Simulation
    ocean::O
    # Land Simulation
    land::L
    # The coupled time step size
    Δt::FT
end
sim = AOLCoupledSimulation(atmos, ocean, land, cpl_Δt)

# step for sims built on OrdinaryDiffEq
function step!(sim::ClimaCoupler.AbstractSimulation, t_stop)
    Δt = t_stop - sim.integrator.t
    step!(sim.integrator, Δt, true)
end

function cpl_run(simulation::AOLCoupledSimulation)
    @info "Run model"
    @unpack atmos, ocean, land, Δt = simulation
    cpl_Δt = Δt
    # coupler stepping
    for t in ((t_start + cpl_Δt):cpl_Δt:t_end)

        ## Atmos
        # pre: reset flux accumulator
        atmos.integrator.u.F_sfc .= 0.0 # reset surface flux to be accumulated

        # run
        # NOTE: use (t - integ_atm.t) here instead of Δt_cpl to avoid accumulating roundoff error in our timestepping.
        step!(atmos, t)

        ## Ocean
        # pre: get accumulated flux from atmos
        ocn_F_sfc = ocean.integrator.p.F_sfc
        Operators.remap!(ocn_F_sfc, atm_to_ocn, atmos.integrator.u.F_sfc ./ cpl_Δt)

        # run
        step!(ocean, t)
        # post: send ocean surface temp to atmos
        Operators.remap!(atmos.integrator.p.T_sfc, ocn_to_atm, ocean.integrator.u.T_sfc)

        ## Land
        # pre: get accumulated flux from atmos
        lnd_F_sfc = land.integrator.p.F_sfc
        Operators.remap!(lnd_F_sfc, atm_to_lnd, atmos.integrator.u.F_sfc ./ cpl_Δt)

        # run
        step!(land, t)
        # post: send ocean surface temp to atmos
        atmos.integrator.p.T_sfc .+= Operators.remap(lnd_to_atm, land.integrator.u.T_sfc)
    end
    @info "Simulation Complete"
end

cpl_run(sim)

# sol = sim.atmos.integrator.sol

# dirname = "sea_breeze_2d"
# path = joinpath(@__DIR__, "output", dirname)
# mkpath(path)

# using JLD2
# save(joinpath(path, "last_sim.jld2"), "coupled_sim", sim)

# post-processing
# import Plots, ClimaCorePlots
# Plots.GRBackend()
# interp = 5

# anim = Plots.@animate for u in sol.u
#     Plots.contourf(u.Yc.ρθ ./ u.Yc.ρ)
# end
# Plots.mp4(anim, joinpath(path, "theta.mp4"), fps = 20)

#  If2c = Operators.InterpolateF2C()
#  anim = Plots.@animate for u in sol.u
#      Plots.contourf(If2c.(u.ρw) ./ u.Yc.ρ)
#  end
#  Plots.mp4(anim, joinpath(path, "vel_w.mp4"), fps = 20)

#  anim = Plots.@animate for u in sol.u
#      Plots.contourf(u.Yc.ρuₕ ./ u.Yc.ρ)
#  end
#  Plots.mp4(anim, joinpath(path, "vel_u.mp4"), fps = 20)

ext: jl

max_stars_repo_name: CliMA/ClimaCoupler

max_issues_repo_name: CliMA/ClimaCoupler

max_forks_repo_name: CliMA/ClimaCoupler


text:
# error handling
# --------------

function find_first_topelevel_scope(bt::Vector{<:Union{Base.InterpreterIP,Ptr{Cvoid}}})
    for (i, ip) in enumerate(bt)
        st = Base.StackTraces.lookup(ip)
        ind = findfirst(st) do frame
            linfo = frame.linfo
            if linfo isa Core.CodeInfo
                linetable = linfo.linetable
                if isa(linetable, Vector) && length(linetable) ≥ 1
                    lin = first(linetable)
                    if isa(lin, Core.LineInfoNode) && lin.method === Symbol("top-level scope")
                        return true
                    end
                end
            else
                return frame.func === Symbol("top-level scope")
            end
        end
        ind === nothing || return i
    end
    return
end

# path utilitiles
# ---------------

function fullpath(path)
    return if isuntitled(path) || isabspath(path)
        path
    else
        basepath(path)
    end |> realpath′
end

isuntitled(path) = occursin(r"Untitled-\d+$", path)
basepath(path) = normpath(joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia", "base", path))

function realpath′(p)
    try
        ispath(p) ? realpath(p) : p
    catch e
        p
    end |> normpath
end

# string utilitiles
# -----------------

# https://github.com/JuliaDebug/Debugger.jl/blob/4cf99c662ab89da0fe7380c1e81461e2428e8b00/src/limitio.jl

mutable struct LimitIO{IO_t <: IO} <: IO
    io::IO_t
    maxbytes::Int
    n::Int
end
LimitIO(io::IO, maxbytes) = LimitIO(io, maxbytes, 0)

struct LimitIOException <: Exception end

function Base.write(io::LimitIO, v::UInt8)
    io.n > io.maxbytes && throw(LimitIOException())
    io.n += write(io.io, v)
end

function sprintlimited(args...; func=show, limit::Int=30, ellipsis::AbstractString="…", color=false)
    io = IOBuffer()
    ioctx = IOContext(LimitIO(io, limit - length(ellipsis)), :limit => true, :color => color, :displaysize => (30, 64))

    try
        Base.invokelatest(func, ioctx, args...)
    catch err
        if err isa LimitIOException
            print(io, ellipsis)
        else
            rethrow(err)
        end
    end

    str = filter(isvalid, String(take!(io)))

    return color ? str : remove_ansi_control_chars(str)
end

function strlimit(str; limit::Int=30, ellipsis::AbstractString="…")
    will_append = length(str) > limit

    io = IOBuffer()
    i = 1
    for c in str
        will_append && i > limit - length(ellipsis) && break
        isvalid(c) || continue

        print(io, c)
        i += 1
    end
    will_append && print(io, ellipsis)

    return String(take!(io))
end

# https://stackoverflow.com/a/33925425/12113178

function remove_ansi_control_chars(str::String)
    replace(str, r"(\x9B|\x1B\[)[0-?]*[ -\/]*[@-~]" => "")
end

function ends_with_semicolon(x)
    return REPL.ends_with_semicolon(split(x, '\n', keepempty=false)[end])
end

const splitlines = Base.Fix2(split, '\n')
const joinlines = Base.Fix2(join, '\n')

# VSCode specific
# ---------------

using Printf

const UNESCAPED = Set(codeunits("0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ-_.!~*'()"))

function encode_uri_component(uri)
    isvalid(uri) || throw(ArgumentError("`encode_uri_component` can only handle valid UTF8 strings."))

    io = IOBuffer()
    for cp in codeunits(uri)
        if cp in UNESCAPED
            print(io, Char(cp))
        else
            print(io, '%')
            @printf(io, "%2X", cp)
        end
    end
    return String(take!(io))
end

vscode_cmd_uri(cmd; cmdargs...) = string("command:", cmd, '?', encode_uri_component(JSON.json(cmdargs)))

# Misc handlers
function cd_to_uri(conn, params::NamedTuple{(:uri,),Tuple{String}})
    cd(params.uri)
    return nothing
end

function activate_uri(conn, params::NamedTuple{(:uri,),Tuple{String}})
    hideprompt(() -> Pkg.activate(params.uri))
    return nothing
end

max_stars_repo_name: FAlobaid/julia-vscode

max_issues_repo_name: FAlobaid/julia-vscode

max_forks_repo_name: FAlobaid/julia-vscode


text:
abstract type AbstractMessageType end

struct NotificationType{TPARAM} <: AbstractMessageType
    method::String
end

struct RequestType{TPARAM,TR} <: AbstractMessageType
    method::String
end

function NotificationType(method::AbstractString, ::Type{TPARAM}) where TPARAM
    return NotificationType{TPARAM}(method)
end

function RequestType(method::AbstractString, ::Type{TPARAM}, ::Type{TR}) where {TPARAM,TR}
    return RequestType{TPARAM,TR}(method)
end

get_param_type(::NotificationType{TPARAM}) where {TPARAM} = TPARAM
get_param_type(::RequestType{TPARAM,TR}) where {TPARAM,TR} = TPARAM
get_return_type(::RequestType{TPARAM,TR}) where {TPARAM,TR} = TR

function send(x::JSONRPCEndpoint, request::RequestType{TPARAM,TR}, params::TPARAM) where {TPARAM,TR}
    res = send_request(x, request.method, params)
    return typed_res(res, TR)::TR
end

# `send_request` must have returned nothing in this case, we pass this on
# so that we get an error in the typecast at the end of `send`
# if that is not the case.
typed_res(res, TR::Type{Nothing}) = res
typed_res(res, TR::Type{<:T}) where {T <: AbstractArray{Any}} = T(res)
typed_res(res, TR::Type{<:AbstractArray{T}}) where T = T.(res)
typed_res(res, TR::Type) = TR(res)

function send(x::JSONRPCEndpoint, notification::NotificationType{TPARAM}, params::TPARAM) where TPARAM
    send_notification(x, notification.method, params)
end

struct Handler
    message_type::AbstractMessageType
    func::Function
end

mutable struct MsgDispatcher
    _handlers::Dict{String,Handler}
    _currentlyHandlingMsg::Bool

    function MsgDispatcher()
        new(Dict{String,Handler}(), false)
    end
end

function Base.setindex!(dispatcher::MsgDispatcher, func::Function, message_type::AbstractMessageType)
    dispatcher._handlers[message_type.method] = Handler(message_type, func)
end

function dispatch_msg(x::JSONRPCEndpoint, dispatcher::MsgDispatcher, msg)
    dispatcher._currentlyHandlingMsg = true
    try
        method_name = msg["method"]
        handler = get(dispatcher._handlers, method_name, nothing)
        if handler !== nothing
            param_type = get_param_type(handler.message_type)
            params = param_type === Nothing ? nothing : param_type <: NamedTuple ? convert(param_type,(;(Symbol(i[1])=>i[2] for i in msg["params"])...)) : param_type(msg["params"])

            res = handler.func(x, params)

            if handler.message_type isa RequestType
                if res isa JSONRPCError
                    send_error_response(x, msg, res.code, res.msg, res.data)
                elseif res isa get_return_type(handler.message_type)
                    send_success_response(x, msg, res)
                else
                    error_msg = "The handler for the '$method_name' request returned a value of type $(typeof(res)), which is not a valid return type according to the request definition."
                    send_error_response(x, msg, -32603, error_msg, nothing)
                    error(error_msg)
                end
            end
        else
            error("Unknown method $method_name.")
        end
    finally
        dispatcher._currentlyHandlingMsg = false
    end
end

is_currently_handling_msg(d::MsgDispatcher) = d._currentlyHandlingMsg


text:
module InplaceQR

include("qr_fact.jl")
include("qr_api.jl")

end # module


text:
module BasisModule
export GaussianBasis, PrimitiveGaussianBasisFunction, ContractedGaussianBasisFunction, computeBasis, normalize!

using ..BaseModule
using ..BasisFunctionsModule
using ..BasisSetModule
using ..GeometryModule
using ..AtomModule
using ..IntegralsModule

abstract Basis

type GaussianBasis <: Basis
  contractedBFs::Array{ContractedGaussianBasisFunction,1}
end

function computeBasis(basSet::BasisSet,geo::Geometry)
  bas=GaussianBasis([])
  for atom in geo.atoms
    for contractedDefinition in basSet.definitions[atom.element]
      for mqn in MQuantumNumbers(contractedDefinition.lQuantumNumber)
	contractedBF = ContractedGaussianBasisFunction([],[])
	for primitiveDefinition in contractedDefinition.primitives
	  exponent=primitiveDefinition.exponent
	  primitiveBF = PrimitiveGaussianBasisFunction(atom.position,exponent,mqn)
	  norm = IntegralsModule.Overlap(primitiveBF,primitiveBF)
	  append!(contractedBF.coefficients,[primitiveDefinition.prefactor/sqrt(norm)])
	  append!(contractedBF.primitiveBFs,[primitiveBF])
	end
	append!(bas.contractedBFs,[contractedBF])
      end
    end
  end
  return bas
end

function normalize!(basis::GaussianBasis)
  for (cgb in basis.contractedBFs)
    normalize!(cgb)
  end
end

function normalize!(cgb::ContractedGaussianBasisFunction)
  N = Overlap(cgb,cgb)
  scale!(cgb.coefficients,1/sqrt(N))
end

end # module

max_stars_repo_name: tkelman/QuantumLab

max_issues_repo_name: tkelman/QuantumLab

max_forks_repo_name: tkelman/QuantumLab


text:
signbit(x::Mag) = false

signbit(x::ArbFloat{P}) where {P} = isfinite(x) ? sign_bit(x) : isneginf(x)

signbit(x::ArbReal{P}) where {P} = isfinite(x) ? sign_bit(x) : isneginf(x)

signbit(x::ArbComplex{P}, ::Type{RealPart}) where {P} = signbit(real(x))
signbit(x::ArbComplex{P}, ::Type{ImagPart}) where {P} = signbit(imag(x))
signbit(x::ArbComplex{P}) where {P} = signbit(x, RealPart)
signbits(x::ArbComplex{P}) where {P} = signbit(x, RealPart), signbit(x, ImagPart)

function sign(x::ArbFloat{P}) where {P}
    thesgn = ccall(@libarb(arf_sgn), Cint, (Ref{ArbFloat},), x)
    return ArbFloat{P}(thesgn)
end

#=
    void arb_sgn(arb_t y, const arb_t x)

    Sets y to the sign function of x.
    The result is [0±1]
       if x contains both zero and nonzero numbers.
=#
function sign(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_sgn), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

#=
    void acb_sgn(acb_t r, const acb_t z, slong prec)

    Sets r to the complex sign of z, defined as 0 if z is exactly zero
        and the projection onto the unit circle z/|z|=exp(iarg(z)) otherwise.

    void acb_csgn(arb_t r, const acb_t z)

    Sets r to the extension of the real sign function taking the value 1
        for z strictly in the right half plane, -1 for z strictly in the left half plane,
       and the sign of the imaginary part when z is on the imaginary axis.
=#
function sign(x::ArbComplex{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(acb_csgn), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x)
    return z
end

function signs(x::ArbComplex{P}) where {P}
    return sign(real(x)), sign(imag(x))
end

function (-)(x::ArbFloat{P}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_neg), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x)
    return z
end

function (-)(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_neg), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

function (-)(x::ArbComplex{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_neg), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}), z, x)
    return z
end

function abs(x::ArbFloat{P}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_abs), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x)
    return z
end

function abs(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_abs), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

function abs(x::ArbComplex{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_abs), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

abs2(x::ArbFloat{P})   where {P} = square( abs(x) )
abs2(x::ArbReal{P})    where {P} = square( abs(x) )
abs2(x::ArbComplex{P}) where {P} = square( abs(x) )

inv(x::ArbFloat{P}) where {P} = ArbFloat{P}( inv(ArbReal{P}(x)) )

function inv(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_inv), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, P)
    return z
end

function inv(x::ArbComplex{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_inv), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -23.82106505746807 <= q <= 20.87497636142748)

                     add_NL_constraint(m, :(sigmoid(-0.7432570628691102 + 0.9028765283094948*sigmoid(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.9177729976456397*sigmoid(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.6091404444372248*sigmoid(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + 0.7075296204914552*sigmoid(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + -0.9115635329002854*sigmoid(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + sigmoid(-0.2671327260424561 + -0.018499351907995898*sigmoid(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + 0.12529725503678657*sigmoid(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + 0.9839396030143157*sigmoid(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + 0.28500560463058955*sigmoid(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.970334521625404*sigmoid(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + sigmoid(0.40040910594674495 + -0.7454849743358074*sigmoid(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + 0.7108986625228169*sigmoid(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.30349624151173904*sigmoid(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.26057109207128715*sigmoid(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.7858005505148937*sigmoid(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + sigmoid(-0.7120753920031948 + -0.005433419997351141*sigmoid(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.8839514190404207*sigmoid(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.2933551592941672*sigmoid(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.5418193939730438*sigmoid(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.8978172795430157*sigmoid(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + sigmoid(0.7950909223530527 + -0.6380481687930124*sigmoid(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.8467501205760102*sigmoid(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + 0.5001198683847332*sigmoid(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.03855272462979853*sigmoid(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.3301524112226937*sigmoid(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using OptimizationProblems
using Optimize

# non convex problems
# Problems = [AMPGO02, AMPGO03, AMPGO04, AMPGO05, AMPGO06, AMPGO08, AMPGO09,
# AMPGO10, AMPGO11, AMPGO12, AMPGO13, AMPGO14, AMPGO15, AMPGO18,
# AMPGO20, AMPGO21, AMPGO22,
# beale, brownbs, broydn7d, brybnd,
# chainwoo, chnrosnb_mod, cosine, cragglvy, curly, dixmaane, dixmaanf, dixmaang,
# dixmaanh, dixmaani, dixmaanj, dixmaank, dixmaanl, dixmaanm, dixmaann, dixmaano, dixmaanp,
# Dus2_1, Dus2_3, Dus2_9, Duscube,
# edensch, eg2, errinros_mod, extrosnb, fletcbv2, fletcbv3_mod, fletchcr, fminsrf2, freuroth, genhumps, genrose, genrose_nash,
# indef_mod, liarwhd, meyer3, morebv, ncb20, ncb20b, noncvxu2, noncvxun, nondia, penalty2, powellsg, schmvett,
# scosine, Shpak1, Shpak2, Shpak3, Shpak4, Shpak5, Shpak6, sinquad, sparsine, sparsqur, srosenbr, tointgss, tquartic, woods]
## only if dimension > 1
Problems = [beale, brownbs, broydn7d, brybnd,
chainwoo, chnrosnb_mod, cosine, cragglvy, curly, dixmaane, dixmaanf, dixmaang,
dixmaanh, dixmaani, dixmaanj, dixmaank, dixmaanl, dixmaanm, dixmaann, dixmaano, dixmaanp,
edensch, eg2, errinros_mod, extrosnb, fletcbv2, fletcbv3_mod, fletchcr, fminsrf2, freuroth, genhumps, genrose, genrose_nash,
indef_mod, liarwhd, meyer3, morebv, ncb20, ncb20b, noncvxu2, noncvxun, nondia, penalty2, powellsg, schmvett,
scosine, sinquad, sparsine, sparsqur, srosenbr, tointgss, tquartic, woods]

# convex problems
# Problems = [arglina, arglinb, arglinc, arwhead, bdqrtic, brownden, cliff, clplatea, clplateb, clplatec, dixon3dq, dqdrtic, dqrtic, engval1, nasty, nondquar,
# palmer1c, palmer1d, palmer2c, palmer3c, palmer4c, palmer5c, palmer5d, palmer6c, palmer7c, palmer8c, power, quartc, tridia, vardim]

# unconstrained OptimizationProblems
# Problems = [AMPGO02, AMPGO03, AMPGO04, AMPGO05, AMPGO06, AMPGO08, AMPGO09,
# AMPGO10, AMPGO11, AMPGO12, AMPGO13, AMPGO14, AMPGO15, AMPGO18,
# AMPGO20, AMPGO21, AMPGO22,
# beale, brownbs, broydn7d, brybnd,
# chainwoo, chnrosnb_mod, cosine, cragglvy, curly, dixmaane, dixmaanf, dixmaang,
# dixmaanh, dixmaani, dixmaanj, dixmaank, dixmaanl, dixmaanm, dixmaann, dixmaano, dixmaanp,
# Dus2_1, Dus2_3, Dus2_9, Duscube,
# edensch, eg2, errinros_mod, extrosnb, fletcbv2, fletcbv3_mod, fletchcr, fminsrf2, freuroth, genhumps, genrose, genrose_nash,
# indef_mod, liarwhd, meyer3, morebv, ncb20, ncb20b, noncvxu2, noncvxun, nondia, penalty2, powellsg, schmvett,
# scosine, Shpak1, Shpak2, Shpak3, Shpak4, Shpak5, Shpak6, sinquad, sparsine, sparsqur, srosenbr, tointgss, tquartic, woods,
# arglina, arglinb, arglinc, arwhead, bdqrtic, brownden, cliff, clplatea, clplateb, clplatec, dixon3dq, dqdrtic, dqrtic, engval1, nasty, nondquar,
# palmer1c, palmer1d, palmer2c, palmer3c, palmer4c, palmer5c, palmer5d, palmer6c, palmer7c, palmer8c, power, quartc, tridia, vardim]

n = 100

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

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

end

stats = lin_two_solvers()
tf = font(13) # titlefont
f = font(12)
pb_type = "non-convex problems"
algo_used = "CG and CR"

p = performance_profile(hcat([p[:, 1] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("Performance profile : #f in the resolution of ", pb_type, " with ", algo_used),
                            titlefont = tf, legendfont = f, guidefont = f, size = (800,800)) #, legend=:bottomright, linealpha = 0.5, linewidth = 2)#, palette=:grays) # Profile for #f

Plots.savefig(p, string("profil_f_", pb_type, ".pdf"))

p = performance_profile(hcat([p[:, 2] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("Performance profile : #g in the resolution of ", pb_type, " with ", algo_used),
                            titlefont = tf, legendfont = f, guidefont = f, size = (800,800), legend=:bottomright) # Profile for #g

Plots.savefig(p, string("profil_g_", pb_type, ".pdf"))

p = performance_profile(hcat([p[:, 3] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("Performance profile : #Hv in the resolution of ", pb_type, " with ", algo_used),
                            titlefont = tf, legendfont = f, guidefont = f, size = (800,800), legend=:bottomright) # Profile for #Hv

Plots.savefig(p, string("profil_Hv_", pb_type, ".pdf"))

p = performance_profile(hcat([sum(p, 2) for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("Performance profile : #f + #g + #Hv in the resolution of ", pb_type, " with ", algo_used),
                            titlefont = tf, legendfont = f, guidefont = f, size = (800,800), legend=:bottomright) # Profile for #f + #g + #Hv

Plots.savefig(p, string("profil_f+g+Hv_", pb_type, ".pdf"))


text:
#TypeDefs
const UniformDict     = Dict{Symbol,Any}
const SymAnyDict      = Dict{Symbol,Any}
const EmptyNamedTuple = NamedTuple{(),Tuple{}}
const VecOrT{T}       = Union{Vector{T},T}

# Gapped Arrays are used in systems
abstract type RenderPassKind end
abstract type RenderTargetKind end

abstract type AbstractGlimpseMesh end

const INSTANCED_MESHES = Dict{Type,AbstractGlimpseMesh}()

struct BasicMesh{D,T,FT} <: AbstractGlimpseMesh
    vertices :: Vector{Point{D,T}}
    faces    :: Vector{FT}
    normals  :: Vector{Vec{D,T}}
end

struct AttributeMesh{AT<:NamedTuple,BM<:BasicMesh} <: AbstractGlimpseMesh
    attributes :: AT
    basic      :: BM
end

mutable struct Diorama <: AbstractLedger
    name::Symbol
    ledger::Ledger
    renderloop_stages::Vector{Stage}
    loop::Union{Task,Nothing}
    reupload::Bool
    function Diorama(name::Symbol, ledger::Ledger, renderloop_stages::Vector{Stage})
        dio = new(name, ledger, renderloop_stages, nothing, true)

        makecurrentdio(dio)
        finalizer(free!, dio)
        return dio
    end
end

include("components.jl")
include("singletons.jl")
include("meshes.jl")
include("system.jl")
include("diorama.jl")

glimpse_call(func::Function) = @tspawnat 2 func()


text:
module MySQLLoader

# https://github.com/JuliaDatabases/MySQL.jl
using MySQL # MySQL.jl 1.1
using .MySQL.DBInterface
using .MySQL.Tables

using Octo: Repo, AdapterBase, DBMS, SQLElement, Structured
using .Repo: SQLKeyword, ExecuteResult, sql_startswith_insert_update_delete
using .AdapterBase: INSERT

# db_dbname
function db_dbname(nt::NamedTuple)::String
    get(nt, :db, "")
end

# db_connect
function db_connect(; kwargs...)
    args = (:hostname, :username, :password)
    (hostname, username, password) = getindex.(Ref(kwargs), args)
    options = filter(kv -> !(kv.first in args), kwargs)
    DBInterface.connect(MySQL.Connection, hostname, username, password; options...)
end

# db_disconnect
function db_disconnect(conn)
    DBInterface.close!(conn)
end

# query
function query(conn, sql::String)
    (Tables.rowtable ∘ DBInterface.execute)(conn, sql)
end

function query(conn, prepared::String, vals::Vector)
    stmt = DBInterface.prepare(conn, prepared)
    (Tables.rowtable ∘ DBInterface.execute)(stmt, vals)
end

# execute

function execute(conn, sql::String)::ExecuteResult
    q = DBInterface.execute(conn, sql)
    sql_startswith_insert_update_delete_then_get_num_affected_rows(q.sql, conn)
end

function execute(conn, prepared::String, vals::Vector)::ExecuteResult
    stmt = DBInterface.prepare(conn, prepared)
    DBInterface.execute(stmt, vals)
    sql_startswith_insert_update_delete_then_get_num_affected_rows(prepared, conn)
end

function execute(conn, prepared::String, nts::Vector{<:NamedTuple})::ExecuteResult
    stmt = DBInterface.prepare(conn, prepared)
    num_affected_rows = 0
    for nt in nts
        DBInterface.execute(stmt, values(nt))
        num_affected_rows += get_num_affected_rows(conn)
    end
    (num_affected_rows=num_affected_rows,)
end

function sql_startswith_insert_update_delete_then_get_num_affected_rows(sql::String, conn)
    if sql_startswith_insert_update_delete(sql)
        num_affected_rows = get_num_affected_rows(conn)
        (num_affected_rows=num_affected_rows,)
    else
        nothing
    end
end

function get_num_affected_rows(conn)::Int
    MySQL.API.affectedrows(conn.mysql)
end

# execute_result
function execute_result(conn, command::SQLKeyword)::NamedTuple
    if INSERT === command
        last_insert_id = MySQL.API.insertid(conn.mysql)
        (id=last_insert_id,)
    else
        NamedTuple()
    end
end

function Base.show(io::IO, mime::MIME"text/plain", element::Union{E,Structured} where E<:SQLElement)
    dbms = DBMS.MySQL()
    AdapterBase._show(io, mime, dbms, element)
end

end # module Octo.Backends.MySQLLoader


text:
using BinaryOptimization
using Test
using LightGraphs
using Random
using LinearAlgebra

function simplified_tsp_qubo(w::Matrix{Float64})
    qubo = tsp_qubo(w)
    println(vars)
    n, _ = indexcapacities(vars(qubo))
    println(n)
    println(length(used_variables(qubo)))
    for i=1:(n-1)
        set_zero!(qubo, [1,i])
    end
    set_one!(qubo, [1,n])
    println(length(used_variables(qubo)))
    qubo
end

function simplified_tsp_hobo(w::Matrix{Float64})
    hobo = tsp_hobo(w)
    println(vars)
    n, logn = indexcapacities(vars(hobo))
    println(n)
    println(length(used_variables(hobo)))
    for i=1:logn
        set_zero!(hobo, [1,i])
    end
    println(length(used_variables(hobo)))
    hobo
end

n = 5
w = rand(n, n)
w += transpose(w) - 2*diagm(diag(w))
qubo = simplified_tsp_qubo(w)
hobo = simplified_tsp_hobo(w)

hamiltonian_qubo = bom_to_hamiltonian(qubo)
hamiltonian_hobo = bom_to_hamiltonian(hobo)
#pritnln(TravelingSalesmanExact.get_optimal_tour(w))
println(minimum(hamiltonian_qubo))
println(minimum(hamiltonian_hobo))
##
n = 5
variables = Variables("test", n)
hobo = variables[1]*variables[3] + 2*variables[1] -3 *variables[2] + 5
println(hobo)
println(used_variables(hobo))

n = 6
w = rand(n, n)
w += transpose(w) - 2*diagm(diag(w))

hobo_tsp = tsp_hobo(w)
println(used_variables(hobo_tsp))

##

n = 6
w = rand(n, n)
w += transpose(w) - 2*diagm(diag(w))

function perm_to_qubovar(v::Vector{Int})
    n = length(v)
    @assert length(unique(v)) == length(v) > 0
    @assert all(1 .<= v .<= n)
    #@assert dim(var) == 2
    #@assert indexcapacities(var)[1] == indexcapacities(var)[2] == n

    result = Set{Vector{Int}}()
    for i=1:n
        push!(result, [i,v[i]])
    end
    result
end

function perm_to_hobovar(v::Vector{Int})
    n = length(v)
    logn = ceil(Int, log(2, n))
    @assert length(unique(v)) == length(v) > 0
    @assert all(1 .<= v .<= n)

    result = Set{Vector{Int}}() # first to 0
    for t=1:n
        bits = reverse(digits(v[t]-1, base=2, pad=logn))
        for (ind, b) = filter(x->x[2]==1, collect(enumerate(bits)))
            push!(result, [t, ind])
        end
    end
    result
end

qubo = tsp_qubo(w, B=1., A=0.)
qubo_test = tsp_qubo(w, B=0., A=1.)
hobo = tsp_hobo(w, B=1., A=0.)
hobo_test = tsp_hobo(w, B=0., A=1.)

function right_value(w::Matrix{T}, v::Vector{Int}) where T<:Real
    n = length(v)
    @assert length(unique(v)) == length(v) > 0
    @assert all(1 .<= v .<= n)
    result = 0.
    for i=1:n
        result += w[v[i],v[i%n+1]]
    end
    result
end

#test qubo part
for _=1:100
    perm = Random.randperm(n)
    var_set_qubo = perm_to_qubovar(perm)
    var_set_hobo = perm_to_hobovar(perm)
    println("$perm $(right_value(w,perm))")
    println("$(qubo(var_set_qubo)) ($(qubo_test(var_set_qubo)))")
    println(var_set_hobo)
    println("$(hobo(var_set_hobo)) ($(hobo_test(var_set_hobo)))")
    @assert qubo(var_set_qubo) ≈ hobo(var_set_hobo)
end

println(hobo(Set{Vector{Int}}()))


text:
func(x) = [0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2]]
Inputs = [0.3,900.0]
fp_simple   = fixed_point(func, Inputs; Algorithm = Simple)
fp_simple.Convergence_ < 1e-10
fp_anderson = fixed_point(func, Inputs; Algorithm = Anderson)
fp_anderson.Convergence_ < 1e-10
fp_aitken   = fixed_point(func, Inputs; Algorithm = Aitken)
fp_aitken.Convergence_ < 1e-10
fp_newton   = fixed_point(func, Inputs; Algorithm = Newton)
fp_newton.Convergence_ < 1e-10
fp_VEA      = fixed_point(func, Inputs; Algorithm = VEA)
fp_VEA.Convergence_ < 1e-10
fp_SEA      = fixed_point(func, Inputs; Algorithm = SEA)
fp_SEA.Convergence_ < 1e-10
fp_MPE      = fixed_point(func, Inputs; Algorithm = MPE)
fp_MPE.Convergence_ < 1e-10
fp_RRE      = fixed_point(func, Inputs; Algorithm = RRE)
fp_RRE.Convergence_ < 1e-10

# Now trying a function with a typeswitch from Int to Float
Inputs = [1,4]
fp_simple   = fixed_point(func, Inputs; Algorithm = Simple)
fp_simple.Convergence_ < 1e-10
fp_anderson = fixed_point(func, Inputs; Algorithm = Anderson)
fp_anderson.Convergence_ < 1e-10


text:
frule(Δargs, T::Type{<:AbstractRNG}, args...) = T(args...), ZeroTangent()

function rrule(T::Type{<:AbstractRNG}, args...)
    function AbstractRNG_pullback(ΔΩ)
        return (NoTangent(), map(Returns(ZeroTangent()), args)...)
    end
    return T(args...), AbstractRNG_pullback
end

@non_differentiable Broadcast.broadcastable(::AbstractRNG)

@non_differentiable Random.randexp(::AbstractRNG)
@non_differentiable Random.randstring(::AbstractRNG)

@non_differentiable rand(::AbstractRNG, ::Random.Sampler)
@non_differentiable rand(::AbstractRNG, ::Integer...)
@non_differentiable rand(::Type{<:Real})
@non_differentiable rand(::Type{<:Real}, ::Tuple)
@non_differentiable rand(::Type{<:Real}, ::Integer...)
@non_differentiable rand(::Integer...)

# There are many different 1-3 arg methods, but not varargs
@non_differentiable rand!(::Any)
@non_differentiable rand!(::Any, ::Any)
@non_differentiable rand!(::Any, ::Any, ::Any)

@non_differentiable randexp(::Any...)

@non_differentiable randexp!(::AbstractArray)
@non_differentiable randexp!(::AbstractRNG, ::AbstractArray)

@non_differentiable randn(::Any...)

@non_differentiable randn!(::AbstractArray)
@non_differentiable randn!(::AbstractRNG, ::AbstractArray)

@non_differentiable copy(::AbstractRNG)
@non_differentiable copy!(::AbstractRNG, ::AbstractRNG)

@non_differentiable Random.default_rng()
@non_differentiable Random.default_rng(::Int)


text:
################################################################################
### Mandatory Interface
################################################################################

### Model
########################################

"""
    hopping_matrix(mc::DQMC, m::Model)

Calculates the hopping matrix \$T_{i\\sigma, j\\sigma '}\$ where \$i, j\$ are
site indices and \$\\sigma , \\sigma '\$ are flavor indices (e.g. spin indices).
The hopping matrix should also contain potential chemical potential terms on the
diagonal.

A matrix element is the hopping amplitude for a hopping process: \$j,\\sigma '
\\rightarrow i,\\sigma\$.

Regarding the order of indices, if `T[i, σ, j, σ']` is your desired 4D hopping
array, then `reshape(T, (n_sites * n_flavors, :))` is the hopping matrix.
"""
hopping_matrix(mc::DQMC, m::Model) = throw(MethodError(hopping_matrix, (mc, m)))

nflavors(m::Model) = throw(MethodError(nflavors, (m,)))

### Field
########################################

"""
    interaction_matrix_exp!(
        mc::DQMC, m::Model, field::AbstractField,
        result::AbstractArray, slice::Int, power::Float64 = 1.0
    )

Calculate the, exponentiated interaction matrix
`exp(- power * delta_tau * V(slice))` and store it in `result::AbstractArray`.

This only includes terms with 4 operators, i.e. not the chemical potential or
any hoppings. By default the calculation will be performed by the appropriate
field type (i.e. by `interaction_matrix_exp!(field, result, slice, power)`)

By default this function will call
`interaction_matrix_exp!(field, result, slice, power)`.

This is a performance critical method and one might consider efficient in-place
(in `result`) construction.
"""
@inline function interaction_matrix_exp!(
        mc, model, field, result, slice, power = +1.0
    )
    interaction_matrix_exp!(field, result, slice, power)
end
function interaction_matrix_exp!(f, A, s, p)
    throw(MethodError(interaction_matrix_exp!, (f, A, s, p)))
end

"""
    propose_local(mc::DQMC, m::Model, field::AbstractField, i::Int, slice::Int)

Propose a local move for lattice site `i` at time slice `slice` for a `field`
holding the current configuration. Returns the Green's function determinant
ratio, the boson energy difference `ΔE_boson = E_boson_new - E_boson`,
and any extra information `passthrough` that might be useful in `accept_local`.

By default this function will call `propose_local(mc, field, i, slice)`.

See also [`accept_local!`](@ref).
"""
@inline propose_local(mc, m, field, i, slice) = propose_local(mc, field, i, slice)
propose_local(mc, f, i, s) = throw(MethodError(propose_local, (mc, f, i, s)))

"""
    accept_local!(
        mc::DQMC, m::Model, field::AbstractField, i::Int, slice::Int,
        detratio, ΔE_boson, passthrough
    )

Accept a local move for site `i` at imaginary time slice `slice` of current
configuration in `field`. Arguments `detratio`, `ΔE_boson` and `passthrough`
correspond to output of `propose_local` for that local move.

By default this function will call
`accept_local!(mc, field, i, slice, detration, ΔE_boson, passthrough)`

See also [`propose_local`](@ref).
"""
@inline function accept_local!(
        mc, m, field, i, slice, detratio, ΔE_boson, passthrough
    )
    accept_local!(mc, field, i, slice, detratio, ΔE_boson, passthrough)
end
function accept_local!(mc, f, i, s, d, ΔE, pt)
    throw(MethodError(accept_local!, (mc, f, i, s, d, ΔE, pt)))
end

Base.rand(f::AbstractField) = throw(MethodError(rand, (f,)))
Random.rand!(f::AbstractField) = throw(MethodError(rand!, (f,)))

nflavors(f::AbstractField) = throw(MethodError(nflavors, (f,)))

# For ConfigRecorder
compress(f::AbstractField) = throw(MethodError(compress, (f, )))
compressed_conf_type(f::AbstractField) = throw(MethodError(compressed_conf_type, (f, )))
decompress(f::AbstractField, c) = throw(MethodError(decompress, (f, c)))
decompress!(f::AbstractField, c) = throw(MethodError(decompress!, (f, c)))

################################################################################
### Optional Interface
################################################################################

### Mixed (Field + Model)
########################################

"""
    greens_eltype(field, model)

Returns the type of the elements of the Green's function matrix. Defaults to
Float64 if both the hopping and interaction matrix contain floats and ComplexF64
otherwise.
"""
function greens_eltype(field::AbstractField, model::Model)
    generalized_eltype(interaction_eltype(field), hopping_eltype(model))
end

"""
    greens_matrix_type(field, model)

Returns the (matrix) type of the greens and most work matrices. Defaults to
`Matrix{greens_eltype(T, m)}`.
"""
greens_matrix_type(f::AbstractField, m::Model) = Matrix{greens_eltype(f, m)}

### Model
########################################

"""
    hopping_eltype(model)

Returns the type of the elements of the hopping matrix. Defaults to `Float64`.
"""
hopping_eltype(::Model) = Float64

"""
    hopping_matrix_type(field, model)

Returns the (matrix) type of the hopping matrix. Defaults to
`Matrix{hopping_eltype(model)}`.
"""
hopping_matrix_type(::AbstractField, m::Model) = Matrix{hopping_eltype(m)}

### Field
########################################

"""
    interaction_eltype(model)

Returns the type of the elements of the interaction matrix. Defaults to `Float64`.
"""
interaction_eltype(::AbstractField) = Float64

"""
    interaction_matrix_type(field, model)

Returns the (matrix) type of the interaction matrix. Defaults to
`Matrix{interaction_eltype(model)}`.
"""
interaction_matrix_type(f::AbstractField, m::Model) = Matrix{interaction_eltype(f)}

"""
    init_interaction_matrix(field::AbstractField, model::Model)

Returns an initial interaction matrix. This only used to allocate a correctly
sized matrix.

By default this uses the matrix type from `interaction_matrix_type` and uses
`max(nflavors(field), nflavors(model)) * length(lattice(model))` as the size.
"""
function init_interaction_matrix(f::AbstractField, m::Model)
    flv = max(nflavors(f), nflavors(m))
    N = length(lattice(m))
    FullT = interaction_matrix_type(f, m)

    if FullT <: BlockDiagonal
        MT = matrix_type(interaction_eltype(f))
        return BlockDiagonal(MT(undef, N, N), MT(undef, N, N))
    elseif FullT <: Diagonal
        VT = vector_type(interaction_eltype(f))
        Diagonal(VT(undef, flv * N))
    else
        FullT(undef, flv * N, flv * N)
    end
end

"""
    energy_boson(mc::DQMC, model::Model, conf)

Calculate bosonic part (non-Green's function determinant part) of energy for
configuration `conf` for Model `m`.

This is required for global and parallel updates as well as boson energy
measurements, but not for local updates. By default calls
`energy_boson(field(mc), conf)`
"""
energy_boson(mc::DQMC, m::Model, c = nothing) = energy_boson(field(mc), c)
energy_boson(f::AbstractField, c = nothing) = throw(MethodError(energy_boson, (f, c)))

conf(f::AbstractField) = f.conf
conf!(f::AbstractField, c) = conf(f) .= c
temp_conf(f::AbstractField) = f.temp_conf


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Compute the Sun velocity vector.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] The Astronomical Almanac for the year 2000 (p. C24).
#
#   [2] http://aa.usno.navy.mil/faq/docs/SunApprox.php
#
#   [3] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorne, CA.
#
#   [4] The Astronomical Almanac for the year 2006.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export sun_velocity_i

"""
    function sun_velocity_i(JD::Number)

Compute the Sun velocity represented in the Mean Equinox of Date (MOD) at the
Julian Day `JD`. The algorithm was obtained by computing the time derivative of
the Sun position in [3, p. 277-279].

"""
function sun_velocity_i(JD::Number)
    # Constants
    cen2s = 1/(36525.0*86400.0)

    # Sun position
    # ==========================================================================

    # Number of Julian centuries from J2000 epoch.
    T_UT1 = (JD-JD_J2000)/36525.0

    # Mean longitude of the Sun [deg].
    λ_m = mod(280.460 + 36000.771T_UT1, 360)

    # Mean anomaly of the Sun [deg].
    #
    # Here, we should use T_TBD (Barycentric Dynamical Time). However, it is
    # sufficient to use T_UT1 because this is a low precision computation [1].
    Ms = mod(357.5291092 + 35999.05034T_UT1, 360)

    # Convert Ms to [rad].
    Ms = deg2rad(Ms)

    # Ecliptic latitude of the Sun [deg].
    λ_e = mod(λ_m + 1.914666471sin(Ms) + 0.019994643sin(2Ms), 360)

    # Obliquity of the ecliptic [deg].
    #
    # Here, we are also approximating T_TBD ≈ T_UT1.
    ϵ = mod(23.439291 - 0.0130042T_UT1, 360)

    # Auxiliary variables.
    sin_Ms  , cos_Ms  = sincos(Ms)
    sin_2Ms , cos_2Ms = sincos(2Ms)
    sin_ϵ   , cos_ϵ   = sincos(deg2rad(ϵ))
    sin_λ_e , cos_λ_e = sincos(deg2rad(λ_e))

    # Distance of the Sun from Earth [m].
    r = ( 1.000140612 - 0.016708617cos_Ms - 0.000139589cos_2Ms )*au2m

    # Sun velocity
    # ==========================================================================

    # Compute the required time derivatives [rad/s].
    ∂λ_m = deg2rad(36000.771)   * cen2s
    ∂Ms  = deg2rad(35999.05034) * cen2s
    ∂ϵ   = deg2rad(-0.0130042)  * cen2s
    ∂r   = (+0.016708617sin_Ms*∂Ms + 2*0.000139589cos_2Ms*∂Ms)*au2m
    ∂λ_e = ∂λ_m + deg2rad(1.914666471cos_Ms*∂Ms + 2*0.019994643cos_2Ms*∂Ms)

    # Compute the Sun velocity vector represented in the Mean Equinox of Date
    # (MOD).
    vs_MOD = SVector{3,Float64}(
               ∂r*cos_λ_e -                            r*sin_λ_e*∂λ_e,
               ∂r*cos_ϵ*sin_λ_e - r*sin_ϵ*sin_λ_e*∂ϵ + r*cos_ϵ*cos_λ_e*∂λ_e,
               ∂r*sin_ϵ*sin_λ_e + r*cos_ϵ*sin_λ_e*∂ϵ + r*sin_ϵ*cos_λ_e*∂λ_e
              )
end


text:
_tostring(io) = replace.(replace.(readlines(seek(io,0)), Ref("\n" => "")), Ref("$Int" => "Int64"))

@testset "256 colors" begin
    ImageInTerminal.use_256()
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), lena)
    @test_reference "lena_show_256" _tostring(io)
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), rgb_line)
    @test_reference "rgbline_show_256" _tostring(io)
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), RGB(0.5,0.1,0.9))
    @test_reference "colorant_show_256" _tostring(io)
end

@testset "24 bit" begin
    ImageInTerminal.use_24bit()
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), lena)
    @test_reference "lena_show_24bit" _tostring(io)
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), rgb_line)
    @test_reference "rgbline_show_24bit" _tostring(io)
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), RGB(0.5,0.1,0.9))
    @test_reference "colorant_show_24bit" _tostring(io)
end

ImageInTerminal.use_256()


text:
function parareal_mpi!(cache::PararealCache, solution::PararealSolution, problem, solver::Parareal; savechunks::Bool = false)
    return error("Missing implementation.")
end


text:
module TestGroupGenerics

using Test
using RLEVectors

# Test function for broadcasting
count = 0
function simple_add(x)
    global count += 1
    x + 100
end

@testset begin

    # compare group
    vec = [1, 1, 2, 2, 4, 4]
    x = RLEVector(vec)
    @test (x .< 3) == RLEVector([true, false], [4, 6])

    # math group
    vec = [1, 1, 2, 2, 4, 4]
    x = RLEVector(vec)
    @test x .+ 5 == 5 .+ x
    @test x .+ 4 == RLEVector([5, 5, 6, 6, 8, 8])
    @test div.(vec, 2) == collect(div.(x, 2))
    @test median(RLEVector([1, 2, 3, 2, 1, 5, 4])) == median([1, 2, 3, 2, 1, 5, 4])
    @test median(RLEVector([1, 2, 3, 2, 1, 5])) == median([1, 2, 3, 2, 1, 5])
    @test median(RLEVector([3, 2, 1])) == 2.0
    @test median(RLEVector([2, 2, 2, 3, 3, 3])) == 2.5 # median is average of end of run and next value
    @test sum(RLEVector([4, 4, 5, 5, 6, 6])) == 30
    @test mean(RLEVector([4, 4, 5, 5, 6, 6])) == 5.0
    @test x .^ 2 == RLEVector([1, 1, 4, 4, 16, 16])
    @test x .^ 3 == RLEVector([1, 1, 8, 8, 64, 64])
    @test x + x == RLEVector([2, 2, 4, 4, 8, 8])
    @test x - x == RLEVector([0, 0, 0, 0, 0, 0])
    rle = RLEVector([4, 4, 9, 9, 16, 16])
    @test sqrt.(rle) == RLEVector([2, 2, 3, 3, 4, 4])

    # math on bools
    vec = [1, 1, 2, 2, 4, 4]
    x = RLEVector(vec)
    @test x .+ true == x .+ 1
    @test x .+ false == x
    vec = [1, 1, 2, 2, 4, 4]
    x = RLEVector(vec)
    @test x .+ true == x .+ 1
    @test x .+ false == x

    # broadcast with different sizes
    x = RLEVector([1, 1, 3])
    @test x .+ 4 == RLEVector([5, 5, 7])
    @test x .+ [4 5 6] == [[5 6 7]; [5 6 7]; [7 8 9]]
    @test x .+ [4, 5, 6] == RLEVector([5, 6, 9])
    @test x .+ [[2 2 2]; [4 4 4]]' == [[3 3 5]; [5 5 7]]'

    # test that broadcast only happens once per element
    x = RLEVector([1,1,1,2,2,3]) # 6 elements, 3 runs
    y = simple_add.(x)
    @test count == 6

    # findmax, findmin
    @test findmin(RLEVector([1, 2, 3, 4, 1, 1])) == findmin([1, 2, 3, 4, 1, 1])
    @test findmax(RLEVector([1, 2, 3, 4, 1, 1])) == findmax([1, 2, 3, 4, 1, 1])

    # indexin
    foo = IntegerRle(collect(1:1000), collect(5:5:5000))
    x = RLEVector([2, 2, 4, 4, 3, 3])
    y = RLEVector([0, 0, 0, 3, 3, 3, 4, 4])
    @test indexin(x, y) == indexin(collect(x), collect(y))
    @test indexin(x, collect(3:11)) == indexin(collect(x), collect(3:11))
    @test indexin([200, 200, 1, 1, 5, 5], foo) == indexin([200, 200, 1, 1, 5, 5], foo)

    @test findall(in(RLEVector(collect(3:10))), RLEVector([1, 1, 2, 2, 3, 3])) == collect(5:6)
    @test findall(in(3:10), RLEVector([1, 1, 2, 2, 3, 3])) == collect(5:6)
    @test findall(in(collect(3:10)), RLEVector([1, 1, 2, 2, 3, 3])) == collect(5:6)
    @test findall(in(RLEVector(collect(1:4))), [3, 4, 5]) == [1, 2]

    # in
    @test in(3, RLEVector([1, 2, 2, 3])) == true
    @test in(4, RLEVector([1, 2, 2, 3])) == false
    @test setdiff(Set([1, 2, 3, 4, 5]), RLEVector([1, 1, 2, 2, 4, 4, 5, 5])) == Set([3])

end # testset

end # module


text:
__precompile__(true)

module DandelionWebSockets

export AbstractWSClient,
       WSClient,
       stop,
       send_text,
       send_binary

export WebSocketHandler,
       WebSocketConnection,
       on_text,
       on_binary,
       state_closed,
       state_closing,
       state_connecting,
       state_open,
       wsconnect
export AbstractClientProtocol

abstract type AbstractWSClient end
abstract type AbstractClientProtocol end
abstract type AbstractPinger end
abstract type AbstractPonger end
abstract type ClosingBehaviour end
abstract type AbstractFrameWriter end

abstract type WebSocketConnection end

# This defines the public interface that the user should implement. These are callbacks called when
# events arrive from this WebSocket library.
abstract type WebSocketHandler end

"Handle a text frame."
on_text(t::WebSocketHandler, ::String) = nothing

"Handle a binary frame."
on_binary(t::WebSocketHandler, ::AbstractVector{UInt8}) = nothing

"The WebSocket was closed."
state_closed(t::WebSocketHandler) = nothing

"The WebSocket is about to close."
state_closing(t::WebSocketHandler) = nothing

"The WebSocket is trying to connect."
state_connecting(t::WebSocketHandler, c::WebSocketConnection) = nothing

"The WebSocket is open and ready to send and receive messages."
state_open(t::WebSocketHandler) = nothing

include("UniformResourceIdentifiers.jl")
include("WebSocketUpgrade.jl")
include("core.jl")
include("protocol_types.jl")
include("network.jl")
include("closing_behaviour.jl")
include("frame_writer.jl")
include("client_protocol.jl")
include("ping.jl")
include("handshake_logic.jl")
include("handshake.jl")
include("Proxy/Proxy.jl")
include("frame_sender.jl")
include("client.jl")
include("reconnect.jl")
include("mock.jl")

end # module


text:
# Example for loop unrolling using invertible networks
# Adapted from Putzky and Welling (2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test
using Flux, JOLI
import Flux.Optimise.update!

# Input
nx = 32
ny = 32
nz = 32
n_in = 4
n_hidden = 8
batchsize = 2
maxiter = 2

# Observed data
nxrec = 20
nyrec = 10
nt = 50
d = randn(Float32, nt, nxrec, nyrec, 1, batchsize)

# Modeling/imaging operator (can be JOLI/JUDI operator or explicit matrix)
J = joMatrix(randn(Float32, nt*nxrec*nyrec, nx*ny*nz))

# Link function
Ψ(η) = identity(η)

# Unrolled loop
L = NetworkLoop3D(n_in, n_hidden, maxiter, Ψ; type="HINT")

# Initializations
η_obs = randn(Float32, nx, ny, nz, 1, batchsize)
s_obs = randn(Float32, nx, ny, nz, n_in-1, batchsize)
η_in = randn(Float32, nx, ny, nz, 1, batchsize)
s_in = randn(Float32, nx, ny, nz, n_in-1, batchsize)

# Forward pass and residual
η_out, s_out = L.forward(η_in, s_in, d, J)
Δη = η_out - η_obs
Δs = s_out - s_obs  # in practice there is no observed s, so Δs=0f0

# Backward pass
Δη_inv, Δs_inv, η_inv, s_inv = L.backward(Δη, Δs, η_out, s_out, d, J)

# Check invertibility
@test isapprox(norm(η_inv - η_in)/norm(η_inv), 0f0, atol=1e-5)
@test isapprox(norm(s_inv - s_in)/norm(s_inv), 0f0, atol=1e-5)

# Update using Flux optimizer
opt = Flux.ADAM()
P = get_params(L)
for p in P
    update!(opt, p.data, p.grad)
end
clear_grad!(L)


text:
#NetwAnalysis show functions
#-------------------------------------------------------------------------------

#==Generate friendly show functions
===============================================================================#

Base.show(io::IO, t::Type{NPType{NP}}) where NP = print(io, "NPType{:$NP}")

_printmatrix(io::IO, m::NetworkParameterMatrix{T}) where {T<:Number} = println(io, m)
function _printmatrix(io::IO, m::NetworkParameterMatrix)
	nports = portcount(m)
	for row = 1:nports, col = 1:nports
		print(io, "  ($row,$col)=")
		if isassigned(m, row, col)
			println(io, m[row, col])
		else
			println(io, "UNDEFINED")
		end
	end
end

function Base.show(io::IO, np::SParameters)
	println(io, "SParameters{$(portcount(np))-port, z0=$(np.z0)}[")
	_printmatrix(io, np.m)
	println(io, "]")
end

function Base.show(io::IO, np::TParameters)
	println(io, "TParameters{$(portcount(np))-port, z0=$(np.z0)}[")
	_printmatrix(io, np.m)
	println(io, "]")
end

function Base.show(io::IO, np::ZParameters)
	println(io, "ZParameters{$(portcount(np))-port}[")
	_printmatrix(io, np.m)
	println(io, "]")
end

function Base.show(io::IO, np::YParameters)
	println(io, "YParameters{$(portcount(np))-port}[")
	_printmatrix(io, np.m)
	println(io, "]")
end

function Base.show(io::IO, np::ABCDParameters)
	println(io, "ABCDParameters[")
	_printmatrix(io, np.m)
	println(io, "]")
end

#Last line


text:
# module IRI2016
# using Dates
# using RemoteFiles
# # import PlotlyJS
# isapple = Sys.isapple()
# islinux = Sys.islinux()

const path_iri_offline = string(@__DIR__, "/offline_2021_01_16")
const _iri2016_so_filename =
            string(
                    @__DIR__, "/iri2016",
                    isapple ? "_macos" : islinux ? "_linux" : "_win",
                    ".so",
                )

const _iri_tar = @RemoteFile(
        "http://irimodel.org/IRI-2016/00_iri.tar",
        file = "00_iri.tar",
        dir = @__DIR__,
        updates = :never
    )

const _ig_rz_dat = @RemoteFile(
        "https://chain-new.chain-project.net/echaim_downloads/ig_rz.dat",
        file = "ig_rz.dat",
        dir = @__DIR__,
        updates = :daily
        )

const _apf107_dat = @RemoteFile(
        "https://chain-new.chain-project.net/echaim_downloads/apf107.dat",
        file = "apf107.dat",
        dir = @__DIR__,
        updates = :daily
        )

const _ccir_ursi_tar = @RemoteFile(
            "http://irimodel.org/COMMON_FILES/00_ccir-ursi.tar",
            file = "00_ccir-ursi.tar",
            dir = @__DIR__,
            updates = :never
        )

function _check_iri2016()
    iri_tar = download_or_cp_file(_iri_tar, path_iri_offline, RemoteFiles.Http())
    run(`tar -xvf $iri_tar --directory $(dirname(iri_tar))`)
    _check_so(_iri2016_so_filename)

end
_check_iri2016()
# IRI_SUB

# hwm14(IYD::Int32, SEC::Float32,
#     ALT::Float32, GLAT::Float32, GLONG::Float32,
#     STL::Float32, F107A::Float32, F107::Float32,
#     AP::Array{Float32,1},
#     W::Array{Float32,1},
# ) = ccall((:hwm14_, _hwm14_so_filename),
#     Nothing,
#     (Ref{Int32},   # IYD
#         Ref{Float32}, # SEC
#         Ref{Float32}, # ALT
#         Ref{Float32}, # GLAT
#         Ref{Float32}, # GLONG
#         Ref{Float32}, # STL
#         Ref{Float32}, # F107A
#         Ref{Float32}, # F107
#         Ptr{Float32}, #  AP
#         Ptr{Float32},),
#     IYD, SEC,
#     ALT, GLAT, GLONG,
#     STL, F107A, F107,
#     AP,
#     W,
# )

function _init_iri()
    download_or_cp_file(_ig_rz_dat, path_iri_offline, )
    download_or_cp_file(_apf107_dat, path_iri_offline, )

    ccir_ursi_tar = path(_ccir_ursi_tar)
    if !isfile(ccir_ursi_tar)
        download_or_cp_file(_ccir_ursi_tar, path_iri_offline, RemoteFiles.Http())
        run(`tar -xvf $ccir_ursi_tar --directory $(dirname(ccir_ursi_tar))`)
    end

    this_pwd = pwd()
    if this_pwd != @__DIR__
        cd(@__DIR__)
        @warn "work directory (iri2016) have changed\n from: $this_pwd\n   to: $(@__DIR__)"
        # c = true
    end

    ccall((:read_ig_rz_, _iri2016_so_filename), Nothing, ())
    ccall((:readapf107_, _iri2016_so_filename), Nothing, ())
    this_pwd
    # if c cd(thispwd) end
end

_iri_sub(JF::Array{Int32,1},JMAG::Int32,
            ALATI::Float32,ALONG::Float32,
            IYYYY::Int32,MMDD::Int32,DHOUR::Float32,
            HEIBEG::Float32,HEIEND::Float32,HEISTP::Float32,
            OUTF::Array{Float32,2},OARR::Array{Float32,1}) =
            ccall((:iri_sub_, _iri2016_so_filename),
                Nothing,
                (Ptr{Int32}, Ref{Int32}, Ref{Float32}, Ref{Float32}, # jf,JMAG, lat lon
                Ref{Int32}, Ref{Int32}, Ref{Float32}, # IYYYY MMDD DHOUR
                Ref{Float32}, Ref{Float32}, Ref{Float32}, # altbeg altend altstp
                Ptr{Float32}, Ptr{Float32},),
                JF,JMAG,ALATI,ALONG,IYYYY,MMDD,DHOUR,HEIBEG,HEIEND,HEISTP,OUTF,OARR
            )

function test_iri()
    # (JF, JMAG, ALATI, ALONG, IYYYY, MMDD, DHOUR, HEIBEG, HEIEND, HEISTP, )

    _init_iri()
    dn = DateTime(2014, 3, 23, 14, 30)

    # jf[3] = 0  # 4 B0,B1 other model-31
    # jf[4] = 0  # 5  foF2 - URSI
    # jf[5] = 0  # 6  Ni - RBV-10 & TTS-03
    # jf[20] = 0  # 21 ion drift not computed
    # jf[22] = 0  # 23 Te_topside (TBT-2011)
    # jf[27] = 0  # 28 spreadF prob not computed
    # jf[28] = 0  # 29 (29,30) => NeQuick
    # jf[29] = 0  # 30
    # # (Brian found a case that stalled IRI when on):
    # jf[32] = 0  # 33 Auroral boundary model off
    # jf[34] = 0  # 35 no foE storm update
    # # Not standard, but outputs same as values as standard so not an issue
    # jf[21] = 0  # 22 ion densities in m^-3 (not %)
    # jf[33] = 0  # 34 turn messages off
    # [3,4,5,20,22,27,28,29,32,34, 21,33] .+ 1

    JF = ones(Int32, 50)
    for ind in [4, 5, 6, 21, 23, 28, 29, 30, 33, 35, 22, 34]
        JF[ind] = 0
    end

    # (datetime.datetime(2014, 3, 23, 14, 30), 40.0, -80.0, 250.0)

    JMAG = Int32(0)
    ALATI = Float32(40.0)
    ALONG = Float32(-80.0)
    IYYYY = Int32(year(dn))
    MMDD = -Int32(dayofyear(dn))
    DHOUR = Float32(hour(dn) * 3600.0 + minute(dn) * 60) / 3600 + 25

    HEIBEG  = Float32(100.0)
    HEIEND = Float32(1000.0)
    HEISTP = Float32(1.0)

    OUTF = Array{Float32,2}(undef, (20, 1000))
    OARR = Array{Float32,1}(undef, 100)
    _iri_sub(JF, JMAG, ALATI, ALONG, IYYYY, MMDD, DHOUR, HEIBEG, HEIEND, HEISTP, OUTF, OARR)
    # @show OUTF[1:10,1]
    # tr = PlotlyJS.scatter(;y=OUTF[1,:])
    # PlotlyJS.savehtml(PlotlyJS.plot(tr, ), "test.html")
    # println(OUTF[1]) 1.0536364e11
    OUTF, OARR
end
# iri()

# end moulde IRI2016
# end


text:
using Miletus, Blpapi
IP = "localhost"; Port = 8194
Session = createSession(IP, Port)
methods(bdp)
Response = bdp(Session, tickers, fields)
ibmLastPrice = Response["IBM US Equity", "PX_Last"]
ibmOpenPrice = Response["IBM US Equity", "PX_Open"]
appleLastPrice = Response["AAPL US Equity", "PX_Last"]
appleOpenPrice = Response["AAPL US Equity", "PX_Open"]
tickers =  ["IBM US Equity", "AAPL US Equity"]
fields = ["PX_Open", "PX_High", "PX_Last", "PX_Close"]

startDate = "20070707"
endDate = "20170707"

history = bdh(session, tickers, fields, startDate, endDate)
using JuliaDB
using Miletus, BusinessDays
using Miletus.TermStructure
using Miletus.DayCounts
import Miletus: Both, Receive, Contract, When, AtObs, value
import Miletus: YieldModel
import BusinessDays: USGovernmentBond
import Base.Dates: today, days, Day, Year

using Miletus
using Gadfly
using Colors
d1 = Dates.today()
d2 = d1 + Dates.Day(120)

m = GeomBMModel(d1, 100.00, 0.05, 0.0, 0.3)
mcm = montecarlo(m, d1:d2, 100_000)

theme=Theme(default_color=Colors.RGBA{Float32}(0.1, 0.1, 0.7, 0.1))
p = plot([layer(x=mcm.dates,y=mcm.paths[i,:],Geom.line,theme) for i = 1:200]...)

o = EuropeanCall(d2, SingleStock(), 100.00)
value(m, o)

value(mcm, o)

oa1 = AsianFixedStrikeCall(d2, SingleStock(), Dates.Month(1), 100.00)

using JLD2

file = "Blpapi"
pwd()
@save "blpapihistory.jld2" file history

pwd()

max_stars_repo_path: Dojo/C++ Programming Practices and Principles - Stroustrup/trialblpapi.jl

max_stars_repo_name: miguelraz/PathToPerformance

max_issues_repo_name: miguelraz/PathToPerformance

max_forks_repo_name: miguelraz/PathToPerformance


text:
# Built in operators, functions, and effects

"Mapping from PDDL data types to Julia types."
const GLOBAL_DATATYPES = Dict{Symbol,Type}(
    :boolean => Bool,
    :integer => Int,
    :numeric => Float64,
)

"Mapping from PDDL comparison operators to Julia functions."
const GLOBAL_PREDICATES = Dict{Symbol,Function}(
    op => eval(op) for op in [:<=, :>=, :<, :>]
)

"Mapping from PDDL evaluation operators to Julia functions."
const GLOBAL_FUNCTIONS = Dict{Symbol,Function}(
    op => eval(op) for op in [:+, :-, :*, :/]
)

"Mapping from PDDL modifiers (i.e. in-place assignments) to Julia functions."
const GLOBAL_MODIFIERS = Dict{Symbol,Symbol}(
    :increase => :+,
    :decrease => :-,
    Symbol("scale-up") => :*,
    Symbol("scale-down") => :/
)

"""
    equiv(a, b)

Equivalence in concrete and abstract domains. Defaults to `a == b`.
"""
equiv(a, b) = a == b
GLOBAL_PREDICATES[:(==)] = equiv

"""
    nequiv(a, b)

Non-equivalence in concrete and abstract domains. Defaults to `a != b`.
"""
nequiv(a, b) = a != b
GLOBAL_PREDICATES[:(!=)] = nequiv

merge!(GLOBAL_FUNCTIONS, GLOBAL_PREDICATES)


text:
"""`meanvalueform` performs the second order extension on the given function `f` with its optional derivative
`f_prime` and initial intervals `x` and  its midpoints `x_mid`."""

function meanvalueform(f::Function, f_prime::Function, x::IntervalBox, x_mid::IntervalBox)

    mvf = f(x_mid) + dot(f_prime(x), (x - x_mid).v)

    return mvf

end

# use automatic differentiation if no derivative function given:

meanvalueform(f::Function, x::IntervalBox, x_mid::IntervalBox) = meanvalueform(f, x->∇(f,x.v), x, x_mid)

# determines the midpoints of x if x_mid is not given:

meanvalueform(f::Function, f_prime, x::IntervalBox) = meanvalueform(f, f_prime, x, IntervalBox(mid(x)))

# use automatic differentiation if no derivative function given and determines the midpoints of x if x_mid is not given:

meanvalueform(f::Function, x::IntervalBox) = meanvalueform(f, x->∇(f,x.v), x, IntervalBox(mid(x)))

"""`smear` performs the election of the variable with the largest value of semar like direction to bisect on the given function `f`
with its optional derivative `f_prime` and initial intervals `x`."""

function smear(f::Function, f_prime::Function, x::IntervalBox)

    smear_magnitude = abs.(f_prime(x)) .* diam.(X)

    return argmax([X.hi for X in ub])

end

# use automatic differentiation if no derivative function given:

smear(f::Function, x::IntervalBox) = smear(f, x->∇(f,x.v), x)

"""`localoptim` performs a local optimization of a given function `f` on intervals `x`."""

function localoptim(f::Function, X::StaticArray)

    xc = X
    gamma = 0.01
    max_iters = 10

    for i in 1:max_iters

        x = xc
        xc = x - gamma * ∇(f, x)

        if f(xc) > f(x)

            break

        end

    end

    return (f(xc), xc)

end

"""`random` generates random numbers over intervals `x`."""

randnum(r::SArray, x::IntervalBox) = (r * (x.hi - x.lo)) + x.lo

function random(X::IntervalBox{N,T}, n::Int64) where {N,T}

    Random_numbers = typeof(mid(X))[ ]

    for i in 1:n

        rand_numb = @SVector rand(N)

        interval_rand_numb = randnum.(rand_numb, X)

        push!(Random_numbers, interval_rand_numb)

    end

    return Random_numbers

end


text:
using LaTeXStrings, Compat
using Base.Test

# write your own tests here
tst1 = L"an equation: $\alpha^2$"
tst1u = Compat.UTF8String(tst1)
@test tst1 == "an equation: \$\\alpha^2\$" == tst1u
@test L"\alpha^2" == "\$\\alpha^2\$"
@test stringmime("text/latex", tst1) == tst1u
buf = IOBuffer()
@compat show(buf, "text/latex", tst1)
@test takebuf_string(buf) == tst1u
@test latexstring(tst1u) == tst1 == LaTeXString(tst1u)

@test ccall(:strlen, Csize_t, (Cstring,), tst1) == ccall(:strlen, Csize_t, (Ptr{UInt8},), tst1) == sizeof(tst1)


text:
# Single run of model with 2 parameters
using Turing, Statistics, Distributions, Random, LinearAlgebra, StatsPlots
using AdvancedMH

# the data generating process
function dgp(θ, n)
    [rand.(Exponential(θ[1]),n) rand.(Poisson(θ[2]),n)]
end

# summary statistics for estimation
function moments(y)
    n = size(y,1)
    sqrt(n) .* [mean(y, dims=1)[:]; std(y, dims=1)[:]]
end

@model function abc(z, S, n, calibrated)
    # create the prior: the product of the following array of marginal priors
    θ  ~ arraydist([LogNormal(1.,1.); LogNormal(1.,1.)])
    # sample from the model, at the trial parameter value, and compute statistics
    y = zeros(n,2)
    zs = zeros(S, size(z,1))
    @inbounds for i = 1:S
        y .= dgp(θ, n)
        zs[i,:] .= moments(y) # simulated summary statistics
    end
    # the asymptotic Gaussian distribution of the statistics
    m = mean(zs, dims=1)[:]
    calibrated ? Σ = Symmetric((1. + 1/S)*cov(zs)) : Σ = I
    z ~ MvNormal(m, Σ)
end;

function main(calibrated=true)
    θ⁰ = [2.; 3.] # true parameters
    n = 100 # sample size
    S = 100 # number of simulation draws
    # get data and statistics
    y = dgp(θ⁰, n)
    z = moments(y)
    # sample chains
    length = 5000
    burnin = 500
    chain = sample(abc(z, S, n, calibrated),
        MH(:θ => AdvancedMH.RandomWalkProposal(MvNormal(zeros(2), 0.25*I))),
        MCMCThreads(), length+burnin, 4)
    chain = chain[burnin+1:end,:,:]
    @show display(chain)
    p = plot(chain)
    chain = Array(chain)
    q1 = quantile(chain[:,1], [0.025, 0.975])
    q2 = quantile(chain[:,2], [0.025, 0.975])
    println()
    println("is first true parameter inside 95% CI?")
    display(q1[1] <= θ⁰[1] && q1[2] >= θ⁰[1])
    println()
    println("is second true parameter inside 95% CI?")
    display(q2[1] <= θ⁰[2] && q2[2] >= θ⁰[2])
    p
end
println("first a well-calibrated run")
p1 = main(true)
println("now a poorly calibrated run")
p2 = main(false)
plot(p1, p2)

max_stars_repo_name: mcreel/ABC-MSM-Turing

max_issues_repo_name: mcreel/ABC-MSM-Turing

max_forks_repo_name: mcreel/ABC-MSM-Turing


text:
module DomainIntegrals

using QuadGK, LinearAlgebra
using GaussQuadrature, FastGaussQuadrature, HCubature

using StaticArrays
using IntervalSets, DomainSets
using CompositeTypes.Display

import Base: convert, checkbounds, length,
    *, ∘

import DomainSets:
    domaintype,
    codomaintype,
    prectype,
    numtype,
    forward_map,
    components

using DomainSets: ×,
    EmptySpace, FullSpace,
    euclideandimension

include("extra/common.jl")

include("gauss.jl")
include("functional/measure.jl")
include("weights.jl")
include("singularity.jl")
include("strategy.jl")
include("integrand.jl")
include("processing/properties.jl")
include("processing/measure.jl")
include("processing/domain.jl")
include("integral.jl")
include("complexplane.jl")

include("functional/forms.jl")

end # module


text:
# smc utilities

function uniqueParticles!(unq::Array{Int64,2},t::Int64,ed_idx::Array{Int64,1},ancestors::Array{Int64,2})

  nunq = zero(Int64)
  anc = ancestors[t,:]
  for n = 1:length(ed_idx)
    if t==1
      root_pos = findfirst( ed_idx .== ed_idx[n] )
    else
      root_pos = findfirst( ( ed_idx .== ed_idx[n] ).*( unq[t-1,anc] .== unq[t-1,anc[n]] ) )
    end

    if root_pos == n
      nunq += one(Int64)
      unq[t,n] = nunq
    else
      unq[t,n] = unq[t,root_pos]
    end
  end

end

function randomWalkProbs!(W::Array{Float64,2},nv::Int64,degrees::Array{Float64,1},eig_pgf::Array{Float64,1},esys_vec::Union{Array{Float64,2},SparseMatrixCSC{Float64,Int64}})
  s = zeros(Float64,1)
  for n = 1:nv # column index
    for m = 1:n # row index
      s[1] = zero(Float64)
      for k = 1:nv
        s[1] += esys_vec[k,m] * esys_vec[k,n] * eig_pgf[k]
      end
      W[m,n] = s[1]*sqrt(degrees[n]/degrees[m])
      n==m ? nothing : W[n,m] = s[1]*sqrt(degrees[m]/degrees[n])
    end

  end
  elMax!(W,zero(Float64)) # for numerical stability
  elMin!(W,one(Float64)) # for numerical stability
end

function randomWalkProbs!(w::Array{Float64,2},vtx_pairs::Array{Int64,2},nv::Int64,degrees::Array{Float64,1},
    eig_pgf::Array{Float64,1},esys_vec::Union{Array{Float64,2},SparseMatrixCSC{Float64,Int64}})
"""
  Calculate random walk probabilities for vertex pairs in `vtx_pairs`
"""

  # w = zeros(Float64,size(vtx_pairs,1),2)
  resetArray!(w)
  for i in 1:size(vtx_pairs,1)
    for k in 1:nv
      w[i,1] += (esys_vec[k,vtx_pairs[i,1]] * esys_vec[k,vtx_pairs[i,2]] * eig_pgf[k])::Float64
    end

    w[i,2] = w[i,1]*(degrees[vtx_pairs[i,1]]/degrees[vtx_pairs[i,2]])^(0.5)::Float64
    w[i,1] *= ((degrees[vtx_pairs[i,2]]/degrees[vtx_pairs[i,1]])^(0.5))::Float64

    elMax!(w,zero(Float64)) # for numerical stability
    elMin!(w,one(Float64)) # for numerical stability

  end

end

function randomWalkProbs(root_vtx::Int64,neighbors::Array{Int64,1},nv::Int64,degrees::Array{Float64,1},eig_pgf::Array{Float64,1},esys_vec::Union{Array{Float64,2},SparseMatrixCSC{Float64,Int64}})::Float64
"""
  Calculate random walk probabilities for vertex pairs in `vtx_pairs`
"""

  w = zero(Float64)
  z = zeros(Float64,1)
  for i in 1:convert(Int64,degrees[root_vtx])
    z[:] = zero(Float64) # reset z
    for k in 1:nv
      z[1] += (esys_vec[k,root_vtx] * esys_vec[k,neighbors[i]] * eig_pgf[k])::Float64
    end
    z[1] *= sqrt(degrees[neighbors[i]]/degrees[root_vtx])::Float64
    w += z[1]::Float64
  end
  w = max(w,zero(Float64)) # for numerical stability
  w = min(w,one(Float64)) # for numerical stability
  return w

end

function fruitlessRWProb(W::Array{Float64,2},nbd_list::Array{SparseVector{Int64,Int64},1},root_vtx::Int64)::Float64

  p = zero(Float64)
  p += W[root_vtx,root_vtx]::Float64

  rnz = nbd_list[root_vtx].nzval
  for i in rnz
    p += W[root_vtx,i]
  end
  return p
end

function fruitlessRWProb(W::Array{Float64,2},nbd_list::Array{Int64,2},root_vtx::Int64)::Float64

  p = zero(Float64)
  p += W[root_vtx,root_vtx]::Float64

  i = 1
  while nbd_list[root_vtx,i] > zero(Float64)
    p += W[root_vtx,nbd_list[root_vtx,i]]
    i += 1
  end

  return p
end

function fruitlessRWProb(W::Array{Float64,2},adj::Array{Int64,2},root_vtx::Int64,nv::Int64)::Float64

  # find neighborhood of root_vtx
  # nbd = find(adj[:,root_vtx])::Array{Float64,1}

  p = zero(Float64)
  p += W[root_vtx,root_vtx]::Float64

  for i = 1:nv
    adj[i,root_vtx]==one(Float64) ? p += W[root_vtx,i] : nothing
  end
  return p
end

function nbPred!(pgf::Array{Float64,1},nv::Int64,a_lambda::Float64,b_lambda::Float64,evals::Union{Array{Float64,1},SparseVector{Float64,Int64}})

    for i=1:nv
      pgf[i] = ((1.0 - b_lambda)^(a_lambda))*( (1.0 .- evals[i])*((1.0 - b_lambda*(1.0 - evals[i])).^(-a_lambda)) )
      isinf(pgf[i]) ? pgf[i] = zero(Float64) : nothing # Moore-Penrose pseudo-inverse
    end

end


text:
import TPT
using Base.Test

using Optim

ahs = TPT.AHSSystem(ρ = 1.0, η = 0.45)
σ = ahs.σ[1]
g = TPT.prdf(ahs)[1,1]

res = optimize(r -> -g(r), 0.9σ, 1.1σ, rel_tol=1e-16)
σ_max = Optim.minimizer(res)

@testset "Unary AHS" begin
  @test σ ≈ σ_max
end


text:
"""
    mutable struct LocalImage
        CreatedAt::Union{Missing, String} = missing
        Digest::String
        ID::String
        Repository::String
        Size::Union{Missing, String} = missing
        Tag::String
        exists::Bool = true
    end

Represents a docker image on the local machine, and stores associated metadata. Should not be
instantiated directly. If `exists` is `true`, then the image is assumed to exit and so should be
visible from utilities such as `docker image ls`.
"""
@with_kw mutable struct LocalImage <: LambdaComponent
  CreatedAt::Union{Missing, String} = missing
  Digest::String
  ID::String
  Repository::String
  Size::Union{Missing, String} = missing
  Tag::String
  exists::Bool = true
end
StructTypes.StructType(::Type{LocalImage}) = StructTypes.Mutable()
Base.:(==)(a::LocalImage, b::LocalImage) = a.ID[1:docker_hash_limit] == b.ID[1:docker_hash_limit]

function get_image_full_name(image::LocalImage)::String
  image.Repository
end

"""
    get_local_image(
      identity::String,
    )::Union{Nothing, LocalImage}

Returns a `LocalImage` object, representing a locally-stored docker image.

The passed `identity` string may be the repository name, or the docker image ID. If the image ID,
it must be at least four characters in length.
"""
function get_local_image(identity::AbstractString)::Union{Nothing, LocalImage}
  all = get_all_local_images()
  index = findfirst(all) do li
    identity_matches = li.Repository == identity
    name_matches = get_lambda_name(li) == identity
    id_matches = begin
       check_len = minimum([length(li.ID), length(identity)])
       (li.ID[begin:check_len] == identity[begin:check_len]) && (check_len >= 4)
    end
    identity_matches || name_matches || id_matches
  end
  isnothing(index) ? nothing : all[index]
end

"""
    get_all_local_images(;
        args::Vector{String} = Vector{String}(),
        jot_generated_only::Bool = true,
      )::Vector{LocalImage}

Returns a vector of `LocalImage`s, representing all locally-stored docker images.

`args` are passed to the call to `docker image ls`, that is used to populate this vector.
`jot_generated_only` specifies whether to filter for jot-generated images only.
"""
function get_all_local_images(;
    args::Vector{String} = Vector{String}(),
    jot_generated_only::Bool = true,
  )::Vector{LocalImage}
  docker_output = readchomp(`docker image ls $args --digests --filter "dangling=false" --format '{{json .}}'`)
  local_images = parse_docker_ls_output(LocalImage, docker_output)
  jot_generated_only ? filter(is_jot_generated, local_images) : local_images
end

"""
    delete!(
      image::LocalImage;
      force::Bool=false,
    )

Deletes a locally-stored docker image. The LocalImage instance continues to exist, but has its
`exists` attribute set to `false`.
"""
function delete!(image::LocalImage; force::Bool=false)
  image.exists || error("Image does not exist")
  args = force ? ["--force"] : []
  run(`docker image rm $(image.ID) $args`)
  image.exists = false
end

function is_lambda(image::LocalImage)::Bool
  occursin(".amazonaws.com/", image.Repository)
end

function get_aws_id(image::LocalImage)::String
  split(image.Repository, '.')[1]
end

function get_aws_region(image::LocalImage)::String
  split(image.Repository, '.')[4]
end

function get_aws_config(image::LocalImage)::AWSConfig
  AWSConfig(get_aws_id(image), get_aws_region(image))
end