text:
using DataFrames

struct Climate
    climate_data::DataFrame
    rainfall_id::String
    et_id::String
    # t_id::Union{String, Nothing} = nothing
end

"""
    subcatchment_data(node::NetworkNode, climate::Climate)

Extract all data for a given node from climate object.
"""
function subcatchment_data(node::NetworkNode, climate::Climate)::DataFrame
    data = climate.climate_data
    cols = filter(x -> occursin(node.name, string(x)), names(data))

    return data[:, vcat(["Date"], cols)]
end

"""
    rainfall_data(node::NetworkNode, climate::Climate)::DataFrame

Extract rainfall data for a given node.
"""
function rainfall_data(node::NetworkNode, climate::Climate)::DataFrame
    data = climate.climate_data
    rain_col = filter(x -> occursin(node.name, x)
                            & occursin(climate.rainfall_id, x),
                            names(data))[1]

    return data[:, rain_col]
end

"""
    climate_values(node::NetworkNode, climate::Climate, timestep::Union{Nothing, Int}=nothing)

Extract climate related data for a given time step.
"""
function climate_values(node::NetworkNode, climate::Climate,
                        timestep::Union{Nothing, Int}=nothing)
    node_name = node.name

    data = climate.climate_data

    # TODO : Catch instances where data is not found (raises BoundsError)
    rain_col = filter(x -> occursin(node_name, x)
                            & occursin(climate.rainfall_id, x),
                            names(data))[1]
    et_col = filter(x -> occursin(node_name, x)
                            & occursin(climate.et_id, x),
                            names(data))[1]

    if isempty(rain_col) | isempty(et_col)
        throw(ArgumentError("No climate data found for $(node_name) at time step: $(timestep)"))
    end

    if isnothing(timestep)
        return select(data, [rain_col, et_col])
    end

    return data[timestep, [rain_col, et_col]]
end

"""
    sim_length(climate::Climate)::Int64

Simulation length is dependent on available climate data.
"""
function sim_length(climate::Climate)::Int64
    return nrow(climate.climate_data)
end

function timesteps(climate::Climate)::Array
    return climate.climate_data[:, :Date]
end

Base.length(climate::Climate) = nrow(climate.climate_data)


text:
# Autogenerated wrapper script for Qt5Declarative_jll for x86_64-w64-mingw32-cxx03
export libqt5qml, libqt5qmlmodels, libqt5qmlworkerscript, libqt5quick, libqt5quickparticles, libqt5quickshapes, libqt5quicktest, libqt5quickwidgets

using Qt5Base_jll
JLLWrappers.@generate_wrapper_header("Qt5Declarative")
JLLWrappers.@declare_library_product(libqt5qml, "Qt5Qml.dll")
JLLWrappers.@declare_library_product(libqt5qmlmodels, "Qt5QmlModels.dll")
JLLWrappers.@declare_library_product(libqt5qmlworkerscript, "Qt5QmlWorkerScript.dll")
JLLWrappers.@declare_library_product(libqt5quick, "Qt5Quick.dll")
JLLWrappers.@declare_library_product(libqt5quickparticles, "Qt5QuickParticles.dll")
JLLWrappers.@declare_library_product(libqt5quickshapes, "Qt5QuickShapes.dll")
JLLWrappers.@declare_library_product(libqt5quicktest, "Qt5QuickTest.dll")
JLLWrappers.@declare_library_product(libqt5quickwidgets, "Qt5QuickWidgets.dll")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Base_jll)
    JLLWrappers.@init_library_product(
        libqt5qml,
        "bin\\Qt5Qml.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlmodels,
        "bin\\Qt5QmlModels.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlworkerscript,
        "bin\\Qt5QmlWorkerScript.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick,
        "bin\\Qt5Quick.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickparticles,
        "bin\\Qt5QuickParticles.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickshapes,
        "bin\\Qt5QuickShapes.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktest,
        "bin\\Qt5QuickTest.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickwidgets,
        "bin\\Qt5QuickWidgets.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
include("robot_operations.jl")
pygui(true)

module Chess
    using HorizonSideRobots
    import Main.inverseSide, Main.isborder, Main.goToTheCorner, Main.followByCount

    export mark_chess

    pathBack = []

    X=0
    Y=0

    BOARD_SIZE = 0

    function init(r::Robot)
        global pathBack
        pathBack = goToTheCorner(r, West, Sud)
    end

    function mark_chess(r::Robot,n::Int)
        global BOARD_SIZE
        BOARD_SIZE = n
        side=Ost

        init(r)
        mark_row(r,side)
        while isborder(r,Nord)==false
            move_decart!(r,Nord)
            side = inverseSide(side)
            mark_row(r,side)
        end

        goToTheCorner(r, West, Sud)
        for item in pathBack[end:-1:1]
            followByCount(r, item[1], item[2])
        end
    end

    function mark_row(r::Robot, side::HorizonSide)
        putmarker_chess!(r)
        while isborder(r,side)==false
            move_decart!(r,side)
            putmarker_chess!(r)
        end
    end

    function putmarker_chess!(r)
        if (mod(X, 2*BOARD_SIZE) in 0:BOARD_SIZE-1) && (mod(Y, 2*BOARD_SIZE) in 0:BOARD_SIZE-1)
            putmarker!(r)
        end
    end

    function move_decart!(r,side)
        global X, Y
        # Id prefer to use switch-case, bruh
        if side==Nord
            Y+=1
        elseif side==Sud
            Y-=1
        elseif side==Ost
            X+=1
        else
            X-=1
        end
        move!(r,side)
    end

end

r=Robot(animate=true)
using .Chess
mark_chess(r, 2)

max_stars_repo_name: TheDarkMoon/mirea-progs

max_issues_repo_name: TheDarkMoon/mirea-progs

max_forks_repo_name: TheDarkMoon/mirea-progs


text:
function get_delta_i(obj::GLONASSOrbitalElements)
    return jcall(obj, "getDeltaI", jdouble, ())
end

function get_delta_t(obj::GLONASSOrbitalElements)
    return jcall(obj, "getDeltaT", jdouble, ())
end

function get_delta_t_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getDeltaTDot", jdouble, ())
end

function get_e(obj::GLONASSOrbitalElements)
    return jcall(obj, "getE", jdouble, ())
end

function get_gamma_n(obj::GLONASSOrbitalElements)
    return jcall(obj, "getGammaN", jdouble, ())
end

function get_iod(obj::GLONASSOrbitalElements)
    return jcall(obj, "getIOD", jint, ())
end

function get_lambda(obj::GLONASSOrbitalElements)
    return jcall(obj, "getLambda", jdouble, ())
end

function get_n4(obj::GLONASSOrbitalElements)
    return jcall(obj, "getN4", jint, ())
end

function get_na(obj::GLONASSOrbitalElements)
    return jcall(obj, "getNa", jint, ())
end

function get_pa(obj::GLONASSOrbitalElements)
    return jcall(obj, "getPa", jdouble, ())
end

function get_time(obj::GLONASSOrbitalElements)
    return jcall(obj, "getTime", jdouble, ())
end

function get_tn(obj::GLONASSOrbitalElements)
    return jcall(obj, "getTN", jdouble, ())
end

function get_x(obj::GLONASSOrbitalElements)
    return jcall(obj, "getX", jdouble, ())
end

function get_x_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getXDot", jdouble, ())
end

function get_x_dot_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getXDotDot", jdouble, ())
end

function get_y(obj::GLONASSOrbitalElements)
    return jcall(obj, "getY", jdouble, ())
end

function get_y_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getYDot", jdouble, ())
end

function get_y_dot_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getYDotDot", jdouble, ())
end

function get_z(obj::GLONASSOrbitalElements)
    return jcall(obj, "getZ", jdouble, ())
end

function get_z_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getZDot", jdouble, ())
end

function get_z_dot_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getZDotDot", jdouble, ())
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("abstract_density.jl")
include("generic_density.jl")
include("const_density.jl")
include("density_product.jl")
include("distribution_density.jl")
include("posterior_density.jl")
include("parameter_mapped_density.jl")
include("external_density.jl")


text:
function Tables.schema(layer::AbstractFeatureLayer)::Tables.Schema
    geom_names, field_names, featuredefn, fielddefns =
        schema_names(layerdefn(layer))
    ngeom = ArchGDAL.ngeom(featuredefn)
    geom_types =
        (IGeometry{gettype(getgeomdefn(featuredefn, i))} for i in 0:ngeom-1)
    field_types =
        (convert(DataType, gettype(fielddefn)) for fielddefn in fielddefns)
    return Tables.Schema(
        (geom_names..., field_names...),
        (geom_types..., field_types...),
    )
end

Tables.istable(::Type{<:AbstractFeatureLayer})::Bool = true
Tables.rowaccess(::Type{<:AbstractFeatureLayer})::Bool = true

function Tables.rows(layer::T)::T where {T<:AbstractFeatureLayer}
    return layer
end

function Tables.getcolumn(row::Feature, i::Int)
    if i > nfield(row)
        return getgeom(row, i - nfield(row) - 1)
    elseif i > 0
        return getfield(row, i - 1)
    else
        return missing
    end
end

function Tables.getcolumn(row::Feature, name::Symbol)
    field = getfield(row, name)
    if !ismissing(field)
        return field
    end
    geom = getgeom(row, name)
    if geom.ptr != C_NULL
        return geom
    end
    return missing
end

function Tables.columnnames(
    row::Feature,
)::NTuple{Int64(nfield(row) + ngeom(row)),Symbol}
    geom_names, field_names = schema_names(getfeaturedefn(row))
    return (geom_names..., field_names...)
end

function schema_names(featuredefn::IFeatureDefnView)
    fielddefns = (getfielddefn(featuredefn, i) for i in 0:nfield(featuredefn)-1)
    field_names = (Symbol(getname(fielddefn)) for fielddefn in fielddefns)
    geom_names = collect(
        Symbol(getname(getgeomdefn(featuredefn, i - 1))) for
        i in 1:ngeom(featuredefn)
    )
    return (geom_names, field_names, featuredefn, fielddefns)
end


text:
# Testing if we can use XGBoost for surrogate modeling
using XGBoost

N, P, k = (400, 50, 3)
active_coefs = 2 + 10*rand(k)
coefs = shuffle(vcat(active_coefs, zeros(P-k)))
X = randn(N, P)
Y = X * coefs
testX = randn(N, P)
testY = testX * coefs

mse(y, yhat) = mean((y .- yhat).^2)
made(y, yhat) = mean(abs( 100.0 * (y .- yhat) ./ y ))

num_round = 100
@time bst = xgboost(X, num_round, label = Y, max_depth = 3)
predY = predict(bst, testX);
mse(testY, predY)
made(testY, predY)


text:
using LMDB_jll

# Error codes (TODO: define more)
const MDB_SUCCESS = Cint(0)
const MDB_NOTFOUND = Cint(-30798)

# Environment flags (TODO: define more)

const MDB_FIXEDMAP = Cuint(0x01)
const MDB_NOSUBDIR = Cuint(0x4000)
const MDB_NOSYNC = Cuint(0x10000)
const MDB_RDONLY = Cuint(0x20000)
const MDB_WRITEMAP = Cuint(0x80000)
const MDB_MAPASYNC = Cuint(0x100000)
const MDB_NOTLS = Cuint(0x200000)
const MDB_NOLOCK = Cuint(0x400000)
const MDB_NORDAHEAD = Cuint(0x800000)
const MDB_NOMEMINIT = Cuint(0x1000000)

# DBI flags
const MDB_CREATE = Cuint(0x40000)

# Cursor operation constants (TODO: define more)
@enum MDBCursorOp::Cint begin
    MDB_FIRST
    MDB_FIRST_DUP
    MDB_GET_BOTH
    MDB_GET_BOTH_RANGE
    MDB_GET_CURRENT
    MDB_GET_MULTIPLE
    MDB_LAST
    MDB_LAST_DUP
    MDB_NEXT
    MDB_NEXT_DUP
    MDB_NEXT_MULTIPLE
    MDB_NEXT_NODUP
    MDB_PREV
    MDB_PREV_DUP
    MDB_PREV_NODUP
    MDB_SET
    MDB_SET_KEY
    MDB_SET_RANGE
    MDB_PREV_MULTIPLE
end

# Common flag combo
const DEFAULT_ROTXN_FLAGS = MDB_RDONLY | MDB_NOTLS

const Cmode_t = Cushort

macro fieldoffset(T, num)
    return :(fieldoffset($T, $(esc(num))))
end

mutable struct MDBValue{T}
    size::Csize_t    # size of the data item
    ptr::Ptr{Cvoid}  # address of the data item
    data::T
    MDBValue() = new{Nothing}(zero(Csize_t), C_NULL)
    MDBValue{T}() where {T} = new{T}(zero(Csize_t), C_NULL)
    MDBValue(v::String) = new{String}(sizeof(v), pointer(v), v)
    MDBValue(v::Array) = new{Array}(sizeof(v), pointer(v), v)
    function MDBValue(v::T) where {T}
        # alloc new object; set size; copy primitive value or reference
        o = new{T}(sizeof(v), C_NULL, v)
        # mutable struct have static pointer address so its safe to just return that
        # primitive types and immutable structs have been copied into o - point to that instead
        o.ptr = isimmutable(v) ? pointer_from_objref(o) + @fieldoffset(MDBValue{0}, 3) : pointer_from_objref(v)
        return o
    end
end
Base.pointer(val::MDBValue) = pointer_from_objref(val)
function Base.:(==)(a::MDBValue, b::MDBValue)
    a !== b || return true
    a.ptr == b.ptr || return true
    a.size == b.size || return true
    arr_a = unsafe_wrap(Array, reinterpret(Ptr{UInt8}, a.ptr). a.size)
    arr_b = unsafe_wrap(Array, reinterpret(Ptr{UInt8}, b.ptr). b.size)
    return arr_a == arr_b
end

struct MDBStat
    ms_psize::Cuint
    ms_depth::Cuint
    ms_branch_pages::Csize_t
    ms_leaf_pages::Csize_t
    ms_overflow_pages::Csize_t
    ms_entries::Csize_t
    MDBStat() = new(zero(Cuint), zero(Cuint), zero(Csize_t), zero(Csize_t), zero(Csize_t), zero(Csize_t))
end

# MDBValue (output) -> T
Base.convert(::Type{String}, val::MDBValue{Nothing}) = unsafe_string(convert(Ptr{UInt8}, val.ptr), val.size)
Base.convert(::Type{Array{T,N}}, val::MDBValue{Nothing}) where {T,N} = unsafe_wrap(Array, convert(Ptr{T}, val.ptr), val.size)
Base.convert(::Type{T}, val::MDBValue{Nothing}) where {T} = unsafe_load(convert(Ptr{T}, val.ptr))
Base.convert(::Type{Any}, val::MDBValue{Nothing}) = val
Base.convert(::Type{MDBValue}, val::MDBValue{Nothing}) = val
Base.convert(::Type{MDBValue{Nothing}}, val::MDBValue{Nothing}) = val

# T -> MDBValue (input)
Base.convert(::Type{MDBValue}, val::MDBValue) = val
Base.convert(::Type{MDBValue}, val) = MDBValue(val)
Base.convert(::Type{T}, val::MDBValue{T}) where {T} = val.data

struct LMDBError <: Exception
    code::Cint
end
Base.show(io::IO, err::LMDBError) = print(io, "$(mdb_strerror(err.code)) ($(err.code))")

macro chkres(res)
    return :(iszero($(esc(res))) || throw(LMDBError($(esc(res)))))
end

function mdb_strerror(err::Cint)
    errstr = ccall((:mdb_strerror, liblmdb), Cstring, (Cint,), err)
    return unsafe_string(errstr)
end

function mdb_env_create()
    env_ref = Ref{Ptr{Cvoid}}(C_NULL)
    @chkres ccall((:mdb_env_create, liblmdb), Cint, (Ptr{Ptr{Cvoid}},), env_ref)
    return env_ref[]
end

function mdb_env_close(env::Ptr{Cvoid})
    ccall((:mdb_env_close, liblmdb), Cvoid, (Ptr{Cvoid},), env)
    return nothing
end

function mdb_env_open(env::Ptr{Cvoid}, path::String, flags::Cuint, mode::Cmode_t)
    @chkres ccall((:mdb_env_open, liblmdb), Cint, (Ptr{Cvoid}, Cstring, Cuint, Cmode_t), env, path, flags, mode)
    return nothing
end

function mdb_version()
    str = ccall((:mdb_version, liblmdb), Cstring, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), C_NULL, C_NULL, C_NULL)
    return unsafe_string(str)
end

function mdb_txn_begin(env::Ptr{Cvoid}, flags::Cuint)
    txn_ref = Ref{Ptr{Cvoid}}(C_NULL)
    @chkres ccall((:mdb_txn_begin, liblmdb), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cuint, Ptr{Ptr{Cvoid}}), env, C_NULL, flags, txn_ref)
    return txn_ref[]
end

function mdb_txn_abort(txn::Ptr{Cvoid})
    ccall((:mdb_txn_abort, liblmdb), Cvoid, (Ptr{Cvoid},), txn)
    return nothing
end

function mdb_txn_reset(txn::Ptr{Cvoid})
    ccall((:mdb_txn_reset, liblmdb), Cvoid, (Ptr{Cvoid},), txn)
    return nothing
end

function mdb_txn_renew(txn::Ptr{Cvoid})
    @chkres ccall((:mdb_txn_renew, liblmdb), Cint, (Ptr{Cvoid},), txn)
    return nothing
end

function mdb_txn_commit(txn::Ptr{Cvoid})
    @chkres ccall((:mdb_txn_commit, liblmdb), Cint, (Ptr{Cvoid},), txn)
    return
end

function mdb_dbi_open(txn::Ptr{Cvoid}, dbname::String, flags::Cuint)
    handle = Cuint[0]
    cdbname = isempty(dbname) ? convert(Cstring, Ptr{UInt8}(C_NULL)) : dbname
    @chkres ccall((:mdb_dbi_open, liblmdb), Cint, (Ptr{Cvoid}, Cstring, Cuint, Ptr{Cuint}), txn, cdbname, flags, handle)
    return handle[1]
end

function mdb_dbi_close(txn::Ptr{Cvoid}, dbi::Cuint)
    ccall((:mdb_dbi_close, liblmdb), Cvoid, (Ptr{Cvoid}, Cuint), txn, dbi)
    return nothing
end

function mdb_env_set_maxreaders(env::Ptr{Cvoid}, readers::Cuint)
    @chkres ccall((:mdb_env_set_maxreaders, liblmdb), Cint, (Ptr{Cvoid}, Cuint), env, readers)
    return nothing
end

function mdb_env_set_mapsize(env::Ptr{Cvoid}, size::Csize_t)
    @chkres ccall((:mdb_env_set_mapsize, liblmdb), Cint, (Ptr{Cvoid}, Csize_t), env, size)
    return nothing
end

function mdb_env_set_maxdbs(env::Ptr{Cvoid}, dbs::Cuint)
    @chkres ccall((:mdb_env_set_maxdbs, liblmdb), Cint, (Ptr{Cvoid}, Cuint), env, dbs)
    return nothing
end

function mdb_env_sync(env::Ptr{Cvoid}, force::Cint)
    @chkres ccall((:mdb_env_sync, liblmdb), Cint, (Ptr{Cvoid}, Cint), env, force)
    return nothing
end

function mdb_put(txn::Ptr{Cvoid}, dbi::Cuint, key::Ptr{Cvoid}, val::Ptr{Cvoid}, flags::Cuint)
    @chkres ccall((:mdb_put, liblmdb), Cint, (Ptr{Cvoid}, Cuint, Ptr{Cvoid}, Ptr{Cvoid}, Cuint), txn, dbi, key, val, flags)
    return nothing
end

function mdb_del(txn::Ptr{Cvoid}, dbi::Cuint, key::Ptr{Cvoid}, val::Ptr{Cvoid})
    res = ccall((:mdb_del, liblmdb), Cint, (Ptr{Cvoid}, Cuint, Ptr{Cvoid}, Ptr{Cvoid}), txn, dbi, key, val)
    (res == MDB_SUCCESS || res == MDB_NOTFOUND) || throw(LMDBError(res))
    return res != MDB_NOTFOUND
end

function mdb_get!(txn::Ptr{Cvoid}, dbi::Cuint, key::Ptr{Cvoid}, val::Ptr{Cvoid})
    res = ccall((:mdb_get, liblmdb), Cint, (Ptr{Cvoid}, Cuint, Ptr{Cvoid}, Ptr{Cvoid}), txn, dbi, key, val)
    (res == MDB_SUCCESS || res == MDB_NOTFOUND) || throw(LMDBError(res))
    return res != MDB_NOTFOUND
end

function mdb_cursor_open(txn::Ptr{Cvoid}, dbi::Cuint)
    handle = Ptr{Cvoid}[C_NULL]
    @chkres ccall((:mdb_cursor_open, liblmdb), Cint, (Ptr{Cvoid}, Cuint, Ptr{Ptr{Cvoid}}), txn, dbi, handle)
    return handle[1]
end

function mdb_cursor_close(cur::Ptr{Cvoid})
    ccall((:mdb_cursor_close, liblmdb), Cvoid, (Ptr{Cvoid},), cur)
    return nothing
end

function mdb_cursor_txn(cur::Ptr{Cvoid})
    return ccall((:mdb_cursor_txn, liblmdb), Ptr{Cvoid}, (Ptr{Cvoid},), cur)
end

function mdb_cursor_get!(cur::Ptr{Cvoid}, key::Ptr{Cvoid}, val::Ptr{Cvoid}, op::MDBCursorOp)
    res = ccall((:mdb_cursor_get, liblmdb), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Cint), cur, key, val, op)
    (res == MDB_SUCCESS || res == MDB_NOTFOUND) || throw(LMDBError(res))
    return res != MDB_NOTFOUND
end

function mdb_env_stat(env::Ptr{Cvoid})
    statref = Ref(MDBStat())
    @chkres ccall((:mdb_env_stat, liblmdb), Cint, (Ptr{Cvoid}, Ptr{MDBStat}), env, statref)
    return statref[]
end

function mdb_stat(txn::Ptr{Cvoid}, dbi::Cuint)
    statref = Ref(MDBStat())
    @chkres ccall((:mdb_stat, liblmdb), Cint, (Ptr{Cvoid}, Cuint, Ptr{MDBStat}), txn, dbi, statref)
    return statref[]
end

function mdb_txn_id(txn::Ptr{Cvoid})
    return ccall((:mdb_stat, liblmdb), Csize_t, (Ptr{Cvoid},), txn)
end

mutable struct Environment
    handle::Ptr{Cvoid}
    path::String
    rotxn::Vector{Ptr{Cvoid}}
    wlock::ReentrantLock
end

const OPENED_ENVS = Dict{String,Environment}()

function Environment(path::String; flags::Cuint=0x80000|0x100000, mode::Cmode_t=0o755, maxdbs=0, mapsize=10485760, maxreaders=126, rotxnflags::Cuint=DEFAULT_ROTXN_FLAGS)
    env = get(OPENED_ENVS, path, nothing)
    if isnothing(env)
        # create all directories needed to host the data
        mkpath(iszero(flags & MDB_NOSUBDIR) ? path : dirname(path), mode=mode)
        rotxn = [C_NULL for i in 1:Threads.nthreads()]
        env = Environment(mdb_env_create(), path, rotxn, ReentrantLock())
        mdb_env_set_maxdbs(env.handle, convert(Cuint, maxdbs))
        mdb_env_set_mapsize(env.handle,  convert(Csize_t, mapsize))
        mdb_env_set_maxreaders(env.handle, convert(Cuint, maxreaders))
        try
            mdb_env_open(env.handle, path, flags, mode)
        catch e
            mdb_env_close(env.handle)
            rethrow(e)
        end
        OPENED_ENVS[path] = env
        # create read-only transaction handles, which can be used quickly (mdb_txn_renew/mdb_txn_reset use-case)
        for i in 1:Threads.nthreads()
            env.rotxn[i] = mdb_txn_begin(env.handle, rotxnflags)
            mdb_txn_reset(env.rotxn[i])
        end
        # register finalizer for lazy cleanup
        finalizer(close, env)
    end
    return env
end

Base.show(io::IO, e::Environment) = print(io, typeof(e), "(", repr(e.path), ")")

Base.isopen(env::Environment) = env.handle != C_NULL

function Base.close(env::Environment)
    if isopen(env)
        for handle in env.rotxn
            handle == C_NULL || mdb_txn_abort(handle)
        end
        @assert get(OPENED_ENVS, env.path, nothing) == env "OPENED_ENVS does not contain env for $(path) - smells like threading issue ..."
        mdb_env_close(env.handle)
        env.handle = C_NULL
        delete!(OPENED_ENVS, env.path)
        return true
    end
    return false
end

Base.flush(env::Environment) = mdb_env_sync(env.handle, one(Cint))
function Base.flush(func::Function, env::Environment)
    try
        func()
    finally
        mdb_env_sync(env.handle, one(Cint))
    end
end

function __init__()
    @info mdb_version()
    # this will catch changes in struct layout
    @assert fieldname(MDBValue{Nothing}, 3) == :data
end


text:
using NavAbilitySDK
using DistributedFactorGraphs
using IncrementalInference

mutable struct NavAbilityDFG{T <: AbstractParams} <: AbstractDFG{T}
    navabilityClient::NavAbilityClient
    solverParams::T
    userId::String
    robotId::String
    sessionId::String
end

function NavAbilityDFG(apiUrl::String, userId::String, robotId::String, sessionId::String)
    navabilityClient = NavAbilityHttpsClient(apiUrl)
    solverParams = SolverParams()
    solverParams.graphinit = false
    return NavAbilityDFG(navabilityClient,solverParams,userId,robotId,sessionId)
end

function Base.show(io::IO, dfg::NavAbilityDFG)
    summary(io, dfg)
    println(io, "\n  UserId: ", dfg.userId)
    println(io, "  RobotId: ", dfg.robotId)
    println(io, "  SessionId: ", dfg.sessionId)
end
Base.show(io::IO, ::MIME"text/plain", dfg::NavAbilityDFG) = show(io, dfg)


text:
using SafeTestsets

@safetestset "Day02" begin
    using AdventOfCode.TestUtils
    using AdventOfCode.Day02

    @testset "run_program!" begin
        using AdventOfCode.Day02: run_program!
        @test run_program!([1,9,10,3,2,3,11,0,99,30,40,50]).mem.vect == [
            3500,9,10,70,
            2,3,11,0,
            99,
            30,40,50
        ]

        @test run_program!([1,0,0,0,99]).mem.vect == [2,0,0,0,99]
        @test run_program!([2,3,0,3,99]).mem.vect == [2,3,0,6,99]
        @test run_program!([2,4,4,5,99,0]).mem.vect == [2,4,4,5,99,9801]
        @test run_program!([1,1,1,4,99,5,6,0,99]).mem.vect == [30,1,1,4,2,5,6,0,99]
    end

    btest("2.1") do
        part1() == 2782414
    end

    btest("2.2") do
        part2(19690720) == 9820
    end
end


text:
function instantiate_function(f, x, ::AbstractADType, p, num_cons = 0)
    grad   = f.grad   === nothing ? nothing : (G,x)->f.grad(G,x,p)
    hess   = f.hess   === nothing ? nothing : (H,x)->f.hess(H,x,p)
    hv     = f.hv     === nothing ? nothing : (H,x,v)->f.hv(H,x,v,p)
    cons   = f.cons   === nothing ? nothing : (x)->f.cons(x,p)
    cons_j = f.cons_j === nothing ? nothing : (res,x)->f.cons_j(res,x,p)
    cons_h = f.cons_h === nothing ? nothing : (res,x)->f.cons_h(res,x,p)

    OptimizationFunction{true,SciMLBase.NoAD,typeof(f.f),typeof(grad),
                         typeof(hess),typeof(hv),typeof(cons),
                         typeof(cons_j),typeof(cons_h)}(f.f,
                         SciMLBase.NoAD(),grad,hess,hv,cons,
                         cons_j,cons_h)
end


text:
# Load in, combine, and plot chains:
using DelimitedFiles
using JLD2
using Statistics
using Printf
using PyPlot
include("nlog_prior.jl")
if !@isdefined(CGS)
  include("../../src/CGS.jl")
  using Main.CGS
end

n = 8
nplanet = n-1
nparam = (n-1)*5+2
nstep = 2000

#@load "T1_hmc_total_02132020.jld2"
@load "../../data/T1_hmc_total_02212020.jld2"

x_opt = [4.598791344049799e-5, 1.5108213441190745, 7257.550484415569, -0.00495015358026971, 0.004671261296998595,
         4.40588192304326e-5, 2.4219505838480635, 7258.587230337616, -0.001188890950465861, 0.0014051713506332184,
         1.305145694403463e-5, 4.049193371144702, 7257.06806273418, -0.0064708229088599216, 0.003529480254615964,
         2.3551037452505444e-5, 6.101020579519802, 7257.827698633612, 0.003217590595460454, -0.003935834596847923,
         3.4925310671129774e-5, 9.207561405004297, 7257.073658382135, -0.009581969253521241, -0.00016501113996023978,
         4.417497904782645e-5, 12.352405680378052, 7257.715102776943, 0.00288626733667452, 0.0016434273040851228,
         1.0931051156428014e-5, 18.772977044007867, 7249.605686691679, -0.004439910317349172, 0.0002702876058394804,
         1.2119892906238703, 0.8222963482043892]

planet = ["b","c","d","e","f","g","h"]
param = ["mu","P","t0","ecos","esin"]

# Now print mass ratios:
xfit = zeros(2,nplanet)
fac = 0.09*MSUN/MEARTH
#println("\$M_p/M_* (M_\\oplus/0.09 M_\\odot)\$ & \$ M_p/M_*\$ & \$\\sigma(M_p/M_*)/(M_p/M_*)\$ & \$ P\$ [day] & \$ t0\$ [JD - 2,450,000] & \$ e\\cos{\\omega}\$ & \$ e\\sin{\\omega}\$  \cr")
for i=1:nplanet
  planetstring = string(planet[i]," &")
  xfit[:,i] = [mean(state_total[(i-1)*5+1,:])*fac,std(state_total[(i-1)*5+1,:])*fac]
  global planetstring = string(planetstring," \$",@sprintf("%6.4f",xfit[1,i])," \\pm ",@sprintf("%6.4f",xfit[2,i]),"\$ &")
  xfit[:,i] = [mean(state_total[(i-1)*5+1,:])*1e5,std(state_total[(i-1)*5+1,:])*1e5]
  global planetstring = string(planetstring," \$",@sprintf("%6.3f",xfit[1,i])," \\pm ",@sprintf("%6.3f",xfit[2,i]),"\$ &")
  planetstring = string(planetstring," \$",@sprintf("%3.2f",xfit[2,i]/xfit[1,i]*100),"\$ &")
  xfit[:,i] = [mean(state_total[(i-1)*5+2,:]),std(state_total[(i-1)*5+2,:])]
  planetstring= string(planetstring," \$",@sprintf("%8.6f",xfit[1,i])," \\pm ",@sprintf("%8.6f",xfit[2,i]),"\$ &")
  xfit[:,i] = [mean(state_total[(i-1)*5+3,:]),std(state_total[(i-1)*5+3,:])]
  planetstring= string(planetstring," \$",@sprintf("%7.5f",xfit[1,i])," \\pm ",@sprintf("%7.5f",xfit[2,i]),"\$ &")
  xfit[:,i] = [mean(state_total[(i-1)*5+4,:]),std(state_total[(i-1)*5+4,:])]
  planetstring= string(planetstring," \$",@sprintf("%6.5f",xfit[1,i])," \\pm ",@sprintf("%6.5f",xfit[2,i]),"\$ &")
  xfit[:,i] = [mean(state_total[(i-1)*5+5,:]),std(state_total[(i-1)*5+5,:])]
  planetstring= string(planetstring," \$",@sprintf("%6.5f",xfit[1,i])," \\pm ",@sprintf("%6.5f",xfit[2,i]),"\$ \\cr ")
  println(planetstring)
end

println("\$ \\log{\\nu} = ",@sprintf("%6.4f",mean(state_total[36,:]))," \\pm ",@sprintf("%6.4f",std(state_total[36,:])),"\$")
println("\$ V_1 e^{1/(2\\nu)} = ",@sprintf("%6.4f",mean(state_total[37,:]))," \\pm ",@sprintf("%6.4f",std(state_total[37,:])),"\$")


text:
function random_labelled_graph()
    n = 15
    Random.seed!(1)
    A = Float64[ rand() < 0.5 ? 0 : rand() for i=1:n, j=1:n]
    for i=1:n
        A[i, 1:i-1] = A[1:i-1, i]
        A[i, i] = 0
    end
    x = rand(n)
    y = rand(n)
    z = rand(n)
    p = graphplot(A,
              nodesize = 0.2,
              node_weights = 1:n,
              nodecolor = range(colorant"yellow", stop=colorant"red", length=n),
              names = 1:n,
              fontsize = 10,
              linecolor = :darkgrey,
              layout_kw = Dict(:x => x, :y => y),
              )
    p, n, A, x, y, z
end

function random_3d_graph()
    n, A, x, y, z = random_labelled_graph()[2:end]
    graphplot(A,
              node_weights = 1:n,
              markercolor = :darkgray,
              dim = 3,
              markersize = 5,
              markershape = :circle,
              linecolor = :darkgrey,
              linealpha = 0.5,
              layout_kw = Dict(:x => x, :y => y, :z => z),
              )
end

function light_graphs()
    g = wheel_graph(10)
    graphplot(g, curves=false)
end

function directed()
    g = [0 1 1;
         0 0 1;
         0 1 0]

    graphplot(g, names=1:3, curvature_scalar=0.1)
end

function edgelabel()
    n = 8
    g = wheel_digraph(n)
    edgelabel_dict = Dict()
    for i in 1:n
        for j in 1:n
            edgelabel_dict[(i, j)] = string("edge ", i, " to ", j)
        end
    end

    graphplot(g, names=1:n, edgelabel=edgelabel_dict, curves=false, nodeshape=:rect)
end

function selfedges()
    g = [1 1 1;
         0 0 1;
         0 0 1]

    graphplot(DiGraph(g), self_edge_size=0.2)
end

function multigraphs()
    graphplot([[1,1,2,2],[1,1,1],[1]], names="node_".*string.(1:3), nodeshape=:circle, self_edge_size=0.25)
end

function arc_chord_diagrams()
    adjmat = Symmetric(sparse(rand(0:1,8,8)))
    plot(
         graphplot(adjmat,
                   method=:chorddiagram,
                   names=[text(string(i), 8) for i in 1:8],
                   linecolor=:black,
                   fillcolor=:lightgray),

         graphplot(adjmat,
                   method=:arcdiagram,
                   markersize=0.5,
                   markershape=:circle,
                   linecolor=:black,
                   markercolor=:black)
    )
end

function marker_properties()
    N = 8
    g = barabasi_albert(N, 1; seed = 42)
    weights = [length(neighbors(g, i)) for i in 1:nv(g)]
    Random.seed!(42)
    graphplot(g, curvature_scalar=0,
              node_weights=weights, nodesize=0.25,
              linecolor=:gray,
              linewidth=2.5,
              nodeshape=:circle,
              node_z=rand(N), markercolor=:viridis,
              nodestrokewidth=1.5,
              markerstrokestyle=:solid,
              markerstrokealpha=1.0,
              markerstrokecolor=:lightgray,
              colorbar=true,
    )
end

function ast_example()
    code = :(
    function mysum(list)
        out = 0
        for value in list
            out += value
        end
        out
    end
    )

    plot(code, fontsize=10, shorten=0.01, axis_buffer=0.15, nodeshape=:rect, size=(1000, 1000))
end

function julia_type_tree()
    plot(AbstractFloat, method=:tree, fontsize=10, nodeshape=:ellipse, size=(1000, 1000))
end

AbstractTrees.children(d::Dict) = [p for p in d]
AbstractTrees.children(p::Pair) = AbstractTrees.children(p[2])
function AbstractTrees.printnode(io::IO, p::Pair)
    str = isempty(AbstractTrees.children(p[2])) ? string(p[1], ": ", p[2]) : string(p[1], ": ")
    print(io, str)
end

function julia_dict_tree()
    d = Dict(:a => 2,:d => Dict(:b => 4,:c => "Hello"),:e => 5.0)

    plot(TreePlot(d), method=:tree, fontsize=10, nodeshape=:ellipse, size=(1000, 1000))
end

function diamond_nodeshape(x_i, y_i, s)
    [(x_i + 0.5s*dx, y_i + 0.5s*dy) for (dx, dy) in [(1,1),(-1,1),(-1,-1),(1,-1)]]
end
function diamond_nodeshape_wh(x_i, y_i, h, w)
    out = Tuple{Float64, Float64}[(-0.5,0),(0,-0.5),(0.5,0),(0,0.5)]
    map(out) do t
        x = t[1]* h
        y = t[2]* w
        (
         x + x_i,
         y + y_i
        )
    end
end

function custom_nodeshapes_single()
    Random.seed!(6)
    g = rand(5,5)
    g[g .> 0.5] .= 0
    for i in 1:5
        g[i,i] = 0
    end
    graphplot(g, nodeshape=diamond_nodeshape)
end

function custom_nodeshapes_various()
    Random.seed!(6)
    g = rand(5,5)
    g[g .> 0.5] .= 0
    for i in 1:5
        g[i,i] = 0
    end
    graphplot(g, nodeshape=[:circle, diamond_nodeshape, diamond_nodeshape_wh, :hexagon, diamond_nodeshape_wh])
end

function funky_edge_and_marker_args()
    Random.seed!(6)
    n = 5
    g = SimpleDiGraph(n)

    add_edge!(g, 1, 2)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 4)
    add_edge!(g, 4, 4)
    add_edge!(g, 4, 5)

    curviness_matrix = zeros(n,n)
    edgewidth_matrix = zeros(n,n)
    edgestyle_dict = Dict()
    for e in edges(g)
        curviness_matrix[e.src,e.dst] = 0.5sin(e.src)
        edgewidth_matrix[e.src,e.dst] = 0.8e.dst
        edgestyle_dict[(e.src,e.dst)] = e.src < 2.0 ? :solid :
                                        e.src > 3.0 ? :dash : :dot
    end
    edge_z_function = (s,d,w) -> curviness_matrix[s,d]

    graphplot(g, names=[" I "," am "," a ","funky","graph"],
        x=[1,2,3,4,5],
        y=[5,4,3,2,1],
        nodesize=0.3,
        size=(1000,1000),
        axis_buffer=0.6,
        fontsize=16,
        self_edge_size=1.3,
        curvature_scalar=curviness_matrix,
        edgestyle=edgestyle_dict,
        edgewidth=edgewidth_matrix,
        edge_z=edge_z_function,
        nodeshape=:circle,
        node_z=[1,2,3,4,5],
        nodestroke_z=[5,4,3,2,1],
        edgecolor=:viridis,
        markercolor=:viridis,
        nodestrokestyle=[:dash,:solid,:dot],
        nodestrokewidth=6,
        linewidth=2,
        colorbar=true,
    )
end


text:
nn = Chain(Dense(21, 21, tanh), Dense(21, 21))
nn1 = FastChain(FastDense(21, 21, tanh), FastDense(21, 21))

X = randn(Float32, 21, 10)
Y = rand(Float32, 21, 10)
sci_train!(nn, (X, Y), ADAM())
sci_train!(nn, Flux.Data.DataLoader((X, Y)), ADAM(); device = cpu, epoch = 1)
sci_train(nn1, (X, Y))

loss(p) = nn1(X, p) |> sum
p1 = initial_params(nn1)
sci_train(loss, p1)

cd(@__DIR__)
model = load_model("model.h5"; mode = :tf)
sci_train!(model, (randn(Float32, 1, 4), randn(Float32, 1, 1)))


text:
# Input data
data = map(x -> parse(Int64, x), split(readline("Day07.input"), ","));

# Part 1
mindistance =
    (F, data) ->
        reduce((acc, v) -> min(acc, v), map(x -> F(x, data), range(extrema(data)...)))

ans = mindistance((x, data) -> sum(map(v -> abs(v - x), data)), data);
println("Part 1: $ans")

# Part 2
ans = mindistance(
    (x, data) -> sum(map(v -> (n = abs(v - x); (n * (n + 1)) >> 1), data)),
    data,
);
println("Part 2: $ans")


text:
# ===========================================================================
# Spacecraft Rendezvous model
# See https://easychair.org/publications/paper/S2V for a description and
# https://gitlab.com/goranf/ARCH-COMP/tree/master/2018/AFF/SpaceEx/Rendezvous
# for a reference model
# ===========================================================================
using SparseArrays, HybridSystems, MathematicalSystems, LazySets, Reachability
using LazySets: HalfSpace  # resolve name-space conflicts with Polyhedra

"""
    spacecraft(; abort_time::Union{Float64, Vector{Float64}}=-1.)

Construct the spacecraft-rendezvous model.

### Input

- `abort_time` -- (optional, default: `-1.`) time of abort; can be either a
                  number (time point) or an interval (time range); an interval
                  is given as a vector `[l, u]` containing the lower (`l`) and
                  the upper (`u`) bound; a negative number is interpreted as a
                  *no-abort* scenario
"""
function spacecraft(; abort_time::Union{Float64, Vector{Float64}}=-1.)
    # variables
    x = 1  # x position (negative!)
    y = 2  # y position (negative!)
    vx = 3  # x velocity
    vy = 4  # y velocity
    t = 5   # time

    # number of variables
    n = 4 + 1

    # flag for activating the "rendezvous abort" scenario
    aborting = true
    if abort_time isa Number
        if abort_time < 0.
            aborting = false
        else
            t_abort_lower = abort_time
            t_abort_upper = abort_time
        end
    else
        @assert length(abort_time) == 2 "abort time must be a point or an interval"
        t_abort_lower = abort_time[1]
        t_abort_upper = abort_time[2]
    end

    # discrete structure (graph)
    automaton = LightAutomaton(aborting ? 3 : 2)

    # common vector of affine dynamics
    b = sparsevec([t], [1.], n)

    # mode 1 ("approaching")
    A = spzeros(n, n)
    A[x, vx] = 1.
    A[y, vy] = 1.
    A[vx, x] = -0.057599765881773
    A[vx, y] = 0.000200959896519766
    A[vx, vx] = -2.89995083970656
    A[vx, vy] = 0.00877200894463775
    A[vy, x] = -0.000174031357370456
    A[vy, y] = -0.0665123984901026
    A[vy, vx] = -0.00875351105536225
    A[vy, vy] = -2.90300269286856
    invariant = HalfSpace(sparsevec([x], [1.], n), -100.)  # x <= -100
    if aborting
        invariant = HPolyhedron([
            invariant,
            HalfSpace(sparsevec([t], [1.], n), t_abort_upper)  # t <= t_abort_upper
           ])
    end
    m_1 = CACS(A, b, invariant)

    # mode 2 ("rendezvous attempt")
    A = spzeros(n, n)
    A[x, vx] = 1.
    A[y, vy] = 1.
    A[vx, x] = -0.575999943070835
    A[vx, y] = 0.000262486079431672
    A[vx, vx] = -19.2299795908647
    A[vx, vy] = 0.00876275931760007
    A[vy, x] = -0.000262486080737868
    A[vy, y] = -0.575999940191886
    A[vy, vx] = -0.00876276068239993
    A[vy, vy] = -19.2299765959399
    invariant = HPolyhedron([
        HalfSpace(sparsevec([x], [-1.], n), 100.),           # x >= -100
        HalfSpace(sparsevec([x], [1.], n), 100.),            # x <= 100
        HalfSpace(sparsevec([y], [-1.], n), 100.),           # y >= -100
        HalfSpace(sparsevec([y], [1.], n), 100.),            # y <= 100
        HalfSpace(sparsevec([x, y], [-1., -1.], n), 141.1),  # x + y >= -141.1
        HalfSpace(sparsevec([x, y], [1., 1.], n), 141.1),    # x + y <= 141.1
        HalfSpace(sparsevec([x, y], [1., -1.], n), 141.1),   # -x + y >= -141.1
        HalfSpace(sparsevec([x, y], [-1., 1.], n), 141.1)    # -x + y <= 141.1
       ])
    if aborting
        addconstraint!(invariant, HalfSpace(sparsevec([t], [1.], n), t_abort_upper))  # t <= t_abort_upper
    end
    m_2 = CACS(A, b, invariant)

    # mode 3 ("aborting")
    A = spzeros(n, n)
    A[x, vx] = 1.
    A[y, vy] = 1.
    A[vx, x] = 0.0000575894721132
    A[vx, vy] = 0.00876276
    A[vy, vx] = -0.00876276
    invariant = Universe(n)
    m_3 = CACS(A, b, invariant)

    # modes
    modes = aborting ? [m_1, m_2, m_3] : [m_1, m_2]

    # transition 1 -> 2
    add_transition!(automaton, 1, 2, 1)
    guard = HPolyhedron([
        HalfSpace(sparsevec([x], [-1.], n), 100.),           # x >= -100
        HalfSpace(sparsevec([x], [1.], n), 100.),            # x <= 100
        HalfSpace(sparsevec([y], [-1.], n), 100.),           # y >= -100
        HalfSpace(sparsevec([y], [1.], n), 100.),            # y <= 100
        HalfSpace(sparsevec([x, y], [-1., -1.], n), 141.1),  # x + y >= -141.1
        HalfSpace(sparsevec([x, y], [1., 1.], n), 141.1),    # x + y <= 141.1
        HalfSpace(sparsevec([x, y], [1., -1.], n), 141.1),   # -x + y >= -141.1
        HalfSpace(sparsevec([x, y], [-1., 1.], n), 141.1)    # -x + y <= 141.1
       ])
    t1 = ConstrainedIdentityMap(n, guard)

    # TODO the SpaceEx model is missing the transition from 2 to 1

    if aborting
        add_transition!(automaton, 1, 3, 2)
        guard = HalfSpace(sparsevec([t], [-1.], n), -t_abort_lower)  # t >= t_abort_lower
        t2 = ConstrainedIdentityMap(n, guard)

        add_transition!(automaton, 2, 3, 3)
        guard = HalfSpace(sparsevec([t], [-1.], n), -t_abort_lower)  # t >= t_abort_lower
        t3 = ConstrainedIdentityMap(n, guard)
    end

    # transition annotations
    resetmaps = aborting ? [t1, t2, t3] : [t1]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode 1
    X0 = Hyperrectangle([-900., -400., 0., 0., 0.],
                        [25., 25., 0., 0., 0.])
    initial_condition = [(1, X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    # safety property in "Mode 2"
    velocity = 0.055 * 60.  # meters per minute
    cx = velocity * cos(π / 8)  # x-coordinate of the octagon's first (ENE) corner
    cy = velocity * sin(π / 8)  # y-coordinate of the octagon's first (ENE) corner
    octagon = [
        HalfSpace(sparsevec([vx], [-1.], n), cx),                # vx >= -cx
        HalfSpace(sparsevec([vx], [1.], n), cx),                 # vx <= cx
        HalfSpace(sparsevec([vy], [-1.], n), cx),                # vy >= -cx
        HalfSpace(sparsevec([vy], [1.], n), cx),                 # vy <= cx
        HalfSpace(sparsevec([vx, vy], [1., 1.], n), cy + cx),    # vx + vy <= cy + cx
        HalfSpace(sparsevec([vx, vy], [1., -1.], n), cy + cx),   # vx - vy <= cy + cx
        HalfSpace(sparsevec([vx, vy], [-1., 1.], n), cy + cx),   # -vx + vy <= cy + cx
        HalfSpace(sparsevec([vx, vy], [-1., -1.], n), cy + cx)   # -vx - vy <= cy + cx
       ]
    tan30 = tan(π/6)
    cone = [
        HalfSpace(sparsevec([x], [-1.], n), 100.),          # x >= -100
        HalfSpace(sparsevec([x, y], [tan30, -1.], n), 0.),  # -x tan(30°) + y >= 0
        HalfSpace(sparsevec([x, y], [tan30, 1.], n), 0.),   # -x tan(30°) - y >= 0
       ]
    property_rendezvous = SafeStatesProperty(HPolytope([octagon; cone]))
    # safety property in "aborting"
    target = HPolytope([
        HalfSpace(sparsevec([x], [-1.], n), 0.2),  # x >= -0.2
        HalfSpace(sparsevec([x], [1.], n), 0.2),   # x <= 0.2
        HalfSpace(sparsevec([y], [-1.], n), 0.2),  # y >= -0.2
        HalfSpace(sparsevec([y], [1.], n), 0.2),   # y <= 0.2
       ])
    property_aborting = BadStatesProperty(target)
    # safety properties
    property = Dict{Int, Property}(2 => property_rendezvous,
                                   3 => property_aborting)

    # default options
    options = Options(:T=>300., :property=>property)

    return (system, options)
end

function run_spacecraft(system, options)
    opC = BFFPSV18(:partition => [1:5], :δ => 0.04)
    opD = LazyDiscretePost(:lazy_R⋂I => true, :lazy_R⋂G => true)
    options[:mode] = "check"
    options[:plot_vars] = [1, 2]
    options[:project_reachset] = true

    solve(system, options, opC, opD)
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
function add_additional_equations(obj::AbstractPropagatorBuilder, arg0::AdditionalEquations)
    return jcall(obj, "addAdditionalEquations", void, (AdditionalEquations,), arg0)
end

function build_propagator(obj::PropagatorBuilder, arg0::Vector{jdouble})
    return jcall(obj, "buildPropagator", Propagator, (Vector{jdouble},), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_attitude_provider(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_frame(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getFrame", Frame, ())
end

function get_initial_orbit_date(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getInitialOrbitDate", AbsoluteDate, ())
end

function get_mu(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getMu", jdouble, ())
end

function get_orbit_type(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getOrbitType", OrbitType, ())
end

function get_orbital_parameters_drivers(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getOrbitalParametersDrivers", ParameterDriversList, ())
end

function get_position_angle(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getPositionAngle", PositionAngle, ())
end

function get_position_scale(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getPositionScale", jdouble, ())
end

function get_propagation_parameters_drivers(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getPropagationParametersDrivers", ParameterDriversList, ())
end

function get_selected_normalized_parameters(obj::AbstractPropagatorBuilder)
    return jcall(obj, "getSelectedNormalizedParameters", Vector{jdouble}, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function reset_orbit(obj::AbstractPropagatorBuilder, arg0::Orbit)
    return jcall(obj, "resetOrbit", void, (Orbit,), arg0)
end

function set_attitude_provider(obj::AbstractPropagatorBuilder, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
# =========================================================================== #
# Compliant julia 1.x

# Using the following packages
using JuMP, GLPK
using LinearAlgebra

include("loadSPP.jl")
include("setSPP.jl")
include("getfname.jl")
include("zDM1_1.jl")
include("DM1_2.jl")

# =========================================================================== #

# Setting the data
#fname = "B:/Cours/Nantes/Metaheuristique/DM1_Metaheuristique/solveSPP-master/Data/didactic.dat"  # path for a standard config on macOS
#cost, matrix, n, m = loadSPP(fname)

# Proceeding to the optimization
#solverSelected = GLPK.Optimizer
#ip, ip_x = setSPP(solverSelected, cost, matrix)
#println("Solving..."); optimize!(ip)

# Displaying the results
#println("z  = ", objective_value(ip))
#sprint("x  = "); println(value.(ip_x))

# =========================================================================== #

# Collecting the names of instances to solve
target = "B:/Cours/Nantes/Metaheuristique/DM1_Metaheuristique/solveSPP-master/Data"            # path for a standard config on macOS
fnames = getfname(target)

#===================================================#
for f in fnames
    println("=================")
    println(f, "\n")
    cost, matrix, n, m = loadSPP(f)
    @time (SOL,crts,z) =Glouton(cost, matrix, n, m)
    @time (SOL2) = kpexchange!(
        SOL, # Notre solution
        1,  # le nombre d'objet à retirer du conteneur
        2, # le nombre d'objet à rajouter dans le conteneur
        n::Int,  # taille de x
        m::Int, # nombre de contraintes
        z,  # solution trouvée
        z, #solution actuelle
        cost, # tableau des coûts de x
        SOL, #meilleur solution trouvée
        crts, # vecteur des contraintes (on suppose que pout tout i crts[i] =< 1)
        matrix # la matrice des contraintes
        )
        println("Solution après amélioration: ", calculz(SOL2[2],cost,n))
        #println("Solution après amelioration: ", SOL)
    #println(f,"\n Solution: ",SOL)
    #println("\n actif: ",a2)
end
cd("../")

# ============================================================================ #
#n = 9
#m = 7
#a  = zeros(Float64, n)
#cost = [10, 5, 8, 6, 9, 13, 11, 4, 6]
#matrix = [1 1 1 0 1 0 1 1 0; 0 1 1 0 0 0 0 1 0;0 1 0 0 1 1 0 1 1;0 0 0 1 0 0 0 0 0;1 0 1 0 1 1 0 0 1;0 1 1 0 0 0 1 0 1; 1 0 0 1 1 0 0 1 1]
#a = Glouton(cost, matrix, n, m)
#println("Solution: ",a)

max_stars_repo_path: Archives/solveSPP- TOULA - ROCHET/main.jl


text:
module Visualize

export visualize,
       external_visualize

using AbstractPlotting
using ..MDTrajAnalyse

const colours = Dict(
    "C"=>:black,
    " C"=>:black,
    "O"=>:red,
    " O"=>:red,
    "H"=>:white,
    " H"=>:white,
    "Ar"=>:orange,
    "N"=>:blue,
    " N"=>:blue
    )

const sizes = Dict(
    "C"=>0.9,
    " C"=>0.9,
    "O"=>0.9,
    " O"=>0.9,
    "H"=>0.6,
    " H"=>0.6,
    "Ar"=>1.0,
    "N"=>0.9,
    " N"=>0.9
    )

function visualize(traj::AbstractTrajectory, frame)
    acol = map( x -> colours[x], traj.names)
    asizes = map( x -> sizes[x], traj.names)
    sc = meshscatter(traj[frame]; markersize=asizes,color=acol)
    return sc
end

function visualize(traj::AbstractTrajectory; step=100)
    acol = map( x -> colours[x], traj.names)
    asizes = map( x -> sizes[x], traj.names)
    fig = Figure()
    ax = Axis3(fig[1,1])
    i = Slider(fig[2,1]; range=1:step:length(traj), startvalue=1)
    x = lift( j -> traj[j][1,:], i.value)
    y = lift( j -> traj[j][2,:], i.value)
    z = lift( j -> traj[j][3,:], i.value)
    meshscatter!(ax,x,y,z; markersize=asizes, color=acol)
    return fig
end

function external_visualize(traj::AbstractTrajectory;
    cmd="vlc", stdout=devnull, stderr=devnull, command="vmd"
    )
    fname= tempname()*".xyz"
    write_xyz(fname, traj)
    run(pipeline(`$(command) $(fname)`, stdout=stdout, stderr=stderr))
    rm(fname)
end

end # module


text:
function error_analysis(sys, e)
    if e isa ModelingToolkit.ExtraVariablesSystemException
        eqs = sys.eqs
        states = sys.states
        t = sys.iv

        rs = [Differential(t)(state) => state for state in states]
        extrastates = [state for state in states]
        extraeqs = [eq for eq in eqs]
        for r in rs
            for eq in extraeqs
                if subsmatch(eq.lhs, r) | subsmatch(eq.rhs, r)
                    setdiff!(extrastates, [r.second])
                    setdiff!(extraeqs, [eq])
                    break
                end
            end
        end
        println("The system of equations is:")
        println(eqs)
        println()
        println("There are $(length(states)) variables and $(length(eqs)) equations.\n")
        println("The variables without time derivatives are:")
        println(extrastates)
        println()
        println("The equations without time derivatives are:")
        println(extraeqs)
        rethrow(e)
    else
        rethrow(e)
    end
end


text:
using Reexport

@reexport module SamplingConstraints
    import Distributions
    # Sampling constraints types
    include("constraint_definitions.jl")

    # Fallbacks when encountering incompatible sampling constraints
    include("fallback_constraints.jl")

    # Truncate values
    include("truncation/truncate.jl")

    # Constrain uncertain values
    include("constrain_uncertainvalue.jl")
    include("constrain_certain_value.jl")
    include("constrain_population.jl")

    # Constrain uncertain datasets
    include("constrain_uncertaindataset.jl")
    include("constrain_uncertainvaluedataset.jl")
    include("constrain_uncertainindexdataset.jl")

    include("ordered_sequences/ordered_sequences.jl")

end # module

"""
    SamplingConstraints

A module defining:
1. Sampling constraints for the data types in the `UncertainValues` and
    `UncertainDatasets` modules.
2. Functions for resampling those data types.
"""
SamplingConstraints


text:
## Activate project environment. NB: This assumes that you are running Julia
## from the project root. If not, you'll need to adjust the activation path to
## the underlying .toml files accordingly. E.g. Use `Pkg.activate("..")` if you
## are running this script directly from the `/julia` sub-directory.
using Pkg
Pkg.activate(".")

## Packages and functions
using MimiPAGE2009
using Mimi
import Mimi: EmpiricalDistribution
using RCall
using DataFrames
using Distributions
using CSVFiles

# Example of how to get 10 draws of the SCC (in 2020 dollars) using the
## MimiPAGE2009 model defaults:
# scc = MimiPAGE2009.compute_scc(year = 2020, n = 10, seed = 1234)

## Extract the empirical TCR distributions from my Bayesian regressions. I'm
## using RCall directly here, because I had trouble with the FstFileFormat.jl
## package.
R"""
library(here)
library(fst)
proj_dir = here()
tcr = read_fst(here('results/main/tcr.fst'))
"""
proj_dir = @rget(proj_dir)
tcr = @rget(tcr)

## Get the list of priors to loop over
priors = unique(tcr.prior)

## Create empty `scc` array (not preallocating types and size, so not the most
## efficient but doesn't make much different here)
scc = []

## Calculate the SCC for each prior type. I'm using 2020 prices and exporting
## the results for each loop to the `results/scc` dir.
for p in priors
    values = tcr[tcr[:prior].== p, :tcr]
    N = length(values)
    probs = repeat([1/N],  N)

    ## We need to redefine the MimiPAGE2009.getsim() function so that it includes
    ## the line `trc_transientresponse = EmpiricalDistribution(values, probs)`
    ## and thus references against my already-calculated values of TCR, rather
    ## than using the PAGE2009 defaults. See:
    ## https://github.com/anthofflab/MimiPAGE2009.jl/issues/222
    function MimiPAGE2009.getsim()
        mcs = @defsim begin

                ############################################################################
                # Define random variables (RVs)
                ############################################################################

                #The folllowing RVs are in more than one component.  For clarity they are
                #set here as opposed to below within the blocks of RVs separated by component
                #so that they are not set more than once.

                save_savingsrate = TriangularDist(10, 20, 15) # components: MarketDamages, NonMarketDamages, GDP, SLRDamages
                tcal_CalibrationTemp = TriangularDist(2.5, 3.5, 3.) # components: MarketDamages, NonMarketDamages

                wincf_weightsfactor["USA"] = TriangularDist(.6, 1, .8) # components: MarketDamages, NonMarketDamages, SLRDamages, Discountinuity
                wincf_weightsfactor["OECD"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor["USSR"] = TriangularDist(.2, .6, .4)
                wincf_weightsfactor["China"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor["SEAsia"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor["Africa"] = TriangularDist(.4, .8, .6)
                wincf_weightsfactor["LatAmerica"] = TriangularDist(.4, .8, .6)

                automult_autonomouschange = TriangularDist(0.5, 0.8, 0.65)  #components: AdaptationCosts, AbatementCosts

                #The following RVs are divided into blocks by component

                # CO2cycle
                air_CO2fractioninatm = TriangularDist(57, 67, 62)
                res_CO2atmlifetime = TriangularDist(50, 100, 70)
                ccf_CO2feedback = TriangularDist(4, 15, 10)
                ccfmax_maxCO2feedback = TriangularDist(30, 80, 50)
                stay_fractionCO2emissionsinatm = TriangularDist(0.25,0.35,0.3)

                # SulphateForcing
                d_sulphateforcingbase = TriangularDist(-0.8, -0.2, -0.4)
                ind_slopeSEforcing_indirect = TriangularDist(-0.8, 0, -0.4)

                # ClimateTemperature
                rlo_ratiolandocean = TriangularDist(1.2, 1.6, 1.4)
                pole_polardifference = TriangularDist(1, 2, 1.5)
                frt_warminghalflife = TriangularDist(10, 65, 30)
                #tcr_transientresponse = TriangularDist(1, 2.8, 1.3)
                #tcr_transientresponse = TriangularDist(5, 8, 6)
                tcr_transientresponse = EmpiricalDistribution(values, probs)

                # SeaLevelRise
                s0_initialSL = TriangularDist(0.1, 0.2, 0.15)
                sltemp_SLtemprise = TriangularDist(0.7, 3., 1.5)
                sla_SLbaselinerise = TriangularDist(0.5, 1.5, 1.)
                sltau_SLresponsetime = TriangularDist(500, 1500, 1000)

                # GDP
                isat0_initialimpactfxnsaturation = TriangularDist(20, 50, 30)

                # MarketDamages
                iben_MarketInitialBenefit = TriangularDist(0, .3, .1)
                W_MarketImpactsatCalibrationTemp = TriangularDist(.2, .8, .5)
                pow_MarketImpactExponent = TriangularDist(1.5, 3, 2)
                ipow_MarketIncomeFxnExponent = TriangularDist(-.3, 0, -.1)

                # NonMarketDamages
                iben_NonMarketInitialBenefit = TriangularDist(0, .2, .05)
                w_NonImpactsatCalibrationTemp = TriangularDist(.1, 1, .5)
                pow_NonMarketExponent = TriangularDist(1.5, 3, 2)
                ipow_NonMarketIncomeFxnExponent = TriangularDist(-.2, .2, 0)

                # SLRDamages
                scal_calibrationSLR = TriangularDist(0.45, 0.55, .5)
                #iben_SLRInitialBenefit = TriangularDist(0, 0, 0) # only usable if lb <> ub
                W_SatCalibrationSLR = TriangularDist(.5, 1.5, 1)
                pow_SLRImpactFxnExponent = TriangularDist(.5, 1, .7)
                ipow_SLRIncomeFxnExponent = TriangularDist(-.4, -.2, -.3)

                # Discountinuity
                rand_discontinuity = Uniform(0, 1)
                tdis_tolerabilitydisc = TriangularDist(2, 4, 3)
                pdis_probability = TriangularDist(10, 30, 20)
                wdis_gdplostdisc = TriangularDist(5, 25, 15)
                ipow_incomeexponent = TriangularDist(-.3, 0, -.1)
                distau_discontinuityexponent = TriangularDist(20, 200, 50)

                # EquityWeighting
                civvalue_civilizationvalue = TriangularDist(1e10, 1e11, 5e10)
                ptp_timepreference = TriangularDist(0.1,2,1)
                emuc_utilityconvexity = TriangularDist(0.5,2,1)

                # AbatementCosts
                AbatementCostParametersCO2_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,75,0)
                AbatementCostParametersCH4_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-25,100,0)
                AbatementCostParametersN2O_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,50,0)
                AbatementCostParametersLin_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,50,0)

                AbatementCostParametersCO2_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,40,20)
                AbatementCostParametersCH4_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
                AbatementCostParametersN2O_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
                AbatementCostParametersLin_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)

                AbatementCostParametersCO2_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)
                AbatementCostParametersCH4_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-8000,-1000,-4000)
                AbatementCostParametersN2O_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-15000,0,-7000)
                AbatementCostParametersLin_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)

                AbatementCostParametersCO2_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)
                AbatementCostParametersCH4_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
                AbatementCostParametersN2O_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
                AbatementCostParametersLin_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)

                AbatementCostParametersCO2_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,700,400)
                AbatementCostParametersCH4_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(3000,10000,6000)
                AbatementCostParametersN2O_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(2000,60000,20000)
                AbatementCostParametersLin_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,600,300)

                AbatementCostParametersCO2_ies_InitialExperienceStockofCutbacks = TriangularDist(100000,200000,150000)
                AbatementCostParametersCH4_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)
                AbatementCostParametersN2O_ies_InitialExperienceStockofCutbacks = TriangularDist(30,80,50)
                AbatementCostParametersLin_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)

                #the following variables need to be set, but set the same in all 4 abatement cost components
                #note that for these regional variables, the first region is the focus region (EU), which is set in the preceding code, and so is always one for these variables

                emitf_uncertaintyinBAUemissfactor["USA"] = TriangularDist(0.8,1.2,1.0)
                emitf_uncertaintyinBAUemissfactor["OECD"] = TriangularDist(0.8,1.2,1.0)
                emitf_uncertaintyinBAUemissfactor["USSR"] = TriangularDist(0.65,1.35,1.0)
                emitf_uncertaintyinBAUemissfactor["China"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["SEAsia"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["Africa"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["LatAmerica"] = TriangularDist(0.5,1.5,1.0)

                q0f_negativecostpercentagefactor["USA"] = TriangularDist(0.75,1.5,1.0)
                q0f_negativecostpercentagefactor["OECD"] = TriangularDist(0.75,1.25,1.0)
                q0f_negativecostpercentagefactor["USSR"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["China"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["SEAsia"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["Africa"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)

                cmaxf_maxcostfactor["USA"] = TriangularDist(0.8,1.2,1.0)
                cmaxf_maxcostfactor["OECD"] = TriangularDist(1.0,1.5,1.2)
                cmaxf_maxcostfactor["USSR"] = TriangularDist(0.4,1.0,0.7)
                cmaxf_maxcostfactor["China"] = TriangularDist(0.8,1.2,1.0)
                cmaxf_maxcostfactor["SEAsia"] = TriangularDist(1,1.5,1.2)
                cmaxf_maxcostfactor["Africa"] = TriangularDist(1,1.5,1.2)
                cmaxf_maxcostfactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)

                q0propmult_cutbacksatnegativecostinfinalyear = TriangularDist(0.3,1.2,0.7)
                qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = TriangularDist(1,1.5,1.3)
                c0mult_mostnegativecostinfinalyear = TriangularDist(0.5,1.2,0.8)
                curve_below_curvatureofMACcurvebelowzerocost = TriangularDist(0.25,0.8,0.45)
                curve_above_curvatureofMACcurveabovezerocost = TriangularDist(0.1,0.7,0.4)
                cross_experiencecrossoverratio = TriangularDist(0.1,0.3,0.2)
                learn_learningrate = TriangularDist(0.05,0.35,0.2)

                # AdaptationCosts
                AdaptiveCostsSeaLevel_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
                AdaptiveCostsSeaLevel_ci_costimpact_eu = TriangularDist(0.0005, 0.002, 0.001)
                AdaptiveCostsEconomic_cp_costplateau_eu = TriangularDist(0.005, 0.02, 0.01)
                AdaptiveCostsEconomic_ci_costimpact_eu = TriangularDist(0.001, 0.008, 0.003)
                AdaptiveCostsNonEconomic_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
                AdaptiveCostsNonEconomic_ci_costimpact_eu = TriangularDist(0.002, 0.01, 0.005)

                cf_costregional["USA"] = TriangularDist(0.6, 1, 0.8)
                cf_costregional["OECD"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["USSR"] = TriangularDist(0.2, 0.6, 0.4)
                cf_costregional["China"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["SEAsia"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["Africa"] = TriangularDist(0.4, 0.8, 0.6)
                cf_costregional["LatAmerica"] = TriangularDist(0.4, 0.8, 0.6)

                ############################################################################
                # Indicate which parameters to save for each model run
                ############################################################################

                save(EquityWeighting.td_totaldiscountedimpacts,
                        EquityWeighting.tpc_totalaggregatedcosts,
                        EquityWeighting.tac_totaladaptationcosts,
                        EquityWeighting.te_totaleffect,
                        co2cycle.c_CO2concentration,
                        TotalForcing.ft_totalforcing,
                        ClimateTemperature.rt_g_globaltemperature,
                        SeaLevelRise.s_sealevel,
                        SLRDamages.rgdp_per_cap_SLRRemainGDP,
                        MarketDamages.rgdp_per_cap_MarketRemainGDP,
                        NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP,
                        Discontinuity.rgdp_per_cap_NonMarketRemainGDP)

                end #def
        return mcs
    end

    ## Main function: Calculate the SCC for this priors's TCR distribution. I'll
    ## go ahead and put into a data frame for convenience.
    scc_p = DataFrame(scc = MimiPAGE2009.compute_scc(year = 2020, n = N, seed = 1234))
    scc_p[:prior] = p

    ## Push DF to parent SCC array
    push!(scc, scc_p)
end

## Concatenate DFs by row
scc = reduce(vcat, scc)

## Write to disk
save(joinpath(proj_dir, "results/scc/scc.csv"), scc)

max_stars_repo_name: grantmcdermott/skeptic-priors

max_issues_repo_name: grantmcdermott/skeptic-priors

max_forks_repo_name: grantmcdermott/sceptic-priors


text:
"A Solver holds the parameters of the differential equation solver."
struct Solver{T,ALG<:Union{OrdinaryDiffEqAlgorithm,Nothing},DT<:Union{T,Nothing},SP<:AbstractSpace}
    tspan::Tuple{T,T}
    algorithm::ALG
    dt::DT
    saveat::DT
    space_saving::AbstractSubsampler
    stiffness::Symbol
    #dense::Bool
    space::SP
end

function Solver{S}(; start_time::S=0.0, stop_time::S, dt::DT, time_save_every=1, space_save_every=1, algorithm::ALG=nothing, stiffness=:auto) where {S, ALG, DT<:Union{S,Nothing}}
    Solver{S,ALG,DT}((start_time, stop_time), algorithm, dt, time_save_every, space_save_every, stiffness)
end
function Solver{T,ALG,DT}(tspan::Tuple{T,T}, algorithm::ALG, simulated_dt::DT,
                          time_save_every::Int, space_save_every::Int, stiffness::Symbol) where {T,ALG,DT<:Union{S,Nothing}}
    Solver{T,ALG,DT}(tspan, algorithm, simulated_dt, time_save_every, [space_save_every], stiffness)
end

# """
#     saved_dt(solver)
#
# Return the time step saved by a solver.
# """
# saved_dt(s::Solver{T}) where T = s.simulated_dt * s.time_save_every
# """
#     initial_value(solver)
#     initial_value(simulation)
#
# Return the model's initial value (defaults to all zeros)
# """
# initial_value(solver::AbstractSolver) = error("undefined.")
# """
#     time_span(solver)
#     time_span(simulation)
#
# Return the time span over which the solver runs.
# """
# time_span(solver::Solver{T}) where T = solver.tspan
# """
#     save_idxs(solver)
#
# Return the indices into space of the values that the solver saves.
# """
# function save_idxs(solver::Solver{T}) where {T}#::Unio{P,n}{Nothing,Array{CartesianIndex}}
#     if all(solver.space_save_every .== 1)
#         return nothing
#     end
#     all_indices = CartesianIndices(initial_value(solver))
#     space_saved_subsample(all_indices, solver)
# end

# """
#     coordinates(model, solver)
#     coordinates(simulation)
#
# Return the spatial coordinates of values saved by `solver`
# """
# function saved_coordinates(solver::Solver)
#     @warn "not subsampling in _coordinates_"
#     coordinates(solver.space)
# #    collect(arr)[[StrideToEnd(i) for i in solver.space_save_every]...]
# end

# """
#     timepoints(solver)
#     timepoints(simulation)
#
# Return the times saved by `solver`.
# """
# function timepoints(solver::Solver{T}) where T
#     start, stop = time_span(solver)
#     start:saved_dt(solver):stop
# end
# saved_dx(solver::Solver)= step(solver.space) .* solver.space_save_every
#
# """
#     space_index_info(model, solver)
#     space_index_info(simulation)
#
# Return IndexInfo for saved space array.
# """
# space_index_info(solver::Solver{T}) where T = IndexInfo(saved_dx(solver), origin_idx(solver))
#
# """
#     time_index_info(solver)
#     time_index_info(simulation)
#
# Return IndexInfo for saved time array.
# """
# time_index_info(solver::Solver{T}) where T = IndexInfo(saved_dt(solver), (1,))
#
# function subsampling_idxs(solver::Solver, space_subsampler, time_subsampler)
#     x_info = get_space_index_info(solver)
#     t_info = get_time_index_info(solver)
#
#     x_dxs = subsampling_idxs(x_info, space_subsampler)
#     t_dxs = subsampling_idxs(t_info, time_subsampler)
#
#     return (x_dxs, 1, t_dxs)
# end

# FIXME bad assumptions
# function subsampling_idxs(simulation::Simulation{T,<:AbstractModel{T}}, x_target::AbstractArray, t_target::AbstractArray) where T
#     return (subsampling_space_idxs(simulation.model, simulation.solver, x_target)...,
#             1,
#             subsampling_time_idxs(simulation.solver, t_target))
# end
# function subsampling_time_idxs(solver::Solver, t_target::AbstractArray)
#     t_solver = time_span(solver)[1]:saved_dt(solver):time_span(solver)[end]
#     subsampling_idxs(t_target, t_solver)
# end
# function subsampling_space_idxs(solver::Solver, x_target::AbstractArray)
#     x_model = coordinates(model, solver)
#     subsampling_idxs(x_target, x_model)
# end

function _solve(problem::ODEProblem, solver::Solver{T,Euler}) where {T}
    # TODO: Calculate save_idxs ACCOUNTING FOR pops
    @show "Solving Euler"
    solve(problem, Euler(), dt=solver.simulated_dt,
            #saveat=saved_dt(solver),
            timeseries_steps=solver.time_save_every,
            save_idxs=save_idxs(solver))
end
function _solve(problem::ODEProblem, solver::Solver{T,Nothing}) where {T}
    @show "Solving with default ALG"
    solve(problem, saveat=saved_dt(solver), timeseries_steps=solver.time_save_every,
        save_idxs=save_idxs(solver), alg_hints=[solver.stiffness])
end


text:
module Example

using Glimmer, Glimmer.FlexUI
using DataFrames, CSV

println("Start [$(splitext(basename(@__FILE__))[1])]")

# ---------------------------------------------------------------
# define the application and some basic properties such as title and initial window size
# ---------------------------------------------------------------
app = App()
prop!(app, :title, "Glimmer Example - Grid Simple Table")
prop!(app, :winInitWidth, 1200)
prop!(app, :winInitHeight, 800)

# ---------------------------------------------------------------
# load the data
# ---------------------------------------------------------------
filename = joinpath(dirname(pathof(Glimmer)), "Data", "mlb_players.csv");
csv = CSV.File(filename)
table = DataFrame(csv)

# ---------------------------------------------------------------
# Define Variables
# ---------------------------------------------------------------
table_var = addVariable!(app, Variable(name="table_var", type="aggrid", value=""))

# setup the grid data and options - will get friendlier in the near future
table_var[] = GridUtils.table2agGrid(table)

# setting some grid options - because of it's complexity, the grid's data and visualization properties
# are all stored in the variable and not the UI control.
gridOption!(table_var, "rowSelection", "single")

# insert a row field
insertGridColumn!(table_var, 1, gridRowIndexCell(
    header="Row",
    width=70,
    cellStyle = Dict(
        "color" => "yellow",
        "background-color" => "pink",
    ),
))

# insert a combobox field in addition to the non-combobox Team field. you can edit either.
insertGridColumn!(table_var, 4, gridSelectCell(
    header="Teams Combobox",
    field="Team",
    selectOptions=[
        Dict(:key=>"BAL", :value=>"Baltimore Orioles"),
        Dict(:key=>"CWS", :value=>"College World Series"),
        Dict(:key=>"ANA", :value=>"Los Angeles Angels"),
        Dict(:key=>"BOS", :value=>"Boston Red Sox"),
        Dict(:key=>"CLE", :value=>"Cleveland Indians"),
        Dict(:key=>"OAK", :value=>"Oakland Athletics"),
        Dict(:key=>"NYY", :value=>"New York Yankees"),
        Dict(:key=>"DET", :value=>"Detroit Tigers"),
        Dict(:key=>"SEA", :value=>"Seattle Mariners"),
        Dict(:key=>"TB", :value=>"Tampa Bay Rays"),
        Dict(:key=>"KC", :value=>"Kansas City Royals"),
        Dict(:key=>"TEX", :value=>"Texas Rangers"),

    ]
))

# insert a button column and defining the event for it. further down we write a function to handle that event per row.
insertGridColumn!(table_var, 5, gridButtonCell(
    header="Name in REPL",
    eventName="myButtonClicked",
))

# define the column default values (note the cellStyle to color a line red according to documantation existance)
gridDefaultColDef!(table_var, gridCell(
    editable=true,
    resizable=true,
    # filter = "agTextColumnFilter",
    # floatingFilter = true,
    # cellStyle="__js__ params => params.data.hasDocs==false ? { color: 'red' } : { color: 'green' }",
))

# ---------------------------------------------------------------
# Define Events
# ---------------------------------------------------------------
# see https://www.ag-grid.com/angular-data-grid/grid-events/ for details

# this is an example of how to convert the grid variable back into a table - not used in this example
on(table_var, :valueChanged) do val
    @show "ON valueChanged"
    # new_table = table_var[];
    # @show new_table
end

on(table_var, :myButtonClicked) do val
    # @show "ON myButtonClicked"
    @info """Button Clicked: name=[$(val["data"]["Name"])]"""
end

on(table_var, :cellClicked) do val
    @show "ON cellClicked"
end

on(table_var, :cellDoubleClicked) do val
    @show "ON cellDoubleClicked"
end

on(table_var, :rowSelected) do val
    @show "ON rowSelected"
end

# when selection changed we update the documentation display for the new selected type
on(table_var, :selectionChanged) do val
    @show "ON selectionChanged"
end

# example - will use in the future
on(table_var, :myGridButtonEvent) do val
    @show "ON myGridButtonEvent" val
end

# ---------------------------------------------------------------
# Define Controls
# ---------------------------------------------------------------
ui = VContainer(
    Card(
        title="This Example demonstrate some the the AGGrid component capabilities",
        subtitle="This is a simple grid showing data loader from a CSV file. Also show an example of a combobox and a user defined button. Grid allow filtering.",
        content=VContainer(
            AGGrid(
                style="width: 100%; height: 500px;",
                variable="table_var",
            ),
        ),
    ),

    Glimmer.exampleSourceAsCard(@__FILE__),     # add the source code of the example as the last control
)
# set the controls for the application
controls!(app, ui)

# ---------------------------------------------------------------
# Run the application
# ---------------------------------------------------------------
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

end # module


text:
using Base64

"""
Storage driver which keeps the data embedded within the TOML file itself.
Useful for small amounts of self-contained data.

## Metadata spec

For Blob:
```
    [datasets.storage]
    driver="TomlDataStorage"
    type="Blob"
    data=\$(base64encode(data))
```

For BlobTree:
```
    [datasets.storage]
    driver="TomlDataStorage"
    type="BlobTree"

        [datasets.storage.data.\$(dirname1)]
        "\$(filename1)" = \$(base64encode(data1))
        "\$(filename2)" = \$(base64encode(data2))

        [datasets.storage.data.\$(dirname2)]
        ...
```
"""
struct TomlDataStorage
    dataset::DataSet
    data::Union{String,Dict{String,Any}}
end

# Get TOML data at `path`, returning nothing if not present
function _getpath(storage::TomlDataStorage, path::RelPath)
    x = storage.data
    for c in path.components
        x = get(x, c, nothing)
        !isnothing(x) || return nothing
    end
    x
end

#--------------------------------------------------
# Storage data interface for trees

Base.isdir(storage::TomlDataStorage, path::RelPath) = _getpath(storage, path) isa Dict
Base.isfile(storage::TomlDataStorage, path::RelPath) = _getpath(storage, path) isa String
Base.ispath(storage::TomlDataStorage, path::RelPath) = !isnothing(_getpath(storage, path))

Base.summary(io::IO, storage::TomlDataStorage) = print(io, "Data.toml")

function Base.readdir(storage::TomlDataStorage, path::RelPath)
    try
        tree = _getpath(storage, path)
        !isnothing(tree) || KeyError(path)
        sort!(collect(keys(tree::AbstractDict)))
    catch
        error("TOML storage requires trees to be as TOML dictionaries")
    end
end

#--------------------------------------------------
# Storage data interface for Blob

function Base.open(func::Function, as_type::Type{IO},
                   storage::TomlDataStorage, path; kws...)
    @context func(@! open(as_type, storage, path; kws...))
end

@! function Base.open(::Type{Vector{UInt8}}, storage::TomlDataStorage, path;
                      write=false, read=!write, kws...)
    if write
        error("Embedded data is read-only from within the DataSets interface")
    end
    try
        str = _getpath(storage, path)
        !isnothing(str) || KeyError(path)
        base64decode(str::AbstractString)
    catch
        error("TOML storage requires data to be as base64 encoded strings")
    end
end

@! function Base.open(::Type{IO}, storage::TomlDataStorage, path; kws...)
    buf = @! open(Vector{UInt8}, storage, path; kws...)
    IOBuffer(buf)
end

# TODO: The following should be factored out and implemented generically
function Base.read(storage::TomlDataStorage, path::RelPath, ::Type{T}) where {T}
    @context begin
        io = @! open(IO, storage, path)
        read(io, T)
    end
end

function Base.read(storage::TomlDataStorage, path::RelPath)
    @context @! open(Vector{UInt8}, storage, path)
end

#-------------------------------------------------------------------------------
# Connect storage backend
function connect_toml_data_storage(f, config, dataset)
    type = config["type"]
    data = get(config, "data", nothing)
    if type == "Blob"
        if !(data isa AbstractString)
            error("TOML data storage requires string data in the \"storage.data\" key")
        end
        f(Blob(TomlDataStorage(dataset, data)))
    elseif type == "BlobTree"
        if !(data isa AbstractDict)
            error("TOML data storage requires a dictionary in the \"storage.data\" key")
        end
        f(BlobTree(TomlDataStorage(dataset, data)))
    else
        throw(ArgumentError("DataSet type $type not supported for data embedded in Data.toml"))
    end
end

add_storage_driver("TomlDataStorage"=>connect_toml_data_storage)


text:
function get_glclim_stst()

    simid = :SimD3
    params = lglob(Dyn, simid, :params, :finite_cg)
    @extract params: Ds ϵs cg

    validator = Dict()
    for (ϵi, ϵ) in enumerate(ϵs)
       for D in Ds

          ## ------------------------------------------------------
          # Dynamic
          simparams = (;D, ϵ, cg)
          validator[simparams] = false
          status = get_status(simid, simparams)
          (status != STST_SIM_STATUS) && continue
          datfile = Dyn.simdat_file(simid, simparams)
          !isfile(datfile) && continue
          S = ldat(datfile)
          dyn_ug = S.ug_av
          dyn_X = S.X

          ## ------------------------------------------------------
          # LP
          net = S.V.net
          Dyn.fix!(net, :z, D)
          Dyn.ub!(net, :ug, cg * D / S.X)
          Dyn.lb!(net, :ug, 0.0)
          Dyn.reducebox!(net)
          Lug, Uug = Dyn.fva(net, :ug)
          UX = cg * D / Lug

          #  glc limited
          abs(Uug - dyn_ug) / abs(Uug + dyn_ug) > 0.1 && continue

          validator[simparams] = true

       end
    end

    return validator
 end


text:
function fill_rect(ctx::GraphicsContext,x::Int=0,y::Int=0,width::Int=32,height::Int=32)
	rectangle(ctx,x,y,width,height)
end
function fill_rect(ctx::GraphicsContext,x::Int=0,y::Int=0,width::Int=32,height::Int=32;color::RGB)
	set_source_rgb(ctx,red(color),green(color),blue(color))
	rectangle(ctx,x,y,width,height)
end
function fill_rect(ctx::GraphicsContext,x::Int=0,y::Int=0,width::Int=32,height::Int=32;color::RGBA)
	set_source_rgba(ctx,red(color),green(color),blue(color),alpha(color))
	rectangle(ctx,x,y,width,height)
end
function fill_text(ctx::GraphicsContext,x::Int=0,y::Int=0,text::String;color::RGB,size::Int)
	set_font_size(ctx,size)
	set_source_rgb(ctx,red(color),green(color),blue(color))
	move_to(ctx,x,y)
	show_text(ctx,text)
end
function fill_matrix(ctx::GraphicsContext,m::Matrix{ARGB32})
	sur=CairoImageSurface(m)
	set_source_surface(ctx,sur)
end


text:
# evaluate.jl

# evaluate and _evaluate functions

for TM in tupleTMs
    # Evaluates the TM1{T,S} at a::T; the computation includes the remainder
    @eval function evaluate(tm::$TM{T,S}, a) where {T,S}
        @assert iscontained(a, tm)

        _order = get_order(tm)

        if $(TM) == TaylorModel1
            Δ = tm.rem
        else
            Δ = tm.rem * a^(_order + 1)
        end

        return tm.pol(a) + Δ
    end

    @eval (tm::$TM{T,S})(a) where {T,S} = evaluate(tm, a)

    @eval evaluate(tm::Vector{$TM{T,S}}, a) where {T,S} = evaluate.(tm, a)

    # _evaluate corresponds to composition: substitute tmf into tmg
    # It **does not** include the remainder
    @eval function _evaluate(tmg::$TM, tmf::$TM)
        _order = get_order(tmf)
        @assert _order == get_order(tmg)

        tmres = zero(tmg.pol[_order]) * tmf
        tmres = tmres + tmg.pol[_order]
        @inbounds for k = _order-1:-1:0
            tmres = tmres * tmf
            tmres = tmres + tmg.pol[k]
        end

        # Returned result does not include remainder of tmg
        return tmres
    end

    @eval (tm::$TM)(x::$TM) = _evaluate(tm, x)

end

# Substitute a TaylorModelN into a TM1; it **does not** include the remainder
function _evaluate(tmg::TaylorModel1{T,S}, tmf::TaylorModelN{N,T,S}) where{N,T,S}
    _order = get_order(tmf)
    @assert _order == get_order(tmg)

    tmres = TaylorModelN(zero(constant_term(tmg.pol)), _order, tmf.x0, tmf.dom)
    @inbounds for k = _order:-1:0
        tmres = tmres * tmf
        tmres = tmres + tmg.pol[k]
    end

    # Returned result does not include remainder of tmg
    return tmres
end

(tm::TaylorModel1)(x::TaylorModelN) = _evaluate(tm, x)

# Evaluates the TMN on an interval, or array with proper dimension;
# the computation includes the remainder
function evaluate(tm::TaylorModelN{N,T,S}, a::IntervalBox{N,S}) where {N,T,S}
    @assert iscontained(a, tm)
    _order = get_order(tm)

    Δ = tm.rem

    return tm.pol(a) + Δ
end

(tm::TaylorModelN{N,T,S})(a::IntervalBox{N,S}) where {N,T,S} = evaluate(tm, a)

function evaluate(tm::TaylorModelN{N,T,S}, a::Array{R,1}) where {N,T,S,R}
    @assert iscontained(a, tm)
    _order = get_order(tm)

    Δ = tm.rem

    return tm.pol(a) + Δ
end

(tm::TaylorModelN{N,T,S})(a::Array{R,1}) where {N,T,S,R} = evaluate(tm, a)

evaluate(tm::Vector{TaylorModelN{N,T,S}}, a::IntervalBox{N,S}) where {N,T,S} =
    IntervalBox( [ tm[i](a) for i in eachindex(tm) ] )

function evaluate(a::Taylor1{TaylorModelN{N,T,S}}, dx::T) where {N, T<:Number, S<:Number}
    @inbounds suma = a[end]*one(dx)
    @inbounds for k in a.order-1:-1:0
        suma = suma*dx + a[k]
    end
    suma
end


text:
using GPUContagion
using Test

@testset "GPUContagion.jl" begin
    # Write your tests here.
end


text:
export LRMat;

struct LRMat{T<:Number}
    # variables
    height::    Int
    width::     Int
    UMat::      Array{T,2}
    VMat::      Array{T,2}
    # global settings
    EPS::       Float64
    MAXRANK::   Int

    function LRMat(D,Eps,MaxRank)
        h = size(D,1);
        w = size(D,2);
        [U,S,V] = svdtrunc(D,Eps,MaxRank);
        new(h,w,U*sqrt(S),V*sqrt(S),Eps,MaxRank);
    end
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
mutable struct NearestNeighboursUCB <: AbstractPolicy
    t::Int64
    batch::Int64
    num_arms::Int64
    initial_batches::Int64
    data::BanditDataset
    Xa::Vector{Matrix{Float64}}
    ra::Vector{Matrix{Float64}}
    θ::Float64
    φ::Float64

    function NearestNeighboursUCB(
        d::Int64, num_arms::Int64, initial_batches::Int64; θ::Float64=1.0, φ::Float64=1.0
    )
        Xa = [Matrix{Float64}(undef, d, 0) for _ in 1:num_arms]
        ra = [Matrix{Float64}(undef, 1, 0) for _ in 1:num_arms]
        data = BanditDataset(d)
        return new(0, 0, num_arms, initial_batches, data, Xa, ra, θ, φ)
    end
end

function nearest_neighbours(x::AbstractVector, X::AbstractMatrix)
    if size(X, 1) != length(x)
        throw(DimensionMismatch("length of x should match size(X, 1)"))
    end
    n = size(X, 2)
    distances = zeros(n)
    @inbounds for i in 1:n
        distances[i] = sqrt(sum((X[:, i] .- x) .^ 2))
    end
    idx = sortperm(distances)
    dst = distances[idx]
    return idx, dst
end

function _choose_k(x::AbstractVector, X::AbstractMatrix, t::Int64, θ::Float64, φ::Float64)
    idx, dst = nearest_neighbours(x, X)
    U = zeros(length(idx))
    @inbounds for k in 1:length(idx)
        radius = dst[k]
        U[k] = sqrt(θ * log(t) / k) + φ * log(t) * radius
    end
    min_k = argmin(U)
    return idx[1:min_k], U[min_k]
end

function (pol::NearestNeighboursUCB)(X::AbstractMatrix)
    n = size(X, 2)
    actions = zeros(Int64, n)
    if pol.batch <= pol.initial_batches
        for i in 1:n
            actions[i] = (pol.t + i) % pol.num_arms + 1
        end
        return actions
    end

    ucb = zeros(pol.num_arms)
    for i in 1:n
        for a in 1:(pol.num_arms)
            idx, U = _choose_k(X[:, i], pol.Xa[a], pol.t + 1, pol.θ, pol.φ)
            ucb[a] = sum(pol.ra[a][idx]) / length(idx) + U
        end
        actions[i] = argmax(ucb)
    end
    return actions
end

function update!(
    pol::NearestNeighboursUCB,
    X::AbstractMatrix,
    a::AbstractVector{<:Int},
    r::AbstractMatrix,
)
    pol.t += size(X, 2)
    pol.batch += 1
    append_data!(pol.data, X, a, r)
    for i in unique(a)
        pol.Xa[i] = hcat(pol.Xa[i], X[:, a .== i])
        pol.ra[i] = hcat(pol.ra[i], r[:, a .== i])
    end
end


text:
function interpolate!(data::Dict{Symbol, Any}, params::Parameters)
  survey_df = surveys2df(data)
  interpolate!(survey_df, params)
end

function interpolate!(survey_df::DataFrame, params::Parameters)
  dtor = params.dtor
  knots = (12.0*survey_df[!, :depth],)
  intp_inclination = interpolate(knots, reverse(dtor*survey_df[!, :inclination]),
    Gridded(Linear()))
  intp_dip = interpolate(knots, reverse(dtor*(90.0 .- survey_df[!, :inclination])),
    Gridded(Linear()))
  intp_heading = interpolate(knots, reverse(dtor*survey_df[!, :heading]),
    Gridded(Linear()))
  intp_holeradius = interpolate(knots, reverse(survey_df[!, :diameter]/2.0),
    Gridded(Linear()))
  for i in 1:length(params.elements)
    params.elements[i].inclination =
      (intp_inclination(params.nodes[i].z...) + intp_inclination(params.nodes[i+1].z...))/2.0
      params.elements[i].heading =
        (intp_heading(params.nodes[i].z...)+intp_heading(params.nodes[i+1].z...))/2.0
  end
  for i in 1:length(params.nodes)
    params.nodes[i].inclination = intp_inclination(params.nodes[i].z...)
    params.nodes[i].heading = intp_heading(params.nodes[i].z...)
    params.nodes[i].dip = intp_dip(params.nodes[i].z...)
    params.nodes[i].holeradius = intp_holeradius(params.nodes[i].z...)
  end
end

function surveys2df(data::Dict{Symbol, Any})
  syms = [:survey_type, :depth, :inclination, :heading, :diameter]
  df = DataFrame(Any[DataType[], Float64[], Float64[], Float64[], Float64[]], syms)
  for surv in data[:surveys]
    push!(df, [typeof(surv)  surv.depth surv.inclination surv.heading surv.diameter])
  end
  df
end


text:
"""
    QuarticNonconvex{T}

The scalar conservation law
``
  \\partial_t u + \\partial_x ( u(t,x)^4 - 10 u(t,x)^2 + 3 u(t,x) ) = 0
``
in one space dimensions using `T` as scalar type.
"""
struct QuarticNonconvex{T} <: ScalarBalanceLaw{T,1} end

function QuarticNonconvex(T=Float64)
  QuarticNonconvex{T}()
end

function Base.show(io::IO, model::QuarticNonconvex{T}) where {T}
  print(io, "Scalar conservation law {T=", T, "} with flux f(u) = u^4 - 10 u^2 + 3 u")
end

"""
    flux(u, model::QuarticNonconvex)

Compute the flux of `u` for `model`.
"""
@inline flux(u, model::QuarticNonconvex) = u^2 * (u^2 - 10) + 3u

"""
    speed(u::Real, model::QuarticNonconvex)

Compute the speed f'(`u`) for `model`.
"""
@inline speed(u, model::QuarticNonconvex) = 4*u^3 - 20*u + 3

function min_max_speed(ul::T, ur::T, model::QuarticNonconvex{T}) where T<:Real
  min_speed = zero(T)
  max_speed = zero(T)

  if ur > 1.2909944487358056284 && -2.5819888974716112568 < ul < 1.2909944487358056284
    min_speed = T(-14.213259316477408379)
  elseif (ur > 1.2909944487358056284 && (ul >= 1.2909944487358056284 || ul <= -2.5819888974716112568)) || (ur <= 1.2909944487358056284 && 2 * ul + ur + sqrt(20 - 3 * max(ur, 1.2909944487358056284)^2) <= 0) || 1.2909944487358056284 + ur <= 0
    min_speed = 3 + 4 * ul * (-5 + ul^2)
  else
    min_speed = 3 + 4 * ur * (-5 + ur^2)
  end

  if (1.2909944487358056284 + ul <= 0 && 1.2909944487358056284 + ur > 0 && ur < 1.2909944487358056284) || (1.2909944487358056284 + ul < 0 && 1.2909944487358056284 <= ur < 2.5819888974716112568)
    max_speed = T(20.213259316477408379)
  elseif (1.2909944487358056284 + ul == 0 && ur == 2.5819888974716112568) || (1.2909944487358056284 + ul >= 0 && (ur ≈ 1.2909944487358056284 || (1.2909944487358056284 < ur < 2.5819888974716112568 && sqrt(20 - 3 * min(ur, 2.5819888974716112568)^2) >= 2 * ul + ur))) || (ur < 1.2909944487358056284 && 1.2909944487358056284 + ul > 0)
    max_speed = 3 + 4 * ul * (-5 + ul^2)
  else
    max_speed = 3 + 4 * ur * (-5 + ur^2)
  end

  min_speed, max_speed
end

function max_abs_speed(ul::T, ur::T, model::QuarticNonconvex{T}) where T<:Real
  min_speed, max_speed = min_max_speed(ul, ur, model)
  max(abs(min_speed), abs(max_speed))
end

################################################################################

"""
    (::EnergyConservativeFlux)(uₗ::T, uᵣ::T, model::QuarticNonconvex{T}) where {T}

Compute the energy (L₂ entropy) conservative flux between `uₗ` and `uᵣ` for `model`.
"""
function (flux::EnergyConservativeFlux)(ul::T, ur::T, model::QuarticNonconvex{T}) where T<:Real
  (ul^4 + ur*ul^3 + ur^2*ul^2 + ur^3*ul + ur^4) / 5 - 10 * (ul^2 + ul*ur + ur^2) / 3 + 3 * (ul + ur) / 2
end


text:
function solve()
    n = parse(Int, readline())
    s = [parse(Int, x) for x in readline()]
    x = readline()

    dp = Set{Int}(0)

    for i in n:-1:1
        if x[i] == 'T'
            dp = Set{Int}(r for r in 0:6 if 10r % 7 ∈ dp || (10r+s[i]) % 7 ∈ dp)
        else
            dp = Set{Int}(r for r in 0:6 if 10r % 7 ∈ dp && (10r+s[i]) % 7 ∈ dp)
        end
    end
    0 ∈ dp ? "Takahashi" : "Aoki"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Documenter, ColorSchemes, Luxor

# before generating the documents, build the list of colorschemes
# and save it in a .md file

function drawscheme(schemename;
        swatchwidth=500,
        swatchheight=20)

    # create swatch
    cols = colorschemes[schemename].colors
    l = length(cols)
    Drawing(swatchwidth, swatchheight, :svg)
    setline(0.5)
    origin()
    t = Tiler(swatchwidth, swatchheight, 1, l, margin=0)
    for (i, c) in enumerate(cols)
        sethue(c)
        box(t, i, :fillstroke)
    end
    finish()

    swatch = svgstring()

    # fill the template
    schemetemplate = """

:$(schemename)

$(swatch)

    """
    return schemetemplate
end

function write_all_schemes_in_category(f, category)
    # start a table for this category
    write(f,
        """
        ## Category: $category
        ```@raw html
        <div id="grid-col">
        """)

    # add row for each scheme
    schemes = filter(s -> occursin(category, colorschemes[s].category), collect(keys(colorschemes))) # hide
    for scheme in sort(schemes)
        write(f, drawscheme(scheme))
    end

    # end table for this category
    write(f, """
        </div>
        ```
        """)
end

open("docs/src/catalogue.md", "w") do f
    write(f, """
    # Catalogue of all colorschemes

    This section lists all the colorschemes provided. Refer
    to the [References](@ref) section for links to the original
    sources for these.
    """)
    # a table for each category
    write_all_schemes_in_category(f, "colorbrewer")
    write_all_schemes_in_category(f, "cmocean")
    write_all_schemes_in_category(f, "scientific")
    write_all_schemes_in_category(f, "matplotlib")
    write_all_schemes_in_category(f, "gnuplot")
    write_all_schemes_in_category(f, "colorcet")
    write_all_schemes_in_category(f, "seaborn")
    write_all_schemes_in_category(f, "tableau")
    write_all_schemes_in_category(f, "cvd")
    write_all_schemes_in_category(f, "flags")
    write_all_schemes_in_category(f, "general")
end

makedocs(
    modules = [ColorSchemes],
    sitename = "ColorSchemes",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/colorschemes-docs.css"]),
    pages = Any[
        "Introduction"      => "index.md",
        # "Gallery"           => "gallery.md",
        "Catalogue"         => "catalogue.md",
        "Basic usage"       => "basics.md",
        "Finding colors"    => "finding.md",
        "Plotting"          => "plotting.md",
        "Images"            => "images.md",
        "References"        => "references.md",
        "Index"             => "functionindex.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaGraphics/ColorSchemes.jl.git",
    target = "build"
)


text:
"""
Run existing CellML model and save output membrane potential Vm
"""

using CellMLToolkit, DifferentialEquations, OrdinaryDiffEq, Sundials
using PyPlot
PyPlot.pygui(true)

function LoadAndRun(model_name::String; show_params=false, show_states=false)
    """
    Load and run model with model name `model_name`
    """
    # path to cellml models
    cellml_path = ".//data//CellML_models//$model_name"

    # CellModel factory function creates an ODESystem
    ml = CellMLToolkit.CellModel(cellml_path; dependency=false)

    # model parameters
    if show_params
        params = list_params(ml)
        @info params
    end

    # model states
    states = list_states(ml)
    if show_states
        @info states
    end

    # find index of membrane voltage in list of state variables
    idx = findfirst([occursin("V", string(s)) for s in states])
    @info "Index of voltage variable:", idx

    # generate and solve ODEProblem
    tspan = (0., 0.4)
    ts = collect(0:1e-5:tspan[2])
    prob = CellMLToolkit.ODEProblem(ml, tspan);
    sol = solve(prob, TRBDF2(), saveat=ts)
    # extract membrane voltage
    Vm = sol[idx, :]

    # plot
    f, ax = plt.subplots(figsize=(10, 5))
    ax.set_ylabel("Voltage (mV)")
    ax.set_xlabel("Time (ms)")
    ax.plot(sol.x, Vm, lw=2)
    f.tight_layout()
    PyPlot.show()
end

model_name = "severi_fantini_charawi_difrancesco_2012.cellml"
LoadAndRun(model_name; show_params=true, show_states=true)

max_stars_repo_name: haganenoneko/EphysAnalysisTools

max_issues_repo_name: haganenoneko/EphysAnalysisTools

max_forks_repo_name: haganenoneko/EphysAnalysisTools


text:
using Flux3D
using Test

include("utils.jl")

@testset "Flux3D" begin

    @testset "Representation" begin
        include("rep.jl")
    end

    @testset "Model" begin
        include("models.jl")
    end

    @testset "Dataset" begin
        include("dataset.jl")
    end

end # testset Flux3D


text:
using PipeLayout.JuncLoc

using JuMP
using SCS

_scs = JuMP.optimizer_with_attributes(
    SCS.Optimizer,
    "eps" => 1e-6,
    "verbose" => 0
)

@testset "solve junction location for three terminals (SOC)" begin
    # equilateral triangle
    nodes = [Node(0,0), Node(40,0), Node(20, sqrt(3)/2*40)]
    demand = [10, 10, -20]
    bounds = fill(Bounds(40, 80), 3)
    diams = [Diameter(t...) for t in [(0.6, 0.5), (0.8, 0.8),(1.0, 1.2)]]

    # single Steiner node with dummy location
    topo = Topology(vcat(nodes, [Node(20, 20)]),
                    [Arc(3,4), Arc(4,1), Arc(4,2)])

    @testset "little flow, smallest diameter, Steiner node in center" begin
        inst = Instance(nodes, demand, bounds, diams, ploss_coeff_nice)
        model, x, y, t, π = JuncLoc.make_soc(
            inst, topo, JuncLoc.SOC(_scs))
        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)

        @test status in (MOI.OPTIMAL, MOI.ITERATION_LIMIT)

        xsol, ysol = JuMP.value.(x), JuMP.value.(y)
        for i=1:3 # fixed terminals
            @test xsol[i] ≈ nodes[i].x atol=0.005
            @test ysol[i] ≈ nodes[i].y atol=0.005
        end
        @test xsol[4] ≈ 20 atol=0.01
        @test ysol[4] ≈ sqrt(3)/6*40 atol=0.01

        tsol = JuMP.value.(t)
        @test sum(tsol[1]) ≈ sum(tsol[2]) atol=0.01
        @test sum(tsol[1]) ≈ sum(tsol[3]) atol=0.01
    end

    @testset "more flow, mixed diameter, Steiner node towards source" begin
        inst = Instance(nodes, 20*demand, bounds, diams, ploss_coeff_nice)
        model, x, y, t, π = JuncLoc.make_soc(
            inst, topo, JuncLoc.SOC(_scs))
        JuMP.optimize!(model)
	    status = JuMP.termination_status(model)

        @test status == MOI.OPTIMAL

        xsol, ysol = JuMP.value.(x), JuMP.value.(y)
        for i=1:3 # fixed terminals
            @test xsol[i] ≈ nodes[i].x atol=0.005
            @test ysol[i] ≈ nodes[i].y atol=0.005
        end
        @test xsol[4] ≈ 20 atol=0.1
        @test ysol[4] >= sqrt(3)/6*40 # move near source

        tsol = JuMP.value.(t)
        @test tsol[2] ≈ tsol[3] atol=0.01
    end

    @testset "using the optimize function to solve" begin
        inst = Instance(nodes, 20*demand, bounds, diams, ploss_coeff_nice)
        result = optimize(inst, topo, JuncLoc.SOC(_scs))
        @test result.status == MOI.OPTIMAL

        sol = result.sol
        toposol = Topology(sol.nodes, topo.arcs)
        L = pipelengths(toposol)
        c = [d.cost for d in diams]
        obj = L' * sol.lsol * c
        @test result.value ≈ obj[1] atol=0.005
    end
end


text:
function get_current_state(obj::ODEStateInterpolator)
    return jcall(obj, "getCurrentState", ODEStateAndDerivative, ())
end

function get_interpolated_state(obj::ODEStateInterpolator, arg0::jdouble)
    return jcall(obj, "getInterpolatedState", ODEStateAndDerivative, (jdouble,), arg0)
end

function get_previous_state(obj::ODEStateInterpolator)
    return jcall(obj, "getPreviousState", ODEStateAndDerivative, ())
end

function is_current_state_interpolated(obj::ODEStateInterpolator)
    return jcall(obj, "isCurrentStateInterpolated", jboolean, ())
end

function is_forward(obj::ODEStateInterpolator)
    return jcall(obj, "isForward", jboolean, ())
end

function is_previous_state_interpolated(obj::ODEStateInterpolator)
    return jcall(obj, "isPreviousStateInterpolated", jboolean, ())
end