text:
# Model LMU-40X-NUV  : https://www.thorlabs.com/newgrouppage9.cfm?objectgroup_id=3271

lmu_40x_nuv  = LaserLab.Objective("LMU-40X-NUV", 5.3mm, 5.1mm, 40.0)
lmm40xf_uvv  = LaserLab.Objective("LMM40XF-UVVV", 5.0mm, 5.1mm, 40.0)
df  = LaserLab.load_df_from_csv("/Users/pabloherrero/sabat/LaserLab/data/objectives", "LMU-40X-NUV.csv", LaserLab.enG)
wl  = 200.0:20:800.0
fdf = LaserLab.dftof(wl, df, "Transmission(%)", 0.0, 0.01)

function lmu_40x_nuv_transmission(λ::Unitful.Length)
    lnm = uconvert(nm, λ)
    fdf(lnm/nm)
end

# response of the objective is very flat from 200 nm to 1000 nm, around 85 %
# the approximation of taking the function to be a constant is sufficient for now.
function lmm40xf_uvv_transmission()
    return 0.85
end

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


text:
"""
    get_query_args_str(args::AbstractDict)

Get the mutation arguments string and a dictionary containing the argument names, which will
not necessarily match the name of the types in the schema.

Uses an internal function to go through `args` dictionary, recursing down through the values
that are  `AbstractDict` or `Vector{<:AbstractDict}`, ensuring that no argument name is
duplicated and enumerating multiple of the same type if in a vector.

# Examples
```jldoctest; setup=:(using GraphQLClient: get_query_args_str, prettify_query)
julia> args = Dict("TopLevelVarA" => "value", "TopLevelVarB" => 2);

julia> str, arg_names = get_query_args_str(args);

julia> println(prettify_query(str))
(
    TopLevelVarA:\$TopLevelVarA
    TopLevelVarB:\$TopLevelVarB
)

julia> args = Dict("A" => "value", "B" => 2, "C" => Dict("D" => 3, "E" => 4));

julia> str, arg_names = get_query_args_str(args);

julia> print(prettify_query(str))
(
    B:\$B
    A:\$A
    C:{
        D:\$D
        E:\$E
    }
)

julia> args = Dict("A" => "value", "B" => 2, "C" => [Dict("D" => 3, "E" => 4), Dict("D" => 5, "E" => 6), Dict("D" => 7, "E" => 8)]);

julia> str, arg_names = get_query_args_str(args);

julia> print(prettify_query(str))
(
    B:\$B
    A:\$A
    C:[
        {
            D:\$D1
            E:\$E1
        }
        {
            D:\$D2
            E:\$E2
        }
        {
            D:\$D3
            E:\$E3
        }
    ]
)
```
"""
function get_query_args_str(args::AbstractDict)
    arg_names = initialise_arg_names(args)
    fieldname_tracker = String[] # Used by internal function

    function recursive_get(arg_values::AbstractDict, arg_names::AbstractDict, enumerate)
        str = ""
        for (fieldname, val) in arg_values
            if isa(val, AbstractDict)
                # Wrap in fieldname:{} and go into next dictionary level
                str *= "$fieldname:{" * recursive_get(val, arg_names[fieldname], enumerate) * "},"
            elseif isa(val, AbstractVector{<:AbstractDict})
                # Wrap in fieldname:[] and to into each element (dictionary) of vector, enumerating argument names
                str *= "$fieldname:["
                for i in eachindex(val)
                    str *= "{" * recursive_get(val[i], arg_names[fieldname][i], enumerate+=1) * "},"
                end
                str *= "],"
            else
                name = fieldname
                str *= "$fieldname:\$$fieldname"
                if is_field_name_used(name, fieldname_tracker)
                    enumerate += 1
                end
                if enumerate > 0
                    name = fieldname * "__" * string(enumerate)
                    str *= string(enumerate)
                end

                # Add to tracker
                push!(fieldname_tracker, name)
                arg_names[fieldname] = name
                str *= ","
            end
        end
        return str
    end
    query_args_str = "(" * recursive_get(args, arg_names, 0) * ")"
    return query_args_str, arg_names
end

"""
    is_field_name_used(fieldname, fieldnamelist)

True if `fieldname` is in `fieldnamelist`.
"""
is_field_name_used(fieldname, fieldnamelist) = fieldname in fieldnamelist

"""
    initialise_arg_names(args::Dict)

Recursively create `arg_names` Dict with same structure (dictionaries and vectors) and keys
as `args`. Values are initialised as empty strings.

# Examples
```jldoctest; setup=:(using GraphQLClient: initialise_arg_names)
julia> args = Dict("int" => 1, "vector" => [1, 2, 3], "dict" => Dict("float" => 1.0), "vec_dict" => [Dict("bool" => true)])
Dict{String, Any} with 4 entries:
  "int"      => 1
  "dict"     => Dict("float"=>1.0)
  "vec_dict" => Dict{String, Bool}[Dict("bool"=>1)]
  "vector"   => [1, 2, 3]

julia> arg_names = initialise_arg_names(args)
Dict{String, Any} with 4 entries:
  "int"      => ""
  "dict"     => Dict{String, Any}("float"=>"")
  "vec_dict" => Dict{String, Any}[Dict("bool"=>"")]
  "vector"   => ""
```
"""
initialise_arg_names(args::Dict) = _initialise_arg_names(args)
function _initialise_arg_names(args::Dict)
    arg_names = Dict{String, Any}()
    for (key, val) in args
        arg_names[key] = _initialise_arg_names(val)
    end
    return arg_names
end
_initialise_arg_names(val) = ""
_initialise_arg_names(val::AbstractVector{<:AbstractDict}) = _initialise_arg_names.(val)

"""
    directly_write_query_args(query_args)

Returns a string which represents `query_args` by inserting values directly into the string.
Strings are bracketed with \" as this is required by GraphQL, but unfortunately this does
not work for ENUMs as we cannot tell if a string is meant to be an enum or not without
introspection. Therefore, any enums must be wrapped by a `GQLEnum`.

# Examples
```jldoctest; setup=:(using GraphQLClient: directly_write_query_args, prettify_query)
julia> query_args = Dict("string"=>"my_string", "dict"=>Dict("bool"=>true,"int"=>1), "vec"=>[Dict("float"=>1.0)])
Dict{String, Any} with 3 entries:
  "dict"   => Dict{String, Integer}("int"=>1, "bool"=>true)
  "string" => "my_string"
  "vec"    => [Dict("float"=>1.0)]

julia> println(prettify_query(directly_write_query_args(query_args)))
dict:{
    int:1
    bool:true
}
string:"my_string"
vec:[
    {
        float:1.0
    }
]

julia> query_args = Dict("string"=>"my_string", "enum" => GQLEnum("my_enum"))
Dict{String, Any} with 2 entries:
  "string" => "my_string"
  "enum"   => GQLEnum("my_enum")

julia> println(prettify_query(directly_write_query_args(query_args)))
string:"my_string"
enum:my_enum
```
"""
function directly_write_query_args(query_args)
    str = ""
    for (key, val) in query_args
        str *= "$key:"
        str *= _get_val_str(val)
        str *= ","
    end
    return rstrip(str, ',')
end

"""
    _get_val_str(val::String)
    _get_val_str(val::Dict)
    _get_val_str(val::Vector)
    _get_val_str(enum::GQLEnum)
    _get_val_str(val)

Return value string for `directly_write_query_args`
"""
_get_val_str(val::String) = "\"$val\""
_get_val_str(val::Dict) = "{" * directly_write_query_args(val) * "}"
_get_val_str(val::Vector) = "[" * join(_get_val_str.(val)) * "]"
_get_val_str(enum::GQLEnum) = "$(enum.value)"
_get_val_str(val) = string(val)

"""
    prettify_query(str)

Prettifies a string, nesting for bracket types and adding new lines for commas.

# Examples
```jldoctest; setup=:(using GraphQLClient)
julia> println(GraphQLClient.prettify_query("[{(a)}]"))
[
    {
        (
            a
        )
    }
]
julia> println(GraphQLClient.prettify_query("[{(a,b),(c,d)}]"))
[
    {
        (
            a
            b
        )
        (
            c
            d
        )
    }
]
```
"""
function prettify_query(str)
    buf = IOBuffer()
    indent = 0
    opening_chars = ('{', '(', '[')
    closing_chars = ('}', ')', ']')
    new_line_chars = ('}', ')', ']', ',')
    previous_char = 'a'
    for char in str
        # if !(char == ',') && (in(char, opening_chars) && in(previous_char, new_line_chars)) || in(previous_char, opening_chars)
        if !(char == ',') && (in(char, opening_chars) && previous_char == ',') || in(previous_char, opening_chars)
            write(buf, repeat(" ", indent))
        end
        if char in opening_chars
            write(buf, "$char\n")
            indent += 4
        elseif char in closing_chars
            previous_char != ',' && write(buf, "\n")
            indent -= 4
            write(buf, repeat(" ", indent))
            write(buf, char)
        elseif char == ','
            write(buf, "\n")
        else
            in(previous_char, new_line_chars) && write(buf, repeat(" ", indent))
            write(buf, char)
        end
        previous_char = char
    end
    String(take!(buf))
end

ext: jl


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Number of rounds", fontsize=14)
ylabel("Cross entropy", fontsize=14)
xticks([10, 100, 200, 500], fontsize=12)
yticks(1.46:0.005:1.475, fontsize=12)

savefig("assets/literate/EX-crabs-xgb-curve1.svg") # hide


text:
#=
 Misc routines that are associated with BED like files
=#

using DataFrames

"""
   load_dataframe_from_url

   Until I can figure out how to use Requests to download dropbox links
"""
function load_dataframe_from_url(url,readtable_options=Dict{Symbol,Char}(:separator=>'\t') )
    tmpfile = tempname()
    curl_cmd=`curl -L $url` |> tmpfile
    readall(curl_cmd)
    df= readtable(tmpfile;readtable_options...)
    return (df,tmpfile)
end

"""
  run_cmd

  To handle compatibility between 0.3 and 0.4
"""
function run_cmd(bedtools_cmd,bedtools_args,outfile)
  if Base.VERSION < v"0.4.0"
        cmd=`$bedtools_cmd $bedtools_args` |> outfile
        readall(cmd)
    else
       cmd=pipeline(`$bedtools_cmd $bedtools_args`,stdout=outfile)
       run(cmd)
    end
end


text:
module FromDigits

export fromdigits, fromdigits_unchecked

"""
Converts from digits to an integer. This is the inverse function of `digits`.

Returns the integer represented by the `digits` in the `base`.
The result is of type `T`, which defaults to `D`.

See also [`fromdigits_unchecked`](@ref).
"""
function fromdigits end

"""
    fromdigits(digits::AbstractVector{D}; base::B=D(10)) where {B<:Integer, D<:Integer}

Returns the integer represented by the `digits` in the `base`.

The result is of type `promote_type(B, D)`.
"""
fromdigits(digits::AbstractVector{D}; base::B=D(10)) where {B<:Integer, D<:Integer} =
    fromdigits(promote_type(B, D), digits; base=base)

"""
    fromdigits(T::Type{<:Integer}, digits::AbstractVector{D}; base::B=D(10)) where {B<:Integer, D<:Integer}

Returns the integer represented by the `digits` in the `base`. The result is of type `T`.
"""
function fromdigits(T::Type{<:Integer}, digits::AbstractVector{D}; base::B=D(10)) where {B<:Integer, D<:Integer}
    all(d -> d >= 0, digits) || all(d -> d <= 0, digits) || throw(ArgumentError("All the digits must have the same sign"))
    x = zero(T)
    for i in reverse(eachindex(digits))
        d = digits[i]
        abs(d) < base || throw(ArgumentError("Invalid digit $d for given base $base"))
        base_x = T(Base.Checked.checked_mul(base, x))
        x = T(Base.Checked.checked_add(base_x, d))
    end
    x
end

"""
Converts from digits to an integer. This is the inverse function of `digits`.

This function does not check the validity of the digits nor overflows in the computations.

See also [`fromdigits`](@ref).
"""
function fromdigits_unchecked end

fromdigits_unchecked(digits::AbstractVector{D}; base::B=D(10)) where {B<:Integer, D<:Integer} =
    fromdigits_unchecked(promote_type(B, D), digits; base=base)

function fromdigits_unchecked(T::Type{<:Integer}, digits::AbstractVector{D}; base::B=(10)) where {B<:Integer, D<:Integer}
    # return foldr((d, x) -> base*x + d, digits; init=zero(T))
    # return evalpoly(base, digits)
    x = zero(T)
    for i in reverse(eachindex(digits))
        x = base*x + digits[i]
    end
    x
end

end


text:
# Autogenerated wrapper script for CSFML_jll for i686-w64-mingw32-cxx03
export libcsfml_audio, libcsfml_graphics, libcsfml_network, libcsfml_system, libcsfml_window

using SFML_jll
JLLWrappers.@generate_wrapper_header("CSFML")
JLLWrappers.@declare_library_product(libcsfml_audio, "csfml-audio-2.dll")
JLLWrappers.@declare_library_product(libcsfml_graphics, "csfml-graphics-2.dll")
JLLWrappers.@declare_library_product(libcsfml_network, "csfml-network-2.dll")
JLLWrappers.@declare_library_product(libcsfml_system, "csfml-system-2.dll")
JLLWrappers.@declare_library_product(libcsfml_window, "csfml-window-2.dll")
function __init__()
    JLLWrappers.@generate_init_header(SFML_jll)
    JLLWrappers.@init_library_product(
        libcsfml_audio,
        "bin\\csfml-audio-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsfml_graphics,
        "bin\\csfml-graphics-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsfml_network,
        "bin\\csfml-network-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsfml_system,
        "bin\\csfml-system-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcsfml_window,
        "bin\\csfml-window-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
###########
# This gives an example of Gauss–Seidel
#
###########

using BlockBandedMatrices, BandedMatrices, ArrayLayouts, FillArrays, LazyBandedMatrices, LazyArrays

function finitedifference_2d(n)
    h = 1/n
    D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2
    D_xx = BandedBlockBandedMatrix(Kron(D², Eye(n)))
    D_yy = BandedBlockBandedMatrix(Kron(Eye(n), D²))
    D_xx + D_yy
end

function _gaussseidel(L, U, b, x=copy(b), y=copy(b), M=5)
    for _=1:M
        @view(y[1:end-1]) .= Mul(U , @view(x[2:end]))
        y[end] = 0
        x .= b .- y
        x .= Ldiv(L, x)
    end
    x
end

function _gaussseidel2(L, U, b, x=copy(b), y=copy(b), M=5)
    for _=1:M
        mul!(@view(y[1:end-1]) , U , @view(x[2:end]))
        y[end] = 0
        y .= b .- y
        x .= L\ y
    end
    x
end

function gaussseidel(A, b, M)
    n = Int(sqrt(length(b)))
    L = LowerTriangular(A)
    U = BandedBlockBandedMatrix(UpperTriangular(@view(A[1:end-1,2:end])), ([fill(n,n-1); n-1], [n-1; fill(n,n-1)]),
                                            (0,1), (0,1))
    x = copy(b)
    _gaussseidel(L,U, b, x, copy(x), M) # 1.6s
    x
end

n = 1000;
    Δt = (1/n^2)/4; @time Δ = finitedifference_2d(n);
    @time A = I - Δt*Δ  # 1m x 1m discretization
    b = randn(n^2)
    L = LowerTriangular(A)
    @time U = BandedBlockBandedMatrix(UpperTriangular(@view(A[1:end-1,2:end])), ([fill(n,n-1); n-1], [n-1; fill(n,n-1)]),
                                            (0,1), (0,1))
    x = copy(b)
    y = copy(x)
    @time _gaussseidel(L,U, b, x, y, 20) # 0.4s

    @time L̃ = sparse(BandedBlockBandedMatrix(L))
    @time Ũ = sparse(U)
    x = copy(b)
    y = copy(x)
    @time _gaussseidel2(L̃,Ũ, b, x, y, 20)

Ã = sparse(A)
@time qr(Ã);
@time Ã \b ;

A = randn(n,n)

n = 1000;
    As = [randn(2k+1,k) for k=1:n];

    @time for A in As
        qr(A)
    end
n = 1000;
    # As = [randn(2k+1,k) for k=1:n];

    @time for A in As
        lu(A)
    end


text:
export NodePosition, Node, InputNode, InteriorNode, OutputNode

using Compat

typealias NodePosition @compat Tuple{Int, Int}

abstract Node

type InputNode <: Node
    index::Integer
    active::Bool
end

InputNode(index::Integer) = InputNode(index, false)

type InteriorNode <: Node
    func::Func
    inputs::Vector{NodePosition}
    active::Bool
    cache::BitString
end

InteriorNode(func::Func, inputs::Vector{NodePosition}) = InteriorNode(func, inputs, false, convert(BitString, 0))

type OutputNode <: Node
    input::NodePosition
    active::Bool
    cache::BitString
end

OutputNode(input::NodePosition) = OutputNode(input, false, convert(BitString, 0))


text:
export BackwardPlanner, BackwardGreedyPlanner, BackwardAStarPlanner

"Heuristic-guided best-first backward search."
@kwdef mutable struct BackwardPlanner <: Planner
    heuristic::Heuristic = GoalCountHeuristic(:backward)
    g_mult::Float64 = 1.0 # Path cost multiplier
    h_mult::Float64 = 1.0 # Heuristic multiplier
    max_nodes::Int = typemax(Int)
    save_search::Bool = false # Flag to save search info
end

"Backward greedy search, with cycle checking."
BackwardGreedyPlanner(heuristic::Heuristic; kwargs...) =
    BackwardPlanner(;heuristic=heuristic, g_mult=0, kwargs...)

"Backward A* search."
BackwardAStarPlanner(heuristic::Heuristic; kwargs...) =
    BackwardPlanner(;heuristic=heuristic, kwargs...)

function solve(planner::BackwardPlanner,
               domain::Domain, state::State, spec::Specification)
    @unpack h_mult, heuristic, save_search = planner
    # Precompute heuristic information
    precompute!(heuristic, domain, state, spec)
    # Convert to backward search goal specification
    spec = BackwardSearchGoal(spec, state)
    state = State(get_goal_terms(spec), PDDL.get_types(state))
    # Initialize search tree and priority queue
    node_id = hash(state)
    search_tree = Dict{UInt,PathNode}(node_id => PathNode(node_id, state, 0))
    est_cost = h_mult * heuristic(domain, state, spec)
    queue = PriorityQueue{UInt,Float64}(node_id => est_cost)
    # Run the search
    status, node_id = search!(planner, domain, spec, search_tree, queue)
    # Reconstruct plan and return solution
    if status != :failure
        plan, traj = reconstruct(node_id, search_tree)
        reverse!(plan); reverse!(traj)
        if save_search
            return PathSearchSolution(status, plan, traj, search_tree, queue)
        else
            return PathSearchSolution(status, plan, traj)
        end
    elseif save_search
        return PathSearchSolution(status, [], [], search_tree, queue)
    else
        return NullSolution()
    end
end

function search!(planner::BackwardPlanner,
                 domain::Domain, spec::BackwardSearchGoal,
                 search_tree::Dict{UInt,PathNode}, queue::PriorityQueue)
    count = 1
    while length(queue) > 0
        # Get state with lowest estimated cost to start state
        node_id = dequeue!(queue)
        node = search_tree[node_id]
        # Return status and current state if search terminates
        if is_goal(spec, domain, node.state)
            return :success, node_id # Start state reached
        elseif count >= planner.max_nodes
            return :max_nodes, node_id # Node budget reached
        end
        count += 1
        # Expand current node
        expand!(planner, node, search_tree, queue, domain, spec)
    end
    return :failure, nothing
end

function expand!(planner::BackwardPlanner, node::PathNode,
                 search_tree::Dict{UInt,PathNode}, queue::PriorityQueue,
                 domain::Domain, spec::BackwardSearchGoal)
    @unpack g_mult, h_mult, heuristic = planner
    state = node.state
    # Iterate over relevant actions
    actions = relevant(state, domain)
    for act in actions
        # Regress (reverse-execute) the action
        next_state = regress(act, state, domain; check=false)
        # Add constraints to regression state
        add_constraints!(spec, state)
        next_id = hash(next_state)
        # Compute path cost
        act_cost = get_cost(spec, domain, state, act, next_state)
        path_cost = node.path_cost + act_cost
        # Update path costs if new path is shorter
        next_node = get!(search_tree, next_id,
                         PathNode(next_id, next_state, Inf))
        cost_diff = next_node.path_cost - path_cost
        if cost_diff > 0
            next_node.parent_id = node.id
            next_node.parent_action = act
            next_node.path_cost = path_cost
            # Update estimated cost from next state to start
            if !(next_id in keys(queue))
                g_val = g_mult * path_cost
                h_val = h_mult * heuristic(domain, next_state, spec)
                enqueue!(queue, next_id, g_val + h_val)
            else
                queue[next_id] -= cost_diff
            end
        end
    end
end


text:
using Plots

"Create a 2D array of cells with 1 infected cell in the middle"
function make_cells(width::Integer=7, height::Integer=7)
    cells = Matrix{Int8}(undef, width, height)
    for i in 1:size(cells)[1]
        for j in 1:size(cells)[2]
            cells[i,j] = 0
        end
    end
    cells[width÷2+1,height÷2+1] = 1
    return cells
end

"Map the cells to colors for plotting"
function to_colors(cell)
    if cell == 0
        return RGB(0.0,0.8,0.0)
    end
    if cell == 1
        return RGB(0.8,0.0,0.0)
    end
end

"""
Simulate interaction between 2 neighbouring cells. If one is infected, it can infect
the other.
"""
function interact(cell, neighbour, infection_rate=0.1)
    new_cell = cell

    # Check if the cell is infected
    if neighbour == 1

        # Infect with probability 0.1
        if rand(1)[1] < infection_rate
            new_cell = 1
        end
    end

    return new_cell
end

"Update the simulation one time step"
function update!(cells, infection_rate::Real=0.1)
    # Create a copy to hold the old state
    old_cells = deepcopy(cells)

    # Loop over pairs of cells in the same row. There are size(cells)[1] columns, and size(cells)[1]-1 pairs.
    for i in 1:size(cells)[1]-1
        # loop over all columns
        for j in 1:size(cells)[2]
            # So the cells are (i+1,j) and (i,j)
            cells[i,j] = interact(cells[i,j], old_cells[i+1,j])
            cells[i+1,j] = interact(cells[i+1,j], old_cells[i,j])
        end
    end

    # Then loop over pairs in the same columns.
    for i in 1:size(cells)[1]
        # loop over all columns
        for j in 1:size(cells)[2]-1
            # So the cells are (i,j+1) and (i,j)
            cells[i,j] = interact(cells[i,j], old_cells[i,j+1])
            cells[i,j+1] = interact(cells[i,j+1], old_cells[i,j])
        end
    end
end

"Count the current number of infections"
function count_infections(cells::Matrix{Int8})
    infections = 0
    for i in 1:size(cells)[1]
        for j in 1:size(cells)[2]
            if cells[i,j].status == infected
                infections += 1
            end
        end
    end
    return infections
end

max_stars_repo_name: tpfau/julia-introduction

max_issues_repo_name: tpfau/julia-introduction

max_forks_repo_name: tpfau/julia-introduction


text:
module SimpleDifferentialOperators

# Pkg Dependencies
using LinearAlgebra, SparseArrays, Parameters

# Includes
include("types.jl")
include("operators.jl")
include("extensionoperators.jl")
include("utilities/extrapolatetoboundary.jl")

# Exports
# Boundary Conditions
export BoundaryCondition,
       HomogeneousBoundaryCondition,
       Reflecting,
       Mixed,
       Absorbing

# Differential Types
export DifferenceMethod,
        ForwardFirstDifference,
        BackwardFirstDifference,
        CentralSecondDifference

# Functions
export DifferentialOperator,
        ExtensionDifferentialOperator,
        L₁₋bc,
        L₁₊bc,
        L₂bc,
        interiornodes,
        L₁₊,
        L₁₋,
        L₂,
        extrapolatetoboundary

end # module


text:
"""
    BinPackingCapacityType

Whether the capacities of a `BinPacking` constraint are fixed:

* either there is no capacity: `NO_CAPACITY_BINPACKING`
* or the capacity values are fixed when creating the set: `FIXED_CAPACITY_BINPACKING`
* or the capacity values are themselves variable: `VARIABLE_CAPACITY_BINPACKING`
"""
@enum BinPackingCapacityType begin
    NO_CAPACITY_BINPACKING
    FIXED_CAPACITY_BINPACKING
    VARIABLE_CAPACITY_BINPACKING
end

"""
    BinPacking(n_bins::Int, n_items::Int, weights::Vector{T})

## Uncapacitated bin packing

Implements an uncapacitated version of the bin-packing problem.

The first `n_bins` variables give the load in each bin, the last `n_items` give
the number of the bin to which the item is assigned to.

The load of a bin is defined as the sum of the sizes of the items put in that
bin.

Also called [`pack`](https://sofdem.github.io/gccat/gccat/Cbin_packing.html).

### Example

    [a, b, c] in BinPacking{NO_CAPACITY_BINPACKING}(1, 2, [2, 3])
    # As there is only one bin, the only solution is to put all the items in
    # that bin.
    # Enforces that:
    # - the bin load is the sum of the weights of the objects in that bin:
    #   a = 2 + 3
    # - the bin number of the two items is 1: b = c = 1

## Fixed-capacity bin packing

Implements a capacitated version of the bin-packing problem where capacities
are constant.

The first `n_bins` variables give the load in each bin, the last `n_items`
give the number of the bin to which the item is assigned to.

The load of a bin is defined as the sum of the sizes of the items put in that
bin.

This constraint is equivalent to `BinPacking` with inequality constraints on
the loads of the bins where the upper bound is a constant. However, there are
more efficient propagators for the combined constraint (bin packing with
maximum load); if such propagators are not available, bridges are available to
make the conversion seamless.

Also called [`bin_packing_capa`](https://sofdem.github.io/gccat/gccat/Cbin_packing_capa.html).

### Example
    [a, b, c] in BinPacking{FIXED_CAPACITY_BINPACKING}(1, 2, [2, 3], [4])
    # As there is only one bin, the only solution is to put all the items in
    # that bin if its capacity is large enough.
    # Enforces that:
    # - the bin load is the sum of the weights of the objects in that bin:
    #   a = 2 + 3
    # - the bin load is at most its capacity: a <= 4 (given in the set)
    # - the bin number of the two items is 1: b = c = 1

## Variable-capacity bin packing

Implements an capacitated version of the bin-packing problem where capacities
are optimisation variables.

The first `n_bins` variables give the load in each bin, the next `n_bins` are
the capacity of each bin, the last `n_items` give the number of the bin to
which the item is assigned to.

The load of a bin is defined as the sum of the sizes of the items put in that
bin.

This constraint is equivalent to `BinPacking` with inequality constraints on
the loads of the bins where the upper bound is any expression. However, there
are more efficient propagators for the combined constraint (bin packing with
maximum load) and for the fixed-capacity version.

Also called [`bin_packing_capa`](https://sofdem.github.io/gccat/gccat/Cbin_packing_capa.html).

## Example
    [a, 2, b, c] in BinPacking{VARIABLE_CAPACITY_BINPACKING}(1, 2, [2, 3])
    # As there is only one bin, the only solution is to put all the items in
    # that bin if its capacity is large enough.
    # Enforces that:
    # - the bin load is the sum of the weights of the objects in that bin:
    #   a = 2 + 3
    # - the bin load is at most its capacity: a <= 2 (given in a variable)
    # - the bin number of the two items is 1: b = c = 1
"""
struct BinPacking{BPCT, T <: Real} <: MOI.AbstractVectorSet
    n_bins::Int
    n_items::Int
    weights::Vector{T}
    capacities::Vector{T}

    function BinPacking{BPCT}(
        n_bins::Int,
        n_items::Int,
        weights::Vector{T},
        capacities::Vector{T},
    ) where {BPCT, T <: Real}
        @assert n_items == length(weights)
        @assert all(weights .>= zero(T))

        if BPCT == FIXED_CAPACITY_BINPACKING
            @assert n_bins == length(capacities)
            @assert all(capacities .>= zero(T))
        end

        @assert n_bins > 0
        @assert n_items > 0

        return new{BPCT, T}(n_bins, n_items, weights, capacities)
    end
end

function BinPacking{NO_CAPACITY_BINPACKING}(n_bins::Int, n_items::Int, weights::Vector{T}) where {T <: Real}
    return BinPacking{NO_CAPACITY_BINPACKING}(n_bins, n_items, weights, T[])
end

function BinPacking(n_bins::Int, n_items::Int, weights::Vector{T}, capacities::Vector{T}) where {T <: Real}
    return BinPacking{FIXED_CAPACITY_BINPACKING}(n_bins, n_items, weights, capacities)
end

function BinPacking{VARIABLE_CAPACITY_BINPACKING}(n_bins::Int, n_items::Int, weights::Vector{T}) where {T <: Real}
    return BinPacking{VARIABLE_CAPACITY_BINPACKING}(n_bins, n_items, weights, T[])
end

MOI.dimension(set::BinPacking{NO_CAPACITY_BINPACKING, T}) where {T} = set.n_bins + set.n_items
MOI.dimension(set::BinPacking{FIXED_CAPACITY_BINPACKING, T}) where {T} = set.n_bins + set.n_items
MOI.dimension(set::BinPacking{VARIABLE_CAPACITY_BINPACKING, T}) where {T} = 2 * set.n_bins + set.n_items

function copy(set::BinPacking{BPCT, T}) where {BPCT, T}
    return BinPacking{BPCT}(
        set.n_bins,
        set.n_items,
        copy(set.weights),
        copy(set.capacities),
    )
end

function Base.:(==)(x::BinPacking{BPCT, T}, y::BinPacking{BPCT, T}) where {BPCT, T}
    return x.n_bins == y.n_bins &&
           x.n_items == y.n_items &&
           x.weights == y.weights &&
           x.capacities == y.capacities
end

# Shortcuts for types.
# const BinPacking{T} = BinPacking{NO_CAPACITY_BINPACKING, T}
const BinPackingNoCapacity{T} = BinPacking{NO_CAPACITY_BINPACKING, T}
const BinPackingFixedCapacity{T} = BinPacking{FIXED_CAPACITY_BINPACKING, T}
const BinPackingVariableCapacity{T} = BinPacking{VARIABLE_CAPACITY_BINPACKING, T}


text:
"""
    simple_normalized_residuals(data::Dict, se_sol::Dict, R::Matrix)
It normalizes the residuals only based on their standard deviation, no sensitivity matrix involved.
R^2 replaces Ω, where Ω = S ⋅ R
Avoids all the matrix calculations but is a "simplified" method
"""
function simple_normalized_residuals(data::Dict, se_sol::Dict, R::Matrix)
    normalized_residuals(data, se_sol, R.^2)
end
"""
Adds the normalized residuals to the solution dictionary and returns the largest normalized residual,
its index (i.e., the measurement it refers to), and whether it exceeds the given threshold `t` or not.
"""
function normalized_residuals(data::Dict, se_sol::Dict, Ω::Matrix; t::Float64=3.0)
    lnr = ("0", 0)
    count = 1
    for (m, meas) in data["meas"]
        if !haskey(se_sol["solution"][string(meas["cmp"])], string(meas["cmp_id"]))
            r = [0.0, 0.0, 0.0]
            se_sol["solution"]["meas"][m]["r"] = [0.0, 0.0, 0.0]
        else
            h_x = se_sol["solution"][string(meas["cmp"])][string(meas["cmp_id"])][string(meas["var"])]
            z = _DST.mean.(meas["dst"])
            r = abs.(h_x-z)
            se_sol["solution"]["meas"][m]["r"] = r
        end
        se_sol["solution"]["meas"][m]["nr"] = [r[i]/sqrt(abs(Ω[count+i-1, count+i-1])) for i in 1:length(meas["dst"])]
        for i in 1:length(meas["dst"])
            if  se_sol["solution"]["meas"][m]["nr"][i] > last(lnr) lnr = (m,  se_sol["solution"]["meas"][m]["nr"][i]) end
        end
        count+=length(meas["dst"])
    end
    return lnr, last(lnr) > t
end
"""
Returns the Measurement Jacobian H
"""
function build_H_matrix(functions::Vector, state::Array)::Matrix{Float64}
    H = Matrix{Float64}(undef, length(functions), length(state))
    for row in 1:length(functions)
        H[row,:] = ForwardDiff.gradient(functions[row], state)
    end
    return H
end
"""
Returns the Gain Matrix G
"""
build_G_matrix(H::Matrix, R::Matrix)::Matrix{Float64} = transpose(H)*inv(R)*H
"""
Returns the Measurement Error Covariance Matrix R
"""
function build_R_matrix(data::Dict)::Matrix{Float64}
    meas_row_order = [m for (m, _) in data["meas"]]
    R_entries = vcat([_DST.std.(data["meas"][mid]["dst"])[1:length(data["meas"][mid]["dst"])] for mid in meas_row_order]...)
    return LinearAlgebra.diagm(R_entries.^2)
end
"""
Returns the Residual Covariance Matrix Ω, where:
Ω = R - H*G^(-1)*H^T
"""
build_omega_matrix(R::Matrix{Float64}, H::Matrix{Float64}, G::Matrix{Float64}) = R - H*inv(G)*transpose(H)
"""
Returns the Residual Covariance Matrix Ω, where:
Ωᵢᵢ = Rᵢᵢ⋅Sᵢᵢ
"""
build_omega_matrix(S::Matrix{Float64}, R::Matrix{Float64}) = diagm((sqrt.(abs.(diag(S).*diag(R)))))
"""
Returns the sensitivity matrix S starting from the hat matrix K
"""
build_S_matrix(K::Matrix{Float64}) = Matrix{Int64}(I, size(K)) - K
"""
Returns the hat matrix K
"""
build_K_matrix(H::Matrix{Float64}, G::Matrix{Float64}, R::Matrix{Float64}) = H*inv(G)*transpose(H)*inv(R)


text:
# Exercise 18-1
struct MyTime
    hour :: Int64
    minute :: Int64
    second :: Int64
end

function timetoint(time::MyTime)
    minutes = time.hour * 60 + time.minute
    seconds = minutes * 60 + time.second
end

function inttotime(seconds::Int64)
    (minutes, second) = divrem(seconds, 60)
    hour, minute = divrem(minutes, 60)
    MyTime(hour, minute, second)
end

import Base.<

function <(t1::MyTime, t2::MyTime)
   return timetoint(t1) < timetoint(t2)
end

# Alternatively
# function <(t1::MyTime, t2::MyTime)
#     return (t1.hour, t1.minute, t1.second) < (t2.hour, t2.minute, t2.second)
# end

# Exercise 18-2
const suit_names = ["♣", "♦", "♥", "♠"]
const rank_names = ["A", "2", "3", "4", "5", "6", "7", "8", "9", "10", "J", "Q", "K"]

struct Card
    suit :: Int64
    rank :: Int64
    function Card(suit::Int64, rank::Int64)
        @assert(1 ≤ suit ≤ 4, "suit is not between 1 and 4")
        @assert(1 ≤ rank ≤ 13, "rank is not between 1 and 13")
        new(suit, rank)
    end
end

function Base.show(io::IO, card::Card)
    print(io, rank_names[card.rank], suit_names[card.suit])
end

function Base.isless(c1::Card, c2::Card)
    (c1.suit, c1.rank) < (c2.suit, c2.rank)
end

struct Deck
    cards :: Array{Card, 1}
end

function Deck()
    deck = Deck(Card[])
    for suit in 1:4
        for rank in 1:13
            push!(deck.cards, Card(suit, rank))
        end
    end
    deck
end

using Random

function Random.shuffle!(deck::Deck)
    shuffle!(deck.cards)
    deck
end

function Base.sort!(deck::Deck)
    sort!(deck.cards)
    deck
end

function ex18_2()
    d = Deck()
    shuffle!(d)
    println(d)
    sort!(d)
    println(d)
end

# ex18_2()

# Exercise 18-3
struct Markov
    order :: Int64
    suffixes :: Dict{Tuple{String,Vararg{String}}, Array{String, 1}}
    prefix :: Array{String, 1}
    function Markov(order::Int64=2)
        new(order, Dict{Tuple{String,Vararg{String}}, Array{String, 1}}(), Array{String, 1}())
    end
end

function processword(markov::Markov, word::String)
    if length(markov.prefix) < markov.order
        push!(markov.prefix, word)
        return
    end
    get!(markov.suffixes, (markov.prefix...,), Array{String, 1}())
    push!(markov.suffixes[(markov.prefix...,)], word)
    popfirst!(markov.prefix)
    push!(markov.prefix, word)
end

function ex18_3()
    m = Markov()
    for x in ["the", "man", "was", "as", "the", "man", "saw", ","]
        processword(m, x)
    end
    println(m)
end

# ex18_3()

# Exercise 18-5
function deal!(deck::Deck, hands::Int, size::Int)
    arr = Array{Hand, 1}()
    for j = 1:hands
        push!(arr, Hand())
    end

    for j = 1:hands
        move!(deck, arr[j], size)
    end
    return arr
end

deck = Deck()
# println(deal!(deck, 3, 5))


text:
"Solve Computing Eccentric Anomaly from Mean Anomally via Kepler's Equation"
module KeplerEqn

export calc_ecc_anom
export calc_ecc_anom_cpu, calc_ecc_anom_gpu
using CUDAnative

"""   ecc_anom_init_guess_danby(M, ecc)
Initial guess for eccentric anomaly given mean anomaly (M) and eccentricity (ecc)
    Based on "The Solution of Kepler's Equations - Part Three"
    Danby, J. M. A. (1987) Journal: Celestial Mechanics, Volume 40, Issue 3-4, pp. 303-312  1987CeMec..40..303D
"""
function ecc_anom_init_guess_danby(M, ecc)
    k = convert(typeof(ecc),0.85)
    if(M<zero(M)) M += 2pi end
    (M<pi) ? M + k*ecc : M - k*ecc;
end

"""   update_ecc_anom_laguerre(E, M, ecc)
Update the current guess (E) for the solution to Kepler's equation given mean anomaly (M) and eccentricity (ecc)
   Based on "An Improved Algorithm due to Laguerre for the Solution of Kepler's Equation"
   Conway, B. A.  (1986) Celestial Mechanics, Volume 39, Issue 2, pp.199-211  1986CeMec..39..199C
"""
function update_ecc_anom_laguerre(E, M, ecc)
  es = ecc*sin(E)
  ec = ecc*cos(E)
  F = (E-es)-M
  Fp = one(E)-ec
  Fpp = es
  n = 5
  root = sqrt(abs((n-1)*((n-1)*Fp*Fp-n*F*Fpp)))
  denom = Fp>zero(E) ? Fp+root : Fp-root
  return E-n*F/denom
end

"Same as update_ecc_anom_laguerre_cpu, except uses GPU math functions"
function update_ecc_anom_laguerre_gpu(E, M, ecc)
  es = ecc*CUDAnative.sin(E)
  ec = ecc*CUDAnative.cos(E)
  F = (E-es)-M
  Fp = one(E)-ec
  Fpp = es
  n = 5
  root = CUDAnative.sqrt(abs((n-1)*((n-1)*Fp*Fp-n*F*Fpp)))
  denom = Fp>zero(E) ? Fp+root : Fp-root
  return E-n*F/denom
end

"Loop to update the current estimate of the solution to Kepler's equation"
function calc_ecc_anom_itterative_laguerre(mean_anom, ecc, tol, max_its)
    M = mod(mean_anom,convert(typeof(mean_anom),2pi))
    E = ecc_anom_init_guess_danby(M,ecc)
    for i in 1:max_its
       E_old = E
       E = update_ecc_anom_laguerre(E_old, M, ecc)
       if abs(E-E_old)<convert(typeof(mean_anom),tol) break end
    end
    return E
end

"Loop to update the current estimate of the solution to Kepler's equation.  Calls GPU update"
function calc_ecc_anom_itterative_laguerre_gpu(mean_anom, ecc, tol, max_its)
    @assert zero(ecc) <= ecc < one(ecc)
    @assert tol*100 <= one(tol)
    M = mod(mean_anom,convert(typeof(mean_anom),2pi))
    E = ecc_anom_init_guess_danby(M,ecc)
    for i in 1:max_its
       E_old = E
       E = update_ecc_anom_laguerre_gpu(E_old, M, ecc)
       if abs(E-E_old)<convert(typeof(mean_anom),tol) break end
    end
    return E
end

const default_max_its_laguerre = 200
const default_ecc_anom_tol = 1e-8
"Calculate eccentric anomaly given mean anomaly and eccentricty (in radians)"
calc_ecc_anom_cpu(mean_anom, ecc) = calc_ecc_anom_itterative_laguerre(mean_anom, ecc, default_ecc_anom_tol, default_max_its_laguerre)
calc_ecc_anom_gpu(mean_anom, ecc) = calc_ecc_anom_itterative_laguerre_gpu(mean_anom, ecc, default_ecc_anom_tol, default_max_its_laguerre)
calc_ecc_anom(mean_anom, ecc) = calc_ecc_anom_cpu(mean_anom, ecc)

end # module KeplerEqn


text:
using IntegralTransforms
abstract type RealFunction <: Function end
abstract type ComplexFunction <: Function end
abstract type AbstractIntegralTransform end
abstract type AbstractTransform <: AbstractIntegralTransform end
abstract type AbstractLaplaceTransform <: AbstractIntegralTransform end
abstract type AbstractForwardLaplaceTransform <: AbstractLaplaceTransform end
abstract type AbstractInverseLaplaceTransform <: AbstractLaplaceTransform end

# struct RealF{Func <: Function, T <: Real} <: Real
#     F :: Func
#     t :: T
# end

struct CompF{Func <: Function, S <: Complex}
    f :: Func
    s :: S
end

"""
     ℒ(f)

Laplace Transform of a real function. Its inverse is the invℒ(f) operator.
- Markdown works in here.
- ``LaTeX``
```math
SameAs LaTeX
```
"""
ℒ(x) = 1

# """
#     +()
# We can add Transformed functions and still get something cool.
# """
# @inline Base.:+(a::RealF, b::RealF) = a.F(a.t) + b.F(b.t)

""" func(x) is 1 """
func(x) = 1

"yay"
func2(x) = 1

"yay"
func3(x) = 1

"yay"
func4(x) = 1

"yay"
func5(x) = 1

"yay"
func6(x) = 1

# @inline Base.:+(a :: RealF, b:: RealF) = a.F(a.t) + b.F(b.t)
# @inline Base.:*(a :: RealF, α :: Real) = α * a.F(a.t)
# @inline Base.:*(a :: Real, α :: RealF) = α * a.F(a.t)
#
# y = RealF(*,5)
#
# x = RealF(*, 6)
#
# @inline Base.:+(a :: ℒ, g :: ℒ) = a + b
# @inline Base:+(a :: ℒ,
#
#
# f = 1/s


text:
struct EventRates{T <: DiseaseStateSequence}
  exposure::Union{Nothing, Vector{Float64}}
  infection::Vector{Float64}
  removal::Union{Nothing, Vector{Float64}}

  function EventRates{T}(n::Int64) where T <: SEIR
    return new{T}(fill(0.0, n), fill(0.0, n), fill(0.0, n))
  end

  function EventRates{T}(n::Int64) where T <: SEI
    return new{T}(fill(0.0, n), fill(0.0, n), nothing)
  end

  function EventRates{T}(n::Int64) where T <: SIR
    return new{T}(nothing, fill(0.0, n), fill(0.0, n))
  end

  function EventRates{T}(n::Int64) where T <: SI
    return new{T}(nothing, fill(0.0, n), nothing)
  end
end

function individuals(x::EventRates{T}) where {
                     T <: DiseaseStateSequence}
  return length(x.infection)
end

function Base.getindex(x::EventRates{T}, new_state::DiseaseState) where T <: DiseaseStateSequence
  if new_state == State_E
    return x.exposure
  elseif new_state == State_I
    return x.infection
  elseif new_state == State_R
    return x.removal
  else
    @error "Unrecognized indexing disease state"
  end
end

function Base.getindex(x::EventRates{T}, new_states::DiseaseStates) where T <: DiseaseStateSequence
  y = x[new_states[1]]
  for i=2:length(new_states)
    y = hcat(y, x[new_states[i]])
  end
  return y
end

function Base.sum(x::EventRates{T}) where T <: DiseaseStateSequence
  return sum([sum(x[i]) for i in convert(DiseaseStates, T)[2:end]])
end


text:
"A macro to make calling KNITRO's KN_* C API a little cleaner"
macro kn_ccall(func, args...)
    f = Base.Meta.quot(Symbol("KN_$(func)"))
    args = [esc(a) for a in args]
    quote
        ccall(($f, libknitro), $(args...))
    end
end

macro kn_get_values(function_name, type)
    name_singular = "KN_" * string(function_name)
    name_plural = "KN_" * string(function_name) * "s"
    fname = Symbol(name_plural)
    fnameshort = Symbol(name_singular)
    # Names of C function in knitro.h
    c_fname = Symbol(name_singular)
    c_fnames = Symbol(name_plural)
    c_fnames_all = Symbol(name_plural * "_all")

    n = if occursin("var", name_singular)
        :(KN_get_number_vars(kc))
    elseif occursin("con", name_singular)
        :(KN_get_number_cons(kc))
    end

    quote
        function $(esc(fname))(kc::Model)
            result = zeros($type, $n)
            $c_fnames_all(kc, result)
            return result
        end
        function $(esc(fname))(kc::Model, index::Vector{Cint})
            result = zeros($type, length(index))
            $c_fnames(kc, length(index), index, result)
            return result
        end
        function $(esc(fname))(kc::Model, index::Integer)
            result = zeros($type, 1)
            $c_fname(kc, index, result)
            return result[1]
        end
        $(esc(fnameshort))(kc::Model, index::Integer) = $(esc(fname))(kc, index)
    end
end

macro kn_get_attribute(function_name, type)
    fname = Symbol("KN_" * string(function_name))
    quote
        function $(esc(fname))(m::Model)
            val = zeros($type, 1)
            ret = $fname(m, val)
            return val[1]
        end
    end
end

"Format output returned by KNITRO as proper Julia string."
function _format_output(output::AbstractString)
    # remove trailing whitespace
    res = strip(output)
    # remove special characters
    return res = strip(res, '\0')
end

"Return the current KNITRO version."
function get_release()
    len = 15
    out = zeros(Cchar, len)
    KN_get_release(len, out)
    return String(strip(String(convert(Vector{UInt8}, out)), '\0'))
end

"Wrapper for KNITRO KN_context."
mutable struct Env
    ptr_env::Ptr{Cvoid}

    function Env()
        ptrptr_env = Ref{Ptr{Cvoid}}()
        res = KN_new(ptrptr_env)
        if res != 0
            error("Fail to retrieve a valid KNITRO KN_context. Error $res")
        end
        return new(ptrptr_env[])
    end

    Env(ptr::Ptr{Cvoid}) = new(ptr)
end

Base.unsafe_convert(ptr::Type{Ptr{Cvoid}}, env::Env) = env.ptr_env::Ptr{Cvoid}
is_valid(env::Env) = env.ptr_env != C_NULL

"""
Free all memory and release any Knitro license acquired by calling KN_new.
"""
function free_env(env::Env)
    if env.ptr_env != C_NULL
        ptrptr_env = Ref{Ptr{Cvoid}}(env.ptr_env)
        KN_free(ptrptr_env)
        env.ptr_env = C_NULL
    end
    return
end

"""
Structure specifying the callback context.

Each evaluation callbacks (for objective, gradient or hessian)
is attached to a unique callback context.

"""
mutable struct CallbackContext
    context::Ptr{Cvoid}
    n::Int
    m::Int
    # Add a dictionnary to store user params.
    userparams::Any

    # Oracle's callbacks are context dependent, so store
    # them inside dedicated CallbackContext.
    eval_f::Function
    eval_g::Function
    eval_h::Function
    eval_rsd::Function
    eval_jac_rsd::Function

    function CallbackContext(ptr_cb::Ptr{Cvoid})
        return new(ptr_cb, 0, 0, nothing)
    end
end

Base.unsafe_convert(ptr::Type{Ptr{Cvoid}}, cb::CallbackContext) = cb.context::Ptr{Cvoid}

mutable struct Model
    # KNITRO context environment.
    env::Env
    # Keep reference to callbacks for garbage collector.
    callbacks::Vector{CallbackContext}
    # Some structures for userParams
    puts_user::Any
    multistart_user::Any
    mip_user::Any
    newpoint_user::Any

    # Solution values.
    # Optimization status. Equal to 1 if problem is unsolved.
    status::Cint
    obj_val::Cdouble
    x::Vector{Cdouble}
    mult::Vector{Cdouble}

    # Special callbacks (undefined by default).
    # (this functions do not depend on callback environments)
    ms_process::Function
    newpt_callback::Function
    mip_callback::Function
    user_callback::Function
    ms_initpt_callback::Function
    puts_callback::Function

    # Constructor.
    function Model()
        model = new(
            Env(),
            CallbackContext[],
            nothing,
            nothing,
            nothing,
            nothing,
            1,
            Inf,
            Cdouble[],
            Cdouble[],
        )
        # Add a destructor to properly delete model.
        finalizer(KN_free, model)
        return model
    end
    # Instantiate a new Knitro instance in current environment `env`.
    function Model(env::Env)
        return new(
            env,
            CallbackContext[],
            nothing,
            nothing,
            nothing,
            nothing,
            1,
            Inf,
            Cdouble[],
            Cdouble[],
        )
    end
end

Base.unsafe_convert(ptr::Type{Ptr{Cvoid}}, kn::Model) = kn.env.ptr_env::Ptr{Cvoid}

"Free solver object."
KN_free(m::Model) = free_env(m.env)

"Create solver object."
KN_new() = Model()

is_valid(m::Model) = is_valid(m.env)
has_callbacks(m::Model) = !isempty(m.callbacks)
register_callback(model::Model, cb::CallbackContext) = push!(model.callbacks, cb)

function Base.show(io::IO, m::Model)
    if is_valid(m)
        println(io, "$(get_release())")
        println(io, "-----------------------")
        println(io, "Problem Characteristics")
        println(io, "-----------------------")
        println(io, "Objective goal:  Minimize")
        println(io, "Objective type:  $(KN_get_obj_type(m))")
        println(
            io,
            "Number of variables:                             $(KN_get_number_vars(m))",
        )
        println(
            io,
            "Number of constraints:                           $(KN_get_number_cons(m))",
        )
        println(
            io,
            "Number of nonzeros in Jacobian:                  $(KN_get_jacobian_nnz(m))",
        )
        println(
            io,
            "Number of nonzeros in Hessian:                   $(KN_get_hessian_nnz(m))",
        )

    else
        println(io, "KNITRO Problem: NULL")
    end
end

#=
    LM license manager
=#
"""
Type declaration for the Artelys License Manager context object.
Applications must not modify any part of the context.
"""
mutable struct LMcontext
    ptr_lmcontext::Ptr{Cvoid}
    # Keep a pointer to instantiated models in order to free
    # memory properly.
    linked_models::Vector{Model}

    function LMcontext()
        ptrref = Ref{Ptr{Cvoid}}()
        res = KN_checkout_license(ptrref)
        if res != 0
            error("KNITRO: Error checkout license")
        end
        lm = new(ptrref[], Model[])
        finalizer(KN_release_license, lm)
        return lm
    end
end
Base.unsafe_convert(ptr::Type{Ptr{Cvoid}}, lm::LMcontext) = lm.ptr_lmcontext::Ptr{Cvoid}

function Env(lm::LMcontext)
    ptrptr_env = Ref{Ptr{Cvoid}}()
    res = KN_new_lm(lm, ptrptr_env)
    if res != 0
        error("Fail to retrieve a valid KNITRO KN_context. Error $res")
    end
    return Env(ptrptr_env[])
end

function attach!(lm::LMcontext, model::Model)
    push!(lm.linked_models, model)
    return
end

# create Model with license manager
function Model(lm::LMcontext)
    model = Model(Env(lm))
    attach!(lm, model)
    return model
end
KN_new_lm(lm::LMcontext) = Model(lm)

function KN_release_license(lm::LMcontext)
    # First, ensure that all linked models are properly freed
    # before releasing license manager!
    KN_free.(lm.linked_models)

    if lm.ptr_lmcontext != C_NULL
        refptr = Ref{Ptr{Cvoid}}(lm.ptr_lmcontext)
        KN_release_license(refptr)
        lm.ptr_lmcontext = C_NULL
    end
end

#=
    VARIABLES
=#
function KN_add_var(m::Model)
    index = Cint[0]
    KN_add_var(m, index)
    return index[1]
end
function KN_add_vars(m::Model, nvars::Int)
    indexes = zeros(Cint, nvars)
    KN_add_vars(m, nvars, indexes)
    return indexes
end

@kn_get_values get_var_lobnd Cdouble
@kn_get_values get_var_upbnd Cdouble
@kn_get_values get_var_eqbnd Cdouble
@kn_get_values get_var_fxbnd Cdouble
@kn_get_values get_var_primal_value Cdouble
@kn_get_values get_var_dual_value Cdouble
@kn_get_values get_var_type Cint

#=
    OBJECTIVE
=#
function KN_add_obj_linear_struct(
    m::Model,
    objIndices::Vector{Cint},
    objCoefs::Vector{Cdouble},
)
    nnz = length(objIndices)
    @assert nnz == length(objCoefs)
    return KN_add_obj_linear_struct(m, nnz, objIndices, objCoefs)
end
function KN_add_obj_linear_struct(m::Model, objindex::Int, objCoefs::Cdouble)
    return KN_add_obj_linear_struct(m, Cint[objindex], [objCoefs])
end

function KN_add_obj_quadratic_struct(
    m::Model,
    indexVars1::Vector{Cint},
    indexVars2::Vector{Cint},
    coefs::Vector{Cdouble},
)
    nnz = length(indexVars1)
    @assert nnz == length(indexVars2) == length(coefs)
    return KN_add_obj_quadratic_struct(m, nnz, indexVars1, indexVars2, coefs)
end

#=
    CONSTRAINTS
=#
function KN_add_cons(m::Model, ncons::Integer)
    indexes = zeros(Cint, ncons)
    KN_add_cons(m, ncons, indexes)
    return indexes
end
function KN_add_con(m::Model)
    index = Cint[0]
    KN_add_con(m, index)
    return index[1]
end

@kn_get_values get_con_lobnd Cdouble
@kn_get_values get_con_upbnd Cdouble
@kn_get_values get_con_eqbnd Cdouble
@kn_get_values get_con_dual_value Cdouble
@kn_get_values get_con_value Cdouble

function KN_add_con_linear_struct(
    m::Model,
    index_cons::Vector{Cint},
    index_vars::Vector{Cint},
    coefs::Vector{Cdouble},
)
    @assert length(index_cons) == length(index_vars) == length(coefs)
    nnz = length(index_cons)
    return KN_add_con_linear_struct(m, nnz, index_cons, index_vars, coefs)
end
function KN_add_con_linear_struct(
    m::Model,
    index_con::Integer,
    index_vars::Vector{Cint},
    coefs::Vector{Cdouble},
)
    @assert length(index_vars) == length(coefs)
    nnz = length(index_vars)
    return KN_add_con_linear_struct_one(m, nnz, index_con, index_vars, coefs)
end
function KN_add_con_linear_struct(
    m::Model,
    index_con::Integer,
    index_var::Integer,
    coef::Cdouble,
)
    nnz = 1
    return KN_add_con_linear_struct_one(m, nnz, index_con, [index_var], [coef])
end

function KN_add_con_quadratic_struct(
    m::Model,
    index_cons::Vector{Cint},
    index_vars1::Vector{Cint},
    index_vars2::Vector{Cint},
    coefs::Vector{Cdouble},
)
    @assert length(index_cons) ==
            length(index_vars1) ==
            length(index_vars2) ==
            length(coefs)
    nnz = length(index_cons)
    return KN_add_con_quadratic_struct(m, nnz, index_cons, index_vars1, index_vars2, coefs)
end
function KN_add_con_quadratic_struct(
    m::Model,
    index_con::Integer,
    index_vars1::Vector{Cint},
    index_vars2::Vector{Cint},
    coefs::Vector{Cdouble},
)
    @assert length(index_vars1) == length(index_vars2) == length(coefs)
    nnz = length(index_vars1)
    return KN_add_con_quadratic_struct_one(
        m,
        nnz,
        index_con,
        index_vars1,
        index_vars2,
        coefs,
    )
end
function KN_add_con_quadratic_struct(
    m::Model,
    index_con::Integer,
    index_var1::Integer,
    index_var2::Integer,
    coef::Cdouble,
)
    nnz = 1
    return KN_add_con_quadratic_struct_one(
        m,
        nnz,
        index_con,
        Cint[index_var1],
        Cint[index_var2],
        [coef],
    )
end

function KN_set_compcons(
    m::Model,
    ccTypes::Vector{Cint},
    indexComps1::Vector{Cint},
    indexComps2::Vector{Cint},
)
    # get number of constraints
    nnc = length(ccTypes)
    @assert nnc == length(indexComps1) == length(indexComps2)
    return KN_set_compcons(m, nnc, ccTypes, indexComps1, indexComps2)
end

#=
    RESIDUALS
=#
function KN_add_rsds(m::Model, ncons::Integer)
    indexes = zeros(Cint, ncons)
    KN_add_rsds(m, ncons, indexes)
    return indexes
end
function KN_add_rsd(m::Model)
    index = Cint[0]
    KN_add_rsd(m, index)
    return index
end

function KN_add_rsd_linear_struct(
    m::Model,
    indexRsds::Vector{Cint},
    indexVars::Vector{Cint},
    coefs::Vector{Cdouble},
)
    nnz = length(indexRsds)
    @assert nnz == length(indexVars) == length(coefs)
    return KN_add_rsd_linear_struct(m, nnz, indexRsds, indexVars, coefs)
end
function KN_add_rsd_linear_struct(
    m::Model,
    indexRsd::Integer,
    indexVar::Vector{Cint},
    coefs::Vector{Cdouble},
)
    nnz = length(indexVar)
    @assert nnz == length(coefs)
    return KN_add_rsd_linear_struct_one(m, nnz, indexRsd, indexVar, coefs)
end

#=
    SOLVE
=#
function KN_solve(m::Model)
    # Check sanity. If model has Julia callbacks, we need to ensure
    # that Knitro is not multithreaded. Otherwise, the code will segfault
    # as we have trouble calling Julia code from multithreaded C
    # code. See issue #93 on https://github.com/jump-dev/KNITRO.jl.
    if has_callbacks(m)
        if KNITRO_VERSION >= v"13.0"
            KN_set_param(m, KN_PARAM_MS_NUMTHREADS, 1)
            KN_set_param(m, KN_PARAM_NUMTHREADS, 1)
            KN_set_param(m, KN_PARAM_MIP_NUMTHREADS, 1)
        else
            KN_set_param(m, "par_numthreads", 1)
            KN_set_param(m, "par_msnumthreads", 1)
        end
    end
    # For KN_solve, we do not return an error if ret is different of 0.
    m.status = KN_solve(m.env)
    return m.status
end

#=
    GETTERS
=#
function KN_get_solution(m::Model)
    # we first check that the model is well defined to avoid segfault
    @assert m.env != C_NULL
    nx = KN_get_number_vars(m)
    nc = KN_get_number_cons(m)

    x = zeros(Cdouble, nx)
    lambda = zeros(Cdouble, nx + nc)
    status = Cint[0]
    obj = Cdouble[0.0]

    KN_get_solution(m, status, obj, x, lambda)

    # Keep solution in cache.
    m.status = status[1]
    m.x = x
    m.mult = lambda
    m.obj_val = obj[1]
    return status[1], obj[1], x, lambda
end

# some wrapper functions for MOI
function get_status(m::Model)
    @assert m.env != C_NULL
    if m.status != 1
        return m.status
    end
    status = Cint[0]
    obj = Cdouble[0.0]
    KN_get_solution(m, status, obj, C_NULL, C_NULL)
    # Keep status in cache.
    m.status = status[1]
    return status[1]
end

function get_objective(m::Model)
    @assert m.env != C_NULL
    if isfinite(m.obj_val)
        return m.obj_val
    end
    status = Cint[0]
    obj = Cdouble[0.0]
    KN_get_solution(m, status, obj, C_NULL, C_NULL)
    # Keep objective value in cache.
    m.obj_val = obj[1]
    return obj[1]
end

function get_solution(m::Model)
    # We first check that the model is well defined to avoid segfault.
    @assert m.env != C_NULL
    if !isempty(m.x)
        return m.x
    end
    nx = KN_get_number_vars(m)
    x = zeros(Cdouble, nx)
    status = Cint[0]
    obj = Cdouble[0.0]
    KN_get_solution(m, status, obj, x, C_NULL)
    # Keep solution in cache.
    m.x = x
    return x
end
get_solution(m::Model, ix::Int) = isempty(m.x) ? get_solution(m)[ix] : m.x[ix]

function get_dual(m::Model)
    # we first check that the model is well defined to avoid segfault
    @assert m.env != C_NULL
    if !isempty(m.mult)
        return m.mult
    end
    nx = KN_get_number_vars(m)
    nc = KN_get_number_cons(m)
    lambda = zeros(Cdouble, nx + nc)
    status = Cint[0]
    obj = Cdouble[0.0]
    KN_get_solution(m, status, obj, C_NULL, lambda)
    # Keep multipliers in cache.
    m.mult = lambda
    return lambda
end
get_dual(m::Model, ix::Int) = isempty(m.mult) ? get_dual(m)[ix] : m.mult[ix]

function KN_get_objgrad_values(m::Model)
    nnz = KN_get_objgrad_nnz(m)
    indexVars = zeros(Cint, nnz)
    objGrad = zeros(Cdouble, nnz)
    KN_get_objgrad_values(m, indexVars, objGrad)
    return indexVars, objGrad
end

function KN_get_jacobian_values(m::Model)
    nnz = KN_get_jacobian_nnz(m)
    jacvars = zeros(Cint, nnz)
    jaccons = zeros(Cint, nnz)
    jaccoef = zeros(Cdouble, nnz)
    KN_get_jacobian_values(m, jacvars, jaccons, jaccoef)
    return jacvars, jaccons, jaccoef
end

function KN_get_rsd_jacobian_values(m::Model)
    nnz = KN_get_rsd_jacobian_nnz(m)
    jacvars = zeros(Cint, nnz)
    jaccons = zeros(Cint, nnz)
    jaccoef = zeros(Cdouble, nnz)
    KN_get_rsd_jacobian_values(m, jacvars, jaccons, jaccoef)
    return jacvars, jaccons, jaccoef
end

function KN_get_hessian_values(m::Model)
    nnz = KN_get_hessian_nnz(m)
    indexVars1 = zeros(Cint, nnz)
    indexVars2 = zeros(Cint, nnz)
    hess = zeros(Cdouble, nnz)
    KN_get_hessian_values(m, indexVars1, indexVars2, hess)
    return indexVars1, indexVars2, hess
end

function KN_get_var_viols(kc::Model, index::Vector{Cint})
    bndInfeas = zeros(Cint, length(index))
    intInfeas = zeros(Cint, length(index))
    viols = zeros(Cdouble, length(index))
    KN_get_var_viols(kc, length(index), index, bndInfeas, intInfeas, viols)
    return bndInfeas, intInfeas, viols
end

function KN_get_con_viols(kc::Model, index::Vector{Cint})
    infeas = zeros(Cint, length(index))
    viols = zeros(Cdouble, length(index))
    KN_get_con_viols(kc, length(index), index, infeas, viols)
    return infeas, viols
end

function KN_get_presolve_error(m::Model)
    @assert m.env != C_NULL
    component, index, error = Cint[0], Cint[0], Cint[0]
    viol = Cdouble[0.0]
    KN_get_presolve_error(m, component, index, error, viol)
    return convert(Bool, component[1]),
    convert(Int64, index[1]),
    convert(Int64, error[1]),
    convert(Float64, viol[1])
end

@kn_get_attribute get_number_vars Cint
@kn_get_attribute get_number_cons Cint
@kn_get_attribute get_obj_value Cdouble
@kn_get_attribute get_obj_type Cint

@kn_get_attribute get_number_iters Cint
@kn_get_attribute get_number_cg_iters Cint
@kn_get_attribute get_abs_feas_error Cdouble
@kn_get_attribute get_rel_feas_error Cdouble
@kn_get_attribute get_abs_opt_error Cdouble
@kn_get_attribute get_rel_opt_error Cdouble
@kn_get_attribute get_objgrad_nnz Cint
@kn_get_attribute get_jacobian_nnz KNLONG
@kn_get_attribute get_rsd_jacobian_nnz KNLONG
@kn_get_attribute get_hessian_nnz KNLONG
@kn_get_attribute get_solve_time_cpu Cdouble
@kn_get_attribute get_solve_time_real Cdouble

@kn_get_attribute get_mip_number_nodes Cint
@kn_get_attribute get_mip_number_solves Cint
@kn_get_attribute get_mip_abs_gap Cdouble
@kn_get_attribute get_mip_rel_gap Cdouble
@kn_get_attribute get_mip_incumbent_obj Cdouble
@kn_get_attribute get_mip_relaxation_bnd Cdouble
@kn_get_attribute get_mip_lastnode_obj Cdouble

#=
PARAMS
=#
function KN_set_param(m::Model, id::Integer, value::Integer)
    return KN_set_int_param(m, id, value)
end
function KN_set_param(m::Model, param::AbstractString, value::Integer)
    return KN_set_int_param_by_name(m, param, value)
end
function KN_set_param(m::Model, id::Integer, value::Cdouble)
    return KN_set_double_param(m, id, value)
end
function KN_set_param(m::Model, param::AbstractString, value::Cdouble)
    return KN_set_double_param_by_name(m, param, value)
end
function KN_set_param(m::Model, id::Integer, value::AbstractString)
    return KN_set_char_param(m, id, value)
end
function KN_set_param(m::Model, param::AbstractString, value::AbstractString)
    return KN_set_char_param_by_name(m, param, value)
end

function KN_get_int_param(m::Model, id::Integer)
    res = Cint[0]
    KN_get_int_param(m, id, res)
    return res[1]
end
function KN_get_int_param(m::Model, param::AbstractString)
    res = Cint[0]
    KN_get_int_param_by_name(m, param, res)
    return res[1]
end

function KN_get_double_param(m::Model, id::Integer)
    res = Cdouble[0.0]
    KN_get_double_param(m, id, res)
    return res[1]
end
function KN_get_double_param(m::Model, param::AbstractString)
    res = Cdouble[0.0]
    KN_get_double_param_by_name(m, param, res)
    return res[1]
end

function KN_get_param_name(m::Model, id::Integer)
    output_size = 128
    res = " "^output_size
    KN_get_param_name(m, id, res, output_size)
    return _format_output(res)
end

function KN_get_param_doc(m::Model, id::Integer)
    output_size = 128
    res = " "^output_size
    KN_get_param_doc(m, id, res, output_size)
    return _format_output(res)
end

function KN_get_param_type(m::Model, id::Integer)
    res = Cint[0]
    KN_get_param_type(m, id, res)
    return res[1]
end

function KN_get_num_param_values(m::Model, id::Integer)
    res = Cint[0]
    KN_get_num_param_values(m, id, res)
    return res[1]
end

function KN_get_param_value_doc(m::Model, id::Integer, value_id::Integer)
    output_size = 128
    res = " "^output_size
    KN_get_param_value_doc(m, id, value_id, res, output_size)
    return _format_output(res)
end

function KN_get_param_id(m::Model, name::AbstractString)
    res = Cint[0]
    KN_get_param_id(m, name, res)
    return res[1]
end

#=
    NAMES
=#
function KN_get_var_names(m::Model, max_length=1024)
    return String[
        KN_get_var_names(m, Cint(id - 1), max_length) for id in 1:KN_get_number_vars(m)
    ]
end
function KN_get_var_names(m::Model, index::Vector{Cint}, max_length=1024)
    return String[KN_get_var_names(m, id, max_length) for id in index]
end
function KN_get_var_names(m::Model, index::Cint, max_length=1024)
    rawname = zeros(Cchar, max_length)
    ret = KN_get_var_name(m, index, max_length, rawname)
    return String(strip(String(convert(Vector{UInt8}, rawname)), '\0'))
end

function KN_get_con_names(m::Model, max_length=1024)
    return String[
        KN_get_con_names(m, Cint(id - 1), max_length) for id in 1:KN_get_number_cons(m)
    ]
end
function KN_get_con_names(m::Model, index::Vector{Cint}, max_length=1024)
    return String[KN_get_con_names(m, id, max_length) for id in index]
end
function KN_get_con_names(m::Model, index::Cint, max_length=1024)
    rawname = zeros(Cchar, max_length)
    ret = KN_get_con_name(m, index, max_length, rawname)
    return String(strip(String(convert(Vector{UInt8}, rawname)), '\0'))
end


text:
module ODBCLoader

using Octo.Repo: ExecuteResult

# https://github.com/JuliaDatabases/ODBC.jl
using ODBC # v0.8.1

const current = Dict{Symbol, Any}(
    :dsn => nothing,
)

current_dsn() = current[:dsn]

# db_connect
function db_connect(; kwargs...)
    if !isempty(kwargs)
        connectionstring = get(kwargs, :dsn, "")
        username = get(kwargs, :username, "")
        password = get(kwargs, :password, "")
        dsn = ODBC.DSN(connectionstring, username, password)
        current[:dsn] = dsn
    end
end

# db_disconnect
function db_disconnect()
    dsn = current_dsn()
    ODBC.disconnect!(dsn)
    current[:dsn] = nothing
end

# query
function query(sql::String)
    dsn = current_dsn()
    q = ODBC.Query(dsn, sql)
    collect(q)
end

function query(prepared::String, vals::Vector)
    dsn = current_dsn()
    stmt = ODBC.prepare(dsn, prepared)
    ODBC.execute!(stmt, vals)
    ExecuteResult()
end

# execute
function execute(sql::String)::ExecuteResult
    dsn = current_dsn()
    ODBC.execute!(dsn, sql)
    ExecuteResult()
end

function execute(prepared::String, vals::Vector)::ExecuteResult
    dsn = current_dsn()
    stmt = ODBC.prepare(dsn, prepared)
    ODBC.execute!(stmt, vals)
    ExecuteResult()
end

function execute(prepared::String, nts::Vector{<:NamedTuple})::ExecuteResult
    dsn = current_dsn()
    stmt = ODBC.prepare(dsn, prepared)
    for tup in nts
        ODBC.execute!(stmt, collect(tup))
    end
    ExecuteResult()
end

end # module Octo.Backends.ODBCLoader


text:
using RiskAdjustedLinearizations, JLD2, Test

output  = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "reference", "compute_psi_output.jld2"), "r")
GAM1    = output["GAM1"]
GAM2    = output["GAM2"]
GAM3    = output["GAM3"]
GAM4    = output["GAM4"]
GAM5    = output["GAM5"]
GAM6    = output["GAM6"]
JV      = output["JV"]
Psi     = output["Psi"]
Psi_det = output["Psi_det"]

Nzy = sum(size(GAM5))
AA = Matrix{Complex{Float64}}(undef, Nzy, Nzy)
BB = similar(AA)
@testset "QZ decomposition for Ψ" begin
    # Out-of-place
    @test RiskAdjustedLinearizations.compute_Ψ(GAM1, GAM2, GAM3, GAM4, GAM5, GAM6, JV)              ≈ Psi
    @test RiskAdjustedLinearizations.compute_Ψ(GAM1, GAM2, GAM3, GAM4, GAM5, GAM6)                  ≈ Psi_det
    @test RiskAdjustedLinearizations.compute_Ψ(GAM1, GAM2, GAM3, GAM4, GAM5, GAM6, zeros(size(JV))) ≈ Psi_det

    # In-place
    @test RiskAdjustedLinearizations.compute_Ψ!(AA, BB, GAM1, GAM2, GAM3, GAM4, GAM5, GAM6, JV)              ≈ Psi
    @test RiskAdjustedLinearizations.compute_Ψ!(AA, BB, GAM1, GAM2, GAM3, GAM4, GAM5, GAM6)                  ≈ Psi_det
    @test RiskAdjustedLinearizations.compute_Ψ!(AA, BB, GAM1, GAM2, GAM3, GAM4, GAM5, GAM6, zeros(size(JV))) ≈ Psi_det
end


text:
import Base: GenericIOBuffer

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tokenize.Tokens.iskeyword, (Tokenize.Tokens.Kind,))
    precompile(Tokenize.Tokens.isliteral, (Tokenize.Tokens.Kind,))
    precompile(Tokenize.Tokens.isoperator, (Tokenize.Tokens.Kind,))
    precompile(Tokenize.Tokens.Token, (Tokenize.Tokens.Kind,Tuple{Int,Int},Tuple{Int,Int},Int,Int,String))
    precompile(Tokenize.Tokens.Token, ())
    precompile(Tokenize.Tokens.kind, (Tokenize.Tokens.Token,))

    precompile(Tokenize.Lexers.is_identifier_start_char, (Char,))
    precompile(Tokenize.Lexers.readchar, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.next_token, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))

    precompile(Tokenize.Lexers.ishex, (Char,))
    precompile(Tokenize.Lexers.isbinary, (Char,))
    precompile(Tokenize.Lexers.isoctal, (Char,))
    precompile(Tokenize.Lexers.iswhitespace, (Char,))
    precompile(Tokenize.Lexers.Lexer, (String,))
    precompile(Tokenize.Lexers.Lexer, (GenericIOBuffer{Array{UInt8, 1}},))
    precompile(Tokenize.Lexers.tokenize, (String,))

    precompile(Tokenize.Lexers.iterate, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, Bool,))
    precompile(Tokenize.Lexers.iterate, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, Bool,))
    precompile(Tokenize.Lexers.startpos, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.startpos!, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},Int))
    precompile(Tokenize.Lexers.start_token!, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))

    precompile(Tokenize.Lexers.lex_greater, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_prime, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_digit, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},Tokenize.Tokens.Kind))
    precompile(Tokenize.Lexers.lex_identifier, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, Char,))
    precompile(Tokenize.Lexers.lex_less, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_forwardslash, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_minus, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_xor, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_equal, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_bar, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_quote, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_plus, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_dot, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_exclaim, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_colon, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_percent, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_comment, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_comment, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},Bool))
    precompile(Tokenize.Lexers.lex_division, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_circumflex, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_backslash, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_star, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))
    precompile(Tokenize.Lexers.lex_amper, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},))

    precompile(Tokenize.Lexers.lex_whitespace, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},Char))

    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, Char,))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, String,))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},typeof( Base.isdigit),))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},typeof( Tokenize.Lexers.iswhitespace),))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},typeof( Tokenize.Lexers.is_identifier_char),))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}},typeof(Tokenize.Lexers.ishex),))
    precompile(Tokenize.Lexers.accept_batch, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, typeof(Tokenize.Lexers.iswhitespace),))
    precompile(Tokenize.Lexers.accept_batch, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, typeof(Tokenize.Lexers.isdigit),))

    precompile(Tokenize.Lexers.accept_batch, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, typeof(Tokenize.Lexers.iswhitespace),))
    precompile(Tokenize.Lexers.accept_batch, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, typeof(Tokenize.Lexers.isdigit),))
    precompile(Tokenize.Lexers.accept, (Tokenize.Lexers.Lexer{GenericIOBuffer{Array{UInt8, 1}}}, Char,))

    precompile(Tokenize.Lexers.readchar, (GenericIOBuffer{Array{UInt8, 1}},))
end


text:
function vanilla_mcmc(;
  model_args=_model_args,
  burn_in::Int=100,
  chain_length::Int=100,
  verbose=false,
)
  samples = similar(1:chain_length, verbose ? Trace : PredatorPreyParameter)
  trace, = generate(predator_prey_model, model_args, data)
  for i=1:(burn_in+chain_length)
    trace, =  mh(trace, select(:predator_prey_parameter))
    if i > burn_in
      samples[i-burn_in] = verbose ? trace : trace[:predator_prey_parameter]
    end
  end
  samples
end

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
"""
    function latticeCompactHexagonal(box,r)

Make a lattice with compact hexagonal structure in box with spheres of radius r.

# Returns
[posX,posY,posZ] Where pos are arrays with the positions of the cells for each component.
"""
function compactHexagonal(box::Array{<:Real,2},r::Number)

    d = size(box)
    if d != (3,2)
        error("Incompatible number of dimensions ", d)
    end

    #Make first dimension
    lineX = Array(box[1,1]:2*r:box[1,2])

    #Make second dimension
    lineY = fill(box[2,1]+r,length(lineX))
    nY = ceil(INT,(box[2,2]-box[2,1])/(4*sin(pi/3)*r))
    areaX = zeros(FLOAT,1)
    areaY = zeros(FLOAT,1)
    dx = 2*r*cos(pi/3)
    dy = 2*r*sin(pi/3)
    for i in 0:nY-1
        append!(areaX,lineX)
        append!(areaX,lineX.+dx)

        append!(areaY,lineY.+dy*2*i)
        append!(areaY,lineY.+dy*2*i.+dy)
    end

    #Make third dimension
    areaZ = fill(box[2,1]+r,length(areaX))
    nZ = ceil(INT,(box[3,2]-box[3,1])/(4*sin(pi/3)*r))
    volumeX = zeros(FLOAT,1)
    volumeY = zeros(FLOAT,1)
    volumeZ = zeros(FLOAT,1)
    dy = 2*r/2/cos(pi/6)
    dz = 2*sqrt(2/3)*r
    for i in 0:nZ-1
        append!(volumeX,areaX)
        append!(volumeX,areaX)

        append!(volumeY,areaY)
        append!(volumeY,areaY.+dy)

        append!(volumeZ,areaZ.+dz*2*i)
        append!(volumeZ,areaZ.+dz*2*i.+dz)
    end

    #Return
    return [volumeX volumeY volumeZ]

end


text:
"""
    is_leap_year(year)

Return `true` if `year` is a leap year in the gregorian calendar.

on every year that is evenly divisible by 4
  except every year that is evenly divisible by 100
    unless the year is also evenly divisible by 400
"""
function is_leap_year(year)
  if mod(year, 4) == 0
    if mod(year, 100) == 0
      if mod(year, 400) == 0
        return true
      else
        return false
      end
    else
      return true
    end
  else
    return false
  end
end


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)

function get_control_or_attribute(prob::XpressProblem, control::Integer)
    if control in INTEGER_CONTROLS_VALUES
        return getintcontrol(prob, Int32(control))
    elseif control in DOUBLE_CONTROLS_VALUES
        return getdblcontrol(prob, Int32(control))
    elseif control in STRING_CONTROLS_VALUES
        return getstrcontrol(prob, Int32(control))
    elseif control in INTEGER_ATTRIBUTES_VALUES
        return getintattrib(prob, Int32(control))
    elseif control in DOUBLE_ATTRIBUTES_VALUES
        return getdblattrib(prob, Int32(control))
    elseif control in STRING_ATTRIBUTES_VALUES
        return getstrattrib(prob, Int32(control))
    else
        error("Unrecognized parameter: $(control).")
    end
end
function get_control_or_attribute(prob::XpressProblem, control::String)
    control_index = get(INTEGER_CONTROLS, control, -1)
    if control_index != -1
        return getintcontrol(prob, control_index)
    end
    control_index = get(DOUBLE_CONTROLS, control, -1)
    if control_index != -1
        return getdblcontrol(prob, control_index)
    end
    control_index = get(STRING_CONTROLS, control, -1)
    if control_index != -1
        return getstrcontrol(prob, control_index)
    end
    control_index = get(INTEGER_ATTRIBUTES, control, -1)
    if control_index != -1
        return getintattrib(prob, control_index)
    end
    control_index = get(DOUBLE_ATTRIBUTES, control, -1)
    if control_index != -1
        return getdblattrib(prob, control_index)
    end
    control_index = get(STRING_ATTRIBUTES, control, -1)
    if control_index != -1
        return getstrattrib(prob, control_index)
    end
    error("Unrecognized control parameter: $(control).")
end

"""
    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: dispatch on Val(control) instead?
    if control in INTEGER_CONTROLS_VALUES
        return getintcontrol(prob, Int32(control))
    elseif control in DOUBLE_CONTROLS_VALUES
        return getdblcontrol(prob, Int32(control))
    elseif control in STRING_CONTROLS_VALUES
        return getstrcontrol(prob, Int32(control))
    else
        error("Unrecognized control parameter: $(control).")
    end
end
function getcontrol(prob::XpressProblem, control::String)
    control_index = get(INTEGER_CONTROLS, control, -1)
    if control_index != -1
        return getintcontrol(prob, control_index)
    end
    control_index = get(DOUBLE_CONTROLS, control, -1)
    if control_index != -1
        return getdblcontrol(prob, control_index)
    end
    control_index = get(STRING_CONTROLS, control, -1)
    if control_index != -1
        return getstrcontrol(prob, control_index)
    end
    error("Unrecognized control parameter: $(control).")
end

# getcontrol(prob::XpressProblem, control::Symbol) = getcontrol(prob, getproperty(Lib, Symbol("XPRS_$(String(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::Integer, val) = setcontrol!(prob, Cint(control), val)
function setcontrol!(prob::XpressProblem, control::Cint, val)
    # TODO: dispatch on Val(control) instead?
    if control in INTEGER_CONTROLS_VALUES
        if isinteger(val)
            setintcontrol(prob, control, Int32(val))
        else
            error("Expected and integer and got $val")
        end
    elseif control in DOUBLE_CONTROLS_VALUES
        setdblcontrol(prob, control, Float64(val))
    elseif control in STRING_CONTROLS_VALUES
        setstrcontrol(prob, control, val)
    else
        error("Unrecognized control parameter: $(control).")
    end
end
function setcontrol!(prob::XpressProblem, control::String, val)
    control_index = get(INTEGER_CONTROLS, control, -1)
    if control_index != -1
        if isinteger(val)
            return setintcontrol(prob, control_index, Int32(val))
        else
            error("Expected and integer and got $val")
        end
    end
    control_index = get(DOUBLE_CONTROLS, control, -1)
    if control_index != -1
        return setdblcontrol(prob, control_index, Float64(val))
    end
    control_index = get(STRING_CONTROLS, control, -1)
    if control_index != -1
        return setstrcontrol(prob, control_index, val)
    end
    error("Unrecognized control parameter: $(control).")
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

# 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:
using Distributions

mdp = SimpleGridWorld(tprob = 1)
hist = simulate(HistoryRecorder(), mdp, RandomPolicy(mdp), GWPos(3,3))

## Test the playback simulator
playback = PlaybackPolicy(collect(action_hist(hist)), RandomPolicy(mdp))
@test all(playback.actions .== action_hist(hist))
@test playback.backup_policy isa RandomPolicy
@test playback.i == 1
@test_throws AssertionError logpdf(playback, hist)

hist2 = simulate(HistoryRecorder(), mdp, playback, GWPos(3,3))
@test hist == hist2

## Test log probability
Distributions.logpdf(p::RandomPolicy, h) = length(h)*log(1. / length(actions(p.problem)))
playback = PlaybackPolicy(collect(action_hist(hist)), RandomPolicy(mdp), logpdfs = -ones(length(hist)))
hist2 = simulate(HistoryRecorder(), mdp, playback, GWPos(3,3))
@test logpdf(playback, hist2) == -1*length(hist2)

playback = PlaybackPolicy([], RandomPolicy(mdp))
@test logpdf(playback, hist2) == length(hist2)*log(0.25)

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
include("../example/example_setup.jl")

using Test, Vecchia, SparseArrays

# Create a vecc object that uses enough block-conditioning points that the
# likelihood evaluation is exact.
const vecc_exact = Vecchia.kdtreeconfig(sim, pts, 64, 10000, kfn)

# Test 1: nll gives the exact likelihood for a vecchia config where the
# conditioning set is every prior point.
vecchia_nll  = nll(vecc_exact, ones(2))
debug_nll    = Vecchia.exact_nll(vecc_exact, ones(2))
@test isapprox(vecchia_nll, debug_nll)

# Test 2: precisionmatrix gives the exact precision when it should....modulo
# some numerical noise, which is a little bit higher than I would expect.
pmat = Vecchia.precisionmatrix(vecc_exact, ones(2))
ptsv = reduce(vcat, vecc_exact.pts)
cmat = [kfn(x,y,ones(2)) for x in ptsv, y in ptsv]
@test maximum(abs, inv(Matrix(pmat))-cmat) < 1e-5

# Test 3: the nll agrees with the one obtained with the Vecchia precision.
datv = reduce(vcat, vecc_exact.data)
pnll = Vecchia.nll_precision(pmat, datv)
@test isapprox(pnll, debug_nll)

# Test 4: the scalarized nll agrees with the non-scalarized one.
scalarized_nll    = Vecchia.nll(vecc_s, ones(2))
nonscalarized_nll = Vecchia.nll(vecc, ones(2))
@test isapprox(scalarized_nll, nonscalarized_nll)

# Test 5: the nll with multiple data sources agrees with the sum of two
# single-data nlls.
const new_data  = range(0.0, 1.0, length=length(sim))
const joint_cfg = Vecchia.kdtreeconfig(hcat(sim, new_data), pts, 64, 3, kfn)
const new_cfg   = Vecchia.kdtreeconfig(new_data, pts, 64, 3, kfn)
@test isapprox(Vecchia.nll(joint_cfg, ones(2)),
               Vecchia.nll(vecc, ones(2)) + Vecchia.nll(new_cfg, ones(2)))

# Test 6: make sure that with multiple data the precision matrices are the same.
const pmat_1 = Vecchia.precisionmatrix(vecc, ones(2))
const pmat_2 = Vecchia.precisionmatrix(joint_cfg, ones(2))
@test isapprox(pmat_1, pmat_2)

# Test 7: confirm that the rchol-based nll is equal to the standard nll.
rchol_nll = Vecchia.nll_rchol(vecc, ones(2))
@test isapprox(rchol_nll, Vecchia.nll(vecc, ones(2)))

# Test 8: the same, but with the "sparse" version of the
const _U = sparse(Vecchia.rchol(vecc, ones(2)))
const _S = Symmetric(_U*_U')
@test isapprox(rchol_nll, Vecchia.nll_precision(_S, datv))


text:
using Swizzles
using Test

@testset "Swizzles" begin
    include("Swizzles.jl")
    include("ShallowArrays.jl")
    include("EachindexArrays.jl")
    include("ExtrudedArrays.jl")
    include("ValArrays.jl")
    include("NamedArrays.jl")
    include("Virtuals.jl")
    include("ArrayifiedArrays.jl")
    include("StylishArrays.jl")
    include("PermutedArrays.jl")
    include("util.jl")
end


text:
using POMDPModels
using POMDPToolbox
using GenerativeModels
using MCTS

mdp = GridWorld()

rmdp = RobustAdapter(mdp)

solver = RobustMCTSSolver(rollout_nature=mdp, c=10.0)

policy = solve(solver, rmdp)

sim = RolloutSimulator(max_steps=100)
simulate(sim, mdp, policy, initial_state(mdp, MersenneTwister()))

solver2 = RobustMCTSSolver(rng=MersenneTwister(), n_iterations=1000, rollout_nature=mdp, c=10.0, depth=100)
policy2 = solve(solver2, rmdp)
mcts_solver = DPWSolver(rng=MersenneTwister(), n_iterations=1000, exploration_constant=10.0, depth=100)
mcts_policy = solve(mcts_solver, mdp)

rng = MersenneTwister(12)
N = 100
# total_rank_difference = 0
# nb_misranked = 0
total_reward_difference = 0.0
total_reward = 0.0
abs_total_reward_difference = 0.0
for i in 1:N
    s = initial_state(mdp, rng)
#     a2 = action(policy2, s)
#     mcts_a = action(mcts_policy, s)
#     print('.')
#     if a2 != mcts_a
#         nb_misranked += 1
#         ranked = MCTS.ranked_actions(policy2, s)
#         total_rank_difference += findfirst(ranked, mcts_a)
#     end
    policy2 = solve(solver2, rmdp)
    mcts_policy = solve(mcts_solver, mdp)
    sim2 = RolloutSimulator(rng=MersenneTwister(i), max_steps=100)
    sim3 = RolloutSimulator(rng=MersenneTwister(i), max_steps=100)
    r2 = simulate(sim2, mdp, policy2, s)
    r3 = simulate(sim3, mdp, mcts_policy, s)
    total_reward += r2
    total_reward_difference += r3-r2
    abs_total_reward_difference += abs(r3-r2)
    print('.')
end

# @show nb_misranked
# @show average_rank_difference = total_rank_difference/N

@show total_reward/N
@show total_reward_difference/N
@show abs_total_reward_difference/N