text:
using Clapeyron, Test

@testset "All tests" begin
    include("test_database.jl")
    include("test_models.jl")
    include("test_methods.jl")
end


text:
# Laguerre and Hermite fail due to linear algebra problems in BigFloat
supports_approximation(::Laguerre{BigFloat}) = false
supports_approximation(::Laguerre{Double64}) = false
supports_approximation(::Hermite{BigFloat}) = false
supports_approximation(::Hermite{Double64}) = false
# It is difficult to do approximation in subsets and operated sets generically
supports_approximation(::Subdictionary) = false
supports_approximation(::OperatedDict) = false
supports_approximation(dict::TensorProductDict) =
    mapreduce(supports_approximation, &, components(dict))
# Monomials and rationals have no associated domain
supports_approximation(::Monomials) = false
supports_approximation(::RationalFunctions) = false

# Make a simple periodic function for Fourier and other periodic sets
suitable_function(::Fourier) =  x -> 1/(10+cos(2*pi*x))
# The function has to be periodic and even symmetric
suitable_function(::CosineSeries) =  x -> 1/(10+cos(2*pi*x))
# The function has to be periodic and odd symmetric
suitable_function(::SineSeries) =  x -> x^3*(1-x)^3
# We use a function that is smooth and decays towards infinity
suitable_function(::Laguerre) = x -> 1/(1000+(2x)^2)
suitable_function(::Hermite) = x -> 1/(1000+(2x)^2)

suitable_interpolation_grid(basis::TensorProductDict) =
    ProductGrid(map(suitable_interpolation_grid, components(basis))...)
suitable_interpolation_grid(basis::SineSeries) = MidpointEquispacedGrid{domaintype(basis)}(length(basis), 0, 1)
suitable_interpolation_grid(basis::WeightedDict) = suitable_interpolation_grid(superdict(basis))
suitable_interpolation_grid(basis::OperatedDict) = suitable_interpolation_grid(src(basis))

ext: jl


text:
"""
Node Genes
"""

@enum NodeType INPUT=-1 HIDDEN=0 OUTPUT=1

mutable struct Node
    id::Int
    dtype::NodeType
    activation::ActivationType
    bias::Float64
    connections::Dict{String, Vector{Connection}}

    function Node(id::Int, dtype::NodeType, activation::ActivationType=sigmoid, bias::Float64=rand())
        connections = Dict{String, Vector{Connection}}("in" => [], "out" => [])
        new(id, dtype,activation, bias, connections)
    end

    function Node(id::Int, dtype::NodeType, activation, bias::Float64, connections)
        # @assert typeof(activation) == ActivationType
        new(id, dtype, activation, bias, connections)
    end
end

Base.copy(n::Node) = Node(n.id, n.dtype, n.activation, n.bias, deepcopy(n.connections))

function get_connection(node::Node, from::Int, to::Int, mode::String)
    for conn in node.connections[mode]
        if conn.input == from && conn.output == to
            return conn
        end
    end
end


text:
using POMDPTools
using Documenter

DocMeta.setdocmeta!(POMDPTools, :DocTestSetup, :(using POMDPTools); recursive=true)

makedocs(;
    modules=[POMDPTools],
    authors="Zachary Sunberg <zachary.sunberg@colorado.edu> and contributors",
    repo="https://github.com/JuliaPOMDP/POMDPTools.jl/blob/{commit}{path}#{line}",
    sitename="POMDPTools.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaPOMDP.github.io/POMDPTools.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaPOMDP/POMDPTools.jl",
    devbranch="main",
)


text:
module NumericExtensions

    warn("""
    A majority of the functionalities implemented in this package has been available & optimized in Julia Base.
    This package is no longer actively developed (except for occasional bugfixes).
    Please consider removing the dependency on NumericExtensions and using the functions in JuliaBase instead.""")

    using ArrayViews
    using NumericFuns
    using Compat

    # size information
    import Base: size, length, ndims, isempty, stride

    # iteration and indexing
    import Base: start, next, done, getindex, setindex!

    # common operations
    import Base: show, copy, similar, pointer, parent, diagm

    # higher-level map & reduction functions to be extended
    import Base: map, map!, reduce, reducedim, mapreduce, foldl, foldr

    # arithmetic functions
    import Base: +, *, \, /, ==

    # reduction functions to be extended
    import Base: sum, sum!, prod, prod!, maximum, maximum!, minimum, minimum!, dot, norm

    # scan functions to be extended
    import Base: cumsum, cumsum!, cummax, cummin, cumprod

    # statistics
    import Base: mean, mean!, var, varm, std, stdm

    # from ArrayViews
    import ArrayViews: view, ellipview

    # from NumericFuns
    import NumericFuns: evaluate, result_type, logsumexp

    export
        # reexport views
        view, ellipview,

        # re-export functor-related methods
        fptype, arithtype, evaluate, result_type,

        # re-export mathfuns
        sqr, rcp, rsqrt, rcbrt, xlogx, xlogy,
        sigmoid, logit, softplus, invsoftplus, logsumexp,

        # shapes
        mapshape, maplength,

        # views
        view, unsafe_view,

        # mathfuns

        logit, xlogx, xlogy, logistic, invlogistic, softplus, invsoftplus,

        # functors
        Functor, @functor1, @functor2, evaluate, fptype,

        # map
        map, map!, map1!, mapdiff, mapdiff!,
        add!, subtract!, multiply!, divide!, negate!, rcp!,
        sqrt!, abs!, sqr!, abs2!, pow!, exp!, log!,
        floor!, ceil!, round!, trunc!,
        absdiff, absdiff!, sqrdiff, sqrdiff!, fma, fma!,

        # vbroadcast
        vbroadcast, vbroadcast!, vbroadcast1!,
        badd, badd!, bsubtract, bsubtract!, bmultiply, bmultiply!, bdivide, bdivide!,

        # diagop
        add_diag!, add_diag, set_diag!, set_diag,

        # reduce
        mapreduce_fdiff, sumfdiff, maxfdiff, minfdiff, meanfdiff,
        sumabs, maxabs, minabs, meanabs, sumsq, meansq,
        sumabsdiff, maxabsdiff, minabsdiff, meanabsdiff,
        sumsqdiff, meansqdiff, sumxlogx, sumxlogy, entropy,

        foldl_fdiff, foldr_fdiff,

        # reducedim
        reducedim!, mapreducedim, mapreducedim!, mapreducedim_fdiff, mapreducedim_fdiff!,
        foldl!, foldr!, sum!, maximum!, minimum!, mean!, dot!,

        sumfdiff!, maxfdiff!, minfdiff!, meanfdiff!,
        sumabs!, maxabs!, minabs!, meanabs!, sumsq!, meansq!,
        sumabsdiff!, maxabsdiff!, minabsdiff!, meanabsdiff!,
        sumsqdiff!, meansqdiff!, sumxlogx!, sumxlogy!, entropy!,

        # norms
        vnorm, vnorm!, vnormdiff, vnormdiff!, normalize, normalize!,

        # scan
        scan, scan!, mapscan, mapscan!,
        cumsum!, cummax!, cummin!,

        # statistics
        var!, varm!, std!, stdm!,
        logsumexp, logsumexp!, softmax, softmax!,

        # weightsum
        wsum, wsum!, wsumfdiff, wsumfdiff!,
        wsumabs, wsumabs!, wsumabsdiff, wsumabsdiff!,
        wsumsq, wsumsq!, wsumsqdiff, wsumsqdiff!,

        # utils
        eachrepeat, sortindexes, groupindexes,

        # benchmark
        BenchmarkTable, nrows, ncolumns, add_row!

    # codes

    include("common.jl")
    include("shapes.jl")
    include("codegen.jl")

    include("mapkernels.jl")
    include("map.jl")

    include("vbroadcast.jl")
    include("diagop.jl")
    include("utils.jl")

    include("rkernels.jl")
    include("reduce.jl")
    include("folddim.jl")
    include("reducedim.jl")

    include("norms.jl")
    include("scan.jl")
    include("statistics.jl")
    include("wsum.jl")

    include("benchmark.jl")

    include("deprecates.jl")
end


text:
module WeightTest

using TestSetup, OnlineStats, Base.Test
O = OnlineStats

@testset "Weighting" begin
    @testset "EqualWeight" begin
        w = EqualWeight()
        O.updatecounter!(w)
        @test O.weight(w) == 1.0
        O.updatecounter!(w, 5)
        @test O.weight(w, 5) == 5 / 6
    end
    @testset "ExponentialWeight" begin
        w = ExponentialWeight(.5)
        O.updatecounter!(w)
        @test O.weight(w) == 0.5
        O.updatecounter!(w, 5)
        @test O.weight(w) == 0.5
        @test ExponentialWeight(5).λ == ExponentialWeight(1/3).λ
    end
    @testset "BoundedEqualWeight" begin
        w = BoundedEqualWeight(.1)
        O.updatecounter!(w)
        @test O.weight(w, 1) == 1.0
        O.updatecounter!(w, 100)
        @test O.weight(w) == 0.1
        @test BoundedEqualWeight(5).λ == BoundedEqualWeight(1/3).λ
    end
    @testset "LearningRate" begin
        w = LearningRate(.6)
        O.updatecounter!(w)
        @test O.weight(w) == 1.0
        O.updatecounter!(w, 5)
        @test O.weight(w, 5) == 2 ^ -.6
        @test O.nups(w) == 2
    end
    @testset "LearningRate2" begin
        w = LearningRate2(0.5)
        O.updatecounter!(w)
        @test O.weight(w) == 1.0
        O.updatecounter!(w, 5)
        @test O.weight(w, 5) == 1 / (1 + .5)
    end
end  # facts
end  # module


text:
using Benchmark
using RingArrays

function print_seperator(file::IOStream, name::AbstractString)
    seperator = "################################################################################"
    println(file)
    println(file, seperator)
    println(file, "# ", name)
    println(file, seperator)
    println(file)
end

open(joinpath(dirname(@__FILE__), "benchmark_result"), "w") do benchmark_file

    num_test = 100
    m_b = 10
    b_s = (10,)
    d_l = 100000
    big_data = rand(Int, d_l)

    print_seperator(benchmark_file, "constructor")

    constructor_ring_bench() = RingArray{Int, 1}(max_blocks=m_b, block_size=b_s, data_length=d_l)
    constructor_normal_array_bench() = Array{Int, 1}()

    println(benchmark_file, benchmark(constructor_ring_bench, "constructor ring", "RingArray{Int, 1}(max_blocks=m_b, block_size=b_s, data_length=d_l)", num_test), "\n")
    println(benchmark_file, benchmark(constructor_normal_array_bench, "constructor normal array", "Array{Int, 1}()", num_test), "\n")
    println(benchmark_file, compare([constructor_ring_bench, constructor_normal_array_bench], num_test), "\n")

    print_seperator(benchmark_file, "load_block")

    function loading_no_gc()
        ring = RingArray{Int, 1}(max_blocks=m_b, block_size=b_s, data_length=d_l)
        for i in 1:b_s[1]:d_l
            prev_write = ring.next_write
            load_block(ring, big_data[i:i+b_s[1]-1])
            let
                view = ring[i:b_s[1]-1]
            end
            ring.num_users[prev_write] = 0
        end
    end

    function loading_with_gc()
        ring = RingArray{Int, 1}(max_blocks=m_b, block_size=b_s, data_length=d_l)
        for i in 1:b_s[1]:d_l
            load_block(ring, big_data[i:i+b_s[1]-1])
            let
                view = ring[i:b_s[1]-1]
            end
        end
    end

    println(benchmark_file, benchmark(loading_no_gc, "loading with no gc", "load_block", num_test), "\n")
    println(benchmark_file, benchmark(loading_with_gc, "loading with gc", "load_block", num_test), "\n")
    println(benchmark_file, compare([loading_no_gc, loading_with_gc], num_test), "\n")

    empty_ring = RingArray{Int, 1}(max_blocks=m_b, block_size=b_s, data_length=d_l)
    full_ring = RingArray{Int, 1}(max_blocks=m_b, block_size=b_s, data_length=d_l)
    end_ring = RingArray{Int, 1}(max_blocks=m_b, block_size=b_s, data_length=d_l)

    for i in 1:b_s[1]:b_s[1]*m_b
        load_block(full_ring, big_data[i:i+b_s[1]-1])
    end

    for i in 1:b_s[1]:d_l
        load_block(end_ring, big_data[i:i+b_s[1]-1])
    end

    print_seperator(benchmark_file, "size")

    size_empty_bench() = size(empty_ring)
    size_full_bench() = size(full_ring)
    size_end_bench() = size(end_ring)
    size_normal_array_bench() = size(big_data)

    println(benchmark_file, benchmark(size_empty_bench, "size empty ring", "size(empty_ring)", num_test), "\n")
    println(benchmark_file, benchmark(size_full_bench, "size full ring", "size(full_ring)", num_test), "\n")
    println(benchmark_file, benchmark(size_end_bench, "size end ring", "size(end_ring)", num_test), "\n")
    println(benchmark_file, benchmark(size_normal_array_bench, "size normal array", "size(big_data)", num_test), "\n")
    println(benchmark_file, compare([size_empty_bench, size_full_bench, size_end_bench, size_normal_array_bench], num_test), "\n")

    print_seperator(benchmark_file, "checkbounds")

    #checkbounds_empty_bench() = checkbounds(empty_ring, empty_ring.range.stop) throws an error
    checkbounds_full_bench() = checkbounds(full_ring, full_ring.range.stop)
    checkbounds_end_bench() = checkbounds(end_ring, end_ring.range.stop)
    checkbounds_normal_array_bench() = checkbounds(big_data, end_ring.range.stop)

    #println(benchmark_file, benchmark(checkbounds_empty_bench, "checkbounds empty ring", "checkbounds(empty_ring, empty_ring.range.stop)", num_test), "\n")
    println(benchmark_file, benchmark(checkbounds_full_bench, "checkbounds full ring", "checkbounds(full_ring, full_ring.range.stop)", num_test), "\n")
    println(benchmark_file, benchmark(checkbounds_end_bench, "checkbounds end ring", "checkbounds(end_ring, end_ring.range.stop)", num_test), "\n")
    println(benchmark_file, benchmark(checkbounds_normal_array_bench, "checkbounds normal array", "checkbounds(big_data, end_ring.range.stop)", num_test), "\n")
    println(benchmark_file, compare([checkbounds_full_bench, checkbounds_end_bench, checkbounds_normal_array_bench], num_test), "\n")

    print_seperator(benchmark_file, "getindex")

    getindex_full_start_bench() = full_ring[full_ring.range.start]
    getindex_full_stop_bench() = full_ring[full_ring.range.stop]
    getindex_end_start_bench() = end_ring[end_ring.range.start]
    getindex_end_stop_bench() = end_ring[end_ring.range.stop]
    getindex_normal_array_bench() = big_data[end_ring.range.stop]

    println(benchmark_file, benchmark(getindex_full_start_bench, "getindex full start ring", "full_ring[full_ring.range.start]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_full_stop_bench, "getindex full stop ring", "full_ring[full_ring.range.stop]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_end_start_bench, "getindex end start ring", "end_ring[end_ring.range.start]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_end_stop_bench, "getindex end stop ring", "end_ring[end_ring.range.stop]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_normal_array_bench, "getindex normal array", "big_data[end_ring.range.stop]", num_test), "\n")
    println(benchmark_file, compare([getindex_full_start_bench, getindex_full_stop_bench, getindex_end_start_bench, getindex_end_stop_bench, getindex_normal_array_bench], num_test), "\n")

    print_seperator(benchmark_file, "getindex range")

    getindex_range_full_all_bench() = full_ring[full_ring.range]
    getindex_range_full_small_bench() = full_ring[full_ring.range.stop:full_ring.range.stop]
    getindex_range_end_all_bench() = end_ring[end_ring.range]
    getindex_range_end_small_bench() = end_ring[end_ring.range.stop:end_ring.range.stop]
    getindex_range_normal_array_all_bench() = big_data[end_ring.range]
    getindex_range_normal_array_small_bench() = big_data[end_ring.range.stop:end_ring.range.stop]

    println(benchmark_file, benchmark(getindex_range_full_all_bench, "getindex range full all ring", "full_ring[full_ring.range]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_range_full_small_bench, "getindex range full small ring", "full_ring[full_ring.range.stop:full_ring.range.stop]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_range_end_all_bench, "getindex range end all ring", "end_ring[end_ring.range]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_range_end_small_bench, "getindex range end small ring", "end_ring[end_ring.range:end_ring.range.stop]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_range_normal_array_all_bench, "getindex range normal all array", "big_data[end_ring.range]", num_test), "\n")
    println(benchmark_file, benchmark(getindex_range_normal_array_small_bench, "getindex range normal small array", "big_data[end_ring.range.stop:end_ring.range.stop]", num_test), "\n")
    println(benchmark_file, compare([getindex_range_full_all_bench, getindex_range_full_small_bench, getindex_range_end_all_bench, getindex_range_end_small_bench, getindex_range_normal_array_all_bench, getindex_range_normal_array_small_bench], num_test), "\n")

    print_seperator(benchmark_file, "overall")

    standard() = 0
    println(benchmark_file, compare(
        [
            standard,
            constructor_ring_bench, constructor_normal_array_bench,
            loading_no_gc, loading_with_gc,
            size_empty_bench, size_full_bench, size_end_bench, size_normal_array_bench,
            checkbounds_full_bench, checkbounds_end_bench, checkbounds_normal_array_bench,
            getindex_full_start_bench, getindex_full_stop_bench, getindex_end_start_bench, getindex_end_stop_bench, getindex_normal_array_bench,
            getindex_range_full_all_bench, getindex_range_full_small_bench, getindex_range_end_all_bench, getindex_range_end_small_bench, getindex_range_normal_array_all_bench, getindex_range_normal_array_small_bench,
        ],
        num_test), "\n")

end


text:
using DataFrames, Dash, DashHtmlComponents, DashCoreComponents, PlotlyJS

df4 = DataFrame(Dict(("Col $(i)" => rand(30)) for i = 1:6))

function create_figure(df4, x_col, y_col, selectedpoints, selectedpoints_local)
    if selectedpoints_local != nothing
        ranges = selectedpoints_local[:range]
        selection_bounds = Dict(
            "x0" => ranges[:x][1],
            "x1" => ranges[:x][2],
            "y0" => ranges[:y][1],
            "x0" => ranges[:y][2],
        )
    else
        selection_bounds = Dict(
            "x0" => minimum(df4[:, x_col][1]),
            "x1" => minimum(df4[:, x_col][2]),
            "y0" => minimum(df4[:, y_col][1]),
            "x0" => minimum(df4[:, x_col][2]),
        )
    end

    return Plot(
        df4,
        x = df4[:, x_col],
        y = y_col,
        mode = "markers+text",
        marker_size = 20,
        text = 1:size(df4)[1],
        customdata = 1:size(df4)[1],
        selectedpoints = selectedpoints,
        unselected = (
            marker = (opacity = 0.3, textfont = (color = "rgba(0,0,0,0)"))
        )
    )
end

app = dash()

app.layout = html_div() do
    html_div(
        dcc_graph(id = "graph1"),
        style = (width = "32%", display = "inline-block"),
    ),
    html_div(
        dcc_graph(id = "graph2"),
        style = (width = "32%", display = "inline-block"),
    ),
    html_div(
        dcc_graph(id = "graph3"),
        style = (width = "32%", display = "inline-block"),
    )
end

callback!(
    app,
    Output("graph1", "figure"),
    Output("graph2", "figure"),
    Output("graph3", "figure"),
    Input("graph1", "selectedData"),
    Input("graph2", "selectedData"),
    Input("graph3", "selectedData"),
) do selection1, selection2, selection3
    selectedpoints = 1:size(df4)[1]

    for selected_data in [selection1, selection2, selection3]
        if selected_data != nothing
            selectedpoints = [p[:customdata] - 1 for p in selected_data.points]
        end
    end

    return create_figure(df4, "Col 1", "Col 2", selectedpoints, selection1),
           create_figure(df4, "Col 3", "Col 4", selectedpoints, selection2),
           create_figure(df4, "Col 5", "Col 6", selectedpoints, selection3)
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
struct FlowSample{T}
    data::Dict{String, Vector{T}}
    params::Dict{String, String}
end

const opt_params = [
("Machine", "\$CYT"),
("Begin Time", "\$BTIM"),
("End Time", "\$ETIM"),
("Date", "\$DATE"),
("File", "\$FIL"),
("Volume run", "\$VOL")
]

function Base.show(io::IO, f::FlowSample)
    spacing = " "^4
    print(io, typeof(f))

    for pair in opt_params
        if haskey(f.params, pair[2])
            print(io, "\n", spacing, "$(pair[1]): $(f.params[pair[2]])")
        end
    end
    print(io, "\n", spacing, "Axes:")
    n_params = parse(Int, f.params["\$PAR"])
    for i in 1:n_params
        print(io, "\n", spacing, spacing, "$(f.params["\$P$(i)N"])")
        if haskey(f.params, "\$P$(i)S")
            print(io, " ($(f.params["\$P$(i)N"]))")
        end
    end
end

# Implement most important parts of Dict interface
Base.length(f::FlowSample)  = length(f.data)
Base.haskey(f::FlowSample, x) = haskey(f.data, x)
Base.getindex(f::FlowSample, key) = f.data[key]
Base.keys(f::FlowSample) = keys(f.data)
Base.values(f::FlowSample) = values(f.data)
Base.iterate(iter::FlowSample) = Base.iterate(iter.data)
Base.iterate(iter::FlowSample, state) = Base.iterate(iter.data, state)


text:
# Julia wrapper for header: /usr/include/scip/cons_orbisack.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrOrbisack(scip)
    ccall((:SCIPincludeConshdlrOrbisack, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPseparateCoversOrbisack(scip, cons, sol, vars1, vars2, nrows, infeasible, ngen)
    ccall((:SCIPseparateCoversOrbisack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{UInt32}, Ptr{Cint}), scip, cons, sol, vars1, vars2, nrows, infeasible, ngen)
end

function SCIPcheckSolutionOrbisack(scip, sol, vars1, vars2, nrows, printreason, feasible)
    ccall((:SCIPcheckSolutionOrbisack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SOL}, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}, Cint, UInt32, Ptr{UInt32}), scip, sol, vars1, vars2, nrows, printreason, feasible)
end

function SCIPcreateConsOrbisack(scip, cons, name, vars1, vars2, nrows, ispporbisack, isparttype, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsOrbisack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}, Cint, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, vars1, vars2, nrows, ispporbisack, isparttype, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicOrbisack(scip, cons, name, vars1, vars2, nrows, ispporbisack, isparttype)
    ccall((:SCIPcreateConsBasicOrbisack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}, Cint, UInt32, UInt32), scip, cons, name, vars1, vars2, nrows, ispporbisack, isparttype)
end


text:
using Flux

export forward, back_propagation, loss, get_gradient, para_update!
export QDNNModel

struct QDNNModel{T}
    layers::Array{L, 1} where {L<:Layer{T}}
    # loss::Function
end

function forward(model::QDNNModel{T}, x::Array{T}) where {T}
    y = x
    Z = [x]
    for layer in model.layers
        y = QDNN.forward(layer, y)
        push!(Z, y)
    end
    return Z
end

function back_propagation(model::QDNNModel{T}, x::Array{T}, L_y::Array{T}) where {T}
    Z = forward(model, x)
    L_x = L_y
    nl = size(model.layers, 1)
    grad = Vector{Vector{T}}[]

    for l = nl:-1:1
        layer = model.layers[l]
        L_xx, L_ww, L_bb = QDNN.back_propagation(layer, Z[l])
        L_w = L_ww * L_x
        if L_bb != 0
            L_b = L_bb * L_x
        else
            L_b = T[]
        end
        pushfirst!(grad, [L_w, L_b])

        L_x = L_xx * L_x
    end

    return grad
end

# function MSE(y1::Vector{T}, y2::Vector{T}) where {T}
#     n = size(y1, 1)
#     return sum((y1-y2).^2)/n, 2/n .* (y1-y2)
# end

function SSE(y1::Vector{T}, y2::Vector{T}) where {T}
    return sum((y1-y2).^2), 2 .* (y1-y2)
end

function loss(qm::QDNNModel{T}, x::Vector{T}, y::Vector{T}) where {T}
    fwd = forward(qm, x)
    y_p = fwd[end]
    l, L_y = SSE(y_p, y)
    return l
end

function loss(qm::QDNNModel{T}, data_x::Vector{Vector{T}}, data_y::Vector{Vector{T}}) where {T}
    nData = size(data_x, 1)
    l = zeros(nData)
    for i = 1:nData
        x = data_x[i]
        y = data_y[i]
        l_x = loss(qm, x, y)
        l[i] = l_x
    end
    return l
end

function get_gradient(qm::QDNNModel{T}, x::Vector{T}, y::Vector{T}) where {T}
    fwd = forward(qm, x)
    y_p = fwd[end]
    l, L_y = SSE(y_p, y)

    grad = back_propagation(qm, x, L_y)
    return grad
end

function get_gradient(qm::QDNNModel{T}, data_x::Vector{Vector{T}}, data_y::Vector{Vector{T}}, nbatch::Integer) where T
    grad = nothing
    nData = size(data_x, 1)
    for j = 1:nbatch
        i = rand(1:nData)
        x = data_x[i]
        y = data_y[i]
        gradi = get_gradient(qm, x, y)
        if grad == nothing
            grad = gradi
        else
            grad += gradi
        end
    end
    grad = grad ./ nbatch
    return grad
end

function para_update!(qm::QDNNModel{T}, grad::Vector{Vector{Vector{T}}}, opt) where {T}
    lys = qm.layers
    nl = size(lys, 1)
    for i = 1:nl
        grad_w = grad[i][1]
        grad_b = grad[i][2]
        Flux.Optimise.update!(opt, lys[i].params, grad_w)
        Flux.Optimise.update!(opt, lys[i].bias, grad_b)
    end
end


text:
####
#### Graph library
####

## This is a collection of some classes of graphs.
## Note that the empty graph is given by Graph(n), so no function
## for the empty graph is provided here.

"""
    cycle_graph(n::Integer)

Return the cycle graph of order `n`.
"""
function cycle_graph(n::Integer)
    g = Graph(n)
    for i in 1:n-1
        add_edge!(g, i, i + 1)
    end
    add_edge!(g, 1, n)
    return g
end

"""
    complete_graph(n::Integer)

Return the complete graph of order `n`.
"""
function complete_graph(n::Integer)
    g = Graph(n)
    for i in 1:n, j in 1:i-1
        add_edge!(g, i, j)
    end
    return g
end

"""
    star_graph(n::Integer)

The star graph of order `n`. This is a
central vertex with an edge to each of
`n-1` other vertices.
"""
function star_graph(n::Integer)
    g = Graph(n)
    for i in 2:n
        add_edge!(g, 1, i)
    end
    return g
end

"""
    star_cycle_graph(n::Integer)

Return the star graph of order `n`.
The star graph is a cycle graph of order `n-1` with
an additional vertex with and edge to all other vertices.
"""
function star_cycle_graph(n::Integer)
    g = cycle_graph(n - 1)
    add_vertex!(g)
    for i in 1:n-1
        add_edge!(g, i, n)
    end
    return g
end

## Note that Solomonoff and Rappaport introduced this graph and
## computed many important properties 10 years before Erd\"os and Renyi
## studied it.
## The Wikipedia article https://en.wikipedia.org/wiki/Random_graph mentions
## even earlier antecedents.
## Here, we use random_graph rather than a proper name.
"""
    random_graph(n::Integer, p::Real)

Return a `Graph` in which a fraction `p` of all possible edges
is present. This implementation is fast rather than exact.
For instance, it samples with replacement. This algorithm is
rather efficient for large graphs and small `p`.
"""
function random_graph(n::Integer, p::Real)
    p < 1 && p >= 0 || error("The value p=$p is not a probability less than 1")
    g = Graph(n)
    num_edges = round(Int, p * n * (n + 1) / 2)
    @show num_edges
    for i in 1:num_edges
        local v0, v1
        while true
            while true
                v0 = rand(1:n)
                v1 = rand(1:n)
                v0 != v1 && break
            end
            result = add_edge!(g, v0, v1)
            result != nothing && break # Retry if edge already added.
        end
    end
    return g
end


text:
# # SIR model for the spread of COVID-19
# ![](covid_evolution.gif)
# This example illustrates how to use `GraphSpace` and how to model agents on an graph
# (network) where the transition probabilities between each node is not constant.
# ## SIR model

# A SIR model tracks the ratio of Susceptible, Infected, and Recovered individuals within a population.
# Here we add one more category of individuals: those who are infected, but do not know it.
# Transmission rate for infected and diagnosed individuals is lower than infected and undetected.
# We also allow a fraction of recovered individuals to catch the disease again, meaning
# that recovering the disease does not bring full immunity.

# ## Model parameters
# Here are the model parameters, some of which have default values.
# * `Ns`: a vector of population sizes per city. The amount of cities is just `C=length(Ns)`.
# * `β_und`: a vector for transmission probabilities β of the infected but undetected per city.
#   Transmission probability is how many susceptible are infected per day by an infected individual.
#   If social distancing is practiced, this number increases.
# * `β_det`: an array for transmission probabilities β of the infected and detected per city.
#   If hospitals are full, this number increases.
# * `infection_period = 30`: how many days before a person dies or recovers.
# * `detection_time = 14`: how many days before an infected person is detected.
# * `death_rate = 0.02`: the probability that the individual will die after the `infection_period`.
# * `reinfection_probability = 0.05`: The probability that a recovered person can get infected again.
# * `migration_rates`: A matrix of migration probability per individual per day from one city to another.
# * `Is = [zeros(C-1)..., 1]`: An array for initial number of infected but undetected people per city.
#   This starts as only one infected individual in the last city.

# Notice that `Ns, β, Is` all need to have the same length, as they are numbers for each
# city. We've tried to add values to the infection parameters similar to the ones you would hear
# on the news about COVID-19.

# The good thing with Agent based models is that you could easily extend the model we
# implement here to also include age as an additional property of each agent.
# This makes ABMs flexible and suitable for research of virus spreading.

# ## Making the model in Agents.jl
# We start by defining the `PoorSoul` agent type and the ABM
cd(@__DIR__) #src
using Agents, Random, DataFrames, LightGraphs
using Distributions: Poisson, DiscreteNonParametric
using DrWatson: @dict
using Plots
gr() # hide

mutable struct PoorSoul <: AbstractAgent
    id::Int
    pos::Int
    days_infected::Int  # number of days since is infected
    status::Symbol  # 1: S, 2: I, 3:R
end

function model_initiation(;
    Ns,
    migration_rates,
    β_und,
    β_det,
    infection_period = 30,
    reinfection_probability = 0.05,
    detection_time = 14,
    death_rate = 0.02,
    Is = [zeros(Int, length(Ns) - 1)..., 1],
    seed = 0,
)

    Random.seed!(seed)
    @assert length(Ns) ==
    length(Is) ==
    length(β_und) ==
    length(β_det) ==
    size(migration_rates, 1) "length of Ns, Is, and B, and number of rows/columns in migration_rates should be the same "
    @assert size(migration_rates, 1) == size(migration_rates, 2) "migration_rates rates should be a square matrix"

    C = length(Ns)
    ## normalize migration_rates
    migration_rates_sum = sum(migration_rates, dims = 2)
    for c in 1:C
        migration_rates[c, :] ./= migration_rates_sum[c]
    end

    properties = @dict(
        Ns,
        Is,
        β_und,
        β_det,
        β_det,
        migration_rates,
        infection_period,
        infection_period,
        reinfection_probability,
        detection_time,
        C,
        death_rate
    )
    space = GraphSpace(complete_digraph(C))
    model = ABM(PoorSoul, space; properties = properties)

    ## Add initial individuals
    for city in 1:C, n in 1:Ns[city]
        ind = add_agent!(city, model, 0, :S) # Susceptible
    end
    ## add infected individuals
    for city in 1:C
        inds = get_node_contents(city, model)
        for n in 1:Is[city]
            agent = model[inds[n]]
            agent.status = :I # Infected
            agent.days_infected = 1
        end
    end
    return model
end
nothing # hide

# We will make a function that starts a model with `C` number of cities,
# and creates the other parameters automatically by attributing some random
# values to them. You could directly use the above constructor and specify all
# `Ns, β`, etc. for a given set of cities.

# All cities are connected with each other, while it is more probable to travel from a city
# with small population into a city with large population.

using LinearAlgebra: diagind

function create_params(;
    C,
    max_travel_rate,
    infection_period = 30,
    reinfection_probability = 0.05,
    detection_time = 14,
    death_rate = 0.02,
    Is = [zeros(Int, C - 1)..., 1],
    seed = 19,
)

    Random.seed!(seed)
    Ns = rand(50:5000, C)
    β_und = rand(0.3:0.02:0.6, C)
    β_det = β_und ./ 10

    Random.seed!(seed)
    migration_rates = zeros(C, C)
    for c in 1:C
        for c2 in 1:C
            migration_rates[c, c2] = (Ns[c] + Ns[c2]) / Ns[c]
        end
    end
    maxM = maximum(migration_rates)
    migration_rates = (migration_rates .* max_travel_rate) ./ maxM
    migration_rates[diagind(migration_rates)] .= 1.0

    params = @dict(
        Ns,
        β_und,
        β_det,
        migration_rates,
        infection_period,
        reinfection_probability,
        detection_time,
        death_rate,
        Is
    )

    return params
end

params = create_params(C = 8, max_travel_rate = 0.01)
model = model_initiation(; params...)

# Alright, let's plot the cities as a graph to get an idea how the model "looks like",
# using the function [`plotabm`](@ref).

using AgentsPlots

plotargs = (node_size = 0.2, method = :circular, linealpha = 0.4)

plotabm(model; plotargs...)

# The node size is proportional to the relative population of each city.
# In principle we could adjust the edge widths to be proportional with the
# migration rates, by doing:

g = model.space.graph
edgewidthsdict = Dict()
for node in 1:nv(g)
    nbs = neighbors(g, node)
    for nb in nbs
        edgewidthsdict[(node, nb)] = params[:migration_rates][node, nb]
    end
end

edgewidthsf(s, d, w) = edgewidthsdict[(s, d)] * 250

plotargs = merge(plotargs, (edgewidth = edgewidthsf,))

plotabm(model; plotargs...)

# In the following we will be coloring each node according to how large percentage of the
# population is infected. So we create a function to give to [`plotabm`](@ref) as a
# second argument

infected_fraction(x) = cgrad(:inferno)[count(a.status == :I for a in x) / length(x)]
plotabm(model; ac = infected_fraction, plotargs...)

# Here this shows all nodes as black, since we haven't run the model yet. Let's change that!

# ## SIR Stepping functions

# Now we define the functions for modelling the virus spread in time

function agent_step!(agent, model)
    migrate!(agent, model)
    transmit!(agent, model)
    update!(agent, model)
    recover_or_die!(agent, model)
end

function migrate!(agent, model)
    nodeid = agent.pos
    d = DiscreteNonParametric(1:(model.C), model.migration_rates[nodeid, :])
    m = rand(d)
    if m ≠ nodeid
        move_agent!(agent, m, model)
    end
end

function transmit!(agent, model)
    agent.status == :S && return
    rate = if agent.days_infected < model.detection_time
        model.β_und[agent.pos]
    else
        model.β_det[agent.pos]
    end

    d = Poisson(rate)
    n = rand(d)
    n == 0 && return

    for contactID in get_node_contents(agent, model)
        contact = model[contactID]
        if contact.status == :S ||
           (contact.status == :R && rand() ≤ model.reinfection_probability)
            contact.status = :I
            n -= 1
            n == 0 && return
        end
    end
end

update!(agent, model) = agent.status == :I && (agent.days_infected += 1)

function recover_or_die!(agent, model)
    if agent.days_infected ≥ model.infection_period
        if rand() ≤ model.death_rate
            kill_agent!(agent, model)
        else
            agent.status = :R
            agent.days_infected = 0
        end
    end
end
nothing # hide

# ## Example animation
model = model_initiation(; params...)

anim = @animate for i in 0:30
    i > 0 && step!(model, agent_step!, 1)
    p1 = plotabm(model; ac = infected_fraction, plotargs...)
    title!(p1, "Day $(i)")
end

gif(anim, "covid_evolution.gif", fps = 5)

# One can really see "explosive growth" in this animation. Things look quite calm for
# a while and then suddenly supermarkets have no toilet paper anymore!

# ## Exponential growth

# We now run the model and collect data. We define two useful functions for
# data collection:
infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x)
nothing # hide

# and then collect data
model = model_initiation(; params...)

to_collect = [(:status, f) for f in (infected, recovered, length)]
data, _ = run!(model, agent_step!, 100; adata = to_collect)
data[1:10, :]

# We now plot how quantities evolved in time to show
# the exponential growth of the virus

N = sum(model.Ns) # Total initial population
x = data.step
p = plot(
    x,
    log10.(data[:, aggname(:status, infected)]),
    label = "infected",
    xlabel = "steps",
    ylabel = "log(count)",
)
plot!(p, x, log10.(data[:, aggname(:status, recovered)]), label = "recovered")
dead = log10.(N .- data[:, aggname(:status, length)])
plot!(p, x, dead, label = "dead")

# The exponential growth is clearly visible since the logarithm of the number of infected increases
# linearly, until everyone is infected.


text:
using Parallax: main

main(keepalive=!Base.isinteractive())


text:
"""
    RemarkableClient(token::String)

Client used for identification. Simply contains the authentification token
"""
struct RemarkableClient
    token::Ref{String}
    function RemarkableClient(token::String=""; path_to_token::String=pwd())
        if isempty(token) && (
            !isfile(joinpath(path_to_token, ".remarkable_jl_token")) &&
            !isfile(joinpath(path_to_token, ".token"))
        )
            @show
            error(
                """
                    You did not give a token or a path to a ".remarkable_jl_token" (or a deprecated ".token") file to `RemarkableClient`.
                    If you don't know what I am talking about run `register()`.
                    If you have a code, e.g "axervi", run `register(code)` to get a token.
                """,
            )
        elseif isfile(joinpath(path_to_token, ".token"))
            @warn ".token files are deprecated, please rename it to a '.remarkable_jl_token' or rerun `register()`"
            token = String(read(joinpath(path_to_token, ".token")))
        elseif isfile(joinpath(path_to_token, ".remarkable_jl_token"))
            token = String(read(joinpath(path_to_token, ".remarkable_jl_token")))
        end
        client = new(Ref(token))
        refresh_token!(client) # Get a new authentification
        discover_storage(client)
        return client
    end
end

set_token!(c::RemarkableClient, tok::String) = c.token[] = tok
token(c) = c.token[]

Base.show(io::IO, c::RemarkableClient) = print(io, "Remarkable Client")

"""
    refresh_token(client::RemarkableClient) -> token

Before running operations the token needs to be refreshed, `refresh_token!`
does just that!
"""
function refresh_token!(client::RemarkableClient; kwargs...)
    @info "Refreshing the auth token"
    body = HTTP.request(client, "POST", AUTH_API * NEW_TOKEN; kwargs...)
    new_token = String(body)
    set_token!(client, new_token)
    return new_token
end

"""
    discover_storage(client) -> String

Check that the storage URL is still the right one and update it if needed
"""
function discover_storage(client::RemarkableClient; kwargs...)
    @info "Discovering storage host"
    body = HTTP.request(
        client,
        "GET",
        SERVICE_DISCOVERY_API * STORAGE_URL;
        query=Dict(
            "environment" => "production",
            "group" => "auth0%7C5a68dc51cb30df1234567890", # Found from https://akeil.de/posts/remarkable-cloud-api/
            # "group" => "auth0|5a68dc51cb30df3877a1d7c4", # Legacy from RemarkableAPI
            "apiVer" => 2,
        ),
        kwargs...,
    )
    data = JSON.parse(String(body))
    if isempty(data) || data["Status"] != "OK"
        error("Service discovery failed")
    end
    STORAGE_API[] = "https://" * data["Host"]
    return STORAGE_API[]
end

"""
    dict_to_query(d::Dict) -> String

From a dictionary creates a string as :
"key1=value1&key2=value2..."
"""
function dict_to_query(d::Dict)
    query = ""
    for (key, value) in d
        query *= "$(key)=$(value)&"
    end
    return query = query[1:(end - 1)]
end

"""
    request(client, verb, url, headers, body; kwargs...) -> body

Workhorse to communicate with the Remarkable servers.
It is behaving like `HTTP.request`, except that it automatically uses
the identification token in every request.
"""
function HTTP.request(
    client::RemarkableClient,
    verb::String,
    url::String,
    headers::Dict=Dict(),
    body="";
    kwargs...,
)
    merge!(headers, Dict("Authorization" => "Bearer $(token(client))"))
    response = HTTP.request(verb, url, headers, body; kwargs...)
    body = response.body
    # if isempty(body)
    #     error("Returned empty body message")
    # end
    return body
end

"""
    storage_request(client, verb, url, obj, kwargs...)

Specific request related to storage. Converts `obj` to a JSON string,
and automatically transform the body in another JSON.
"""
function storage_request(
    client::RemarkableClient, verb::String, url::String, obj::Dict; kwargs...
)
    body = request_json(client, verb, STORAGE_API[] * ITEM * url, obj; kwargs...)
    data = first(JSON.parse(String(body)))
    if !data["Success"]
        error(data["Message"])
    end
    return data
end

function request_json(
    client::RemarkableClient,
    verb::String,
    url::String,
    data::Dict,
    headers::Dict=Dict();
    kwargs...,
)
    return HTTP.request(client, verb, url, headers, JSON.json([data]); kwargs...)
end


text:
module Defs

using StaticArrays

export unit
function unit end

export bitsign
bitsign(b::Bool) = b ? -1 : 1
bitsign(i::Integer) = bitsign(isodd(i))

export sort_perm
"""
    sort_perm

Sort and count permutations.
"""
sort_perm(xs::SVector{0}) = xs, 0
sort_perm(xs::SVector{1}) = xs, 0
function sort_perm(xs::SVector{D}) where {D}
    @assert D > 1
    xs1 = xs[StaticArrays.SUnitRange(1, D - 1)]
    xend = xs[end]
    rs1, s1 = sort_perm(xs1)
    i = findfirst(>(xend), rs1)
    i === nothing && (i = D)
    rs = SVector{D}(j < i ? rs1[j] : j == i ? xend : rs1[j - 1] for j in 1:D)
    s = s1 + D - i
    # @assert issorted(rs)
    return rs, s
end

function sort_perm(xs::AbstractVector)
    D = length(xs)
    D <= 1 && return xs, 0
    xs1 = @view xs[1:(D - 1)]
    xend = xs[end]
    rs1, s1 = sort_perm(xs1)
    i = findfirst(>(xend), rs1)
    i === nothing && (i = D)
    rs = [j < i ? rs1[j] : j == i ? xend : rs1[j - 1] for j in 1:D]
    s = s1 + D - i
    # @assert issorted(rs)
    return rs, s
end

# Using mergesort, which is probably overkill and slower
# sort_perm(xs::SVector{0}) = xs, false
# sort_perm(xs::SVector{1}) = xs, false
# function sort_perm(xs::SVector{2})
#     if xs[1] <= xs[2]
#         xs, false
#     else
#         reverse(xs), true
#     end
# end
# function sort_perm(xs::SVector{D}) where D
#     @assert D >= 2
#     D1 = prevpow(2, D)
#     @assert 1 <= D1 < D
#     D2 = D - D1
#     @assert 1 <= D2
#     xs1 = xs[StaticArrays.SUnitRange(1, D1)]
#     xs2 = xs[StaticArrays.SUnitRange(D1+1, D2)]
#     rs1, s1 = sort_perm(xs1)
#     rs2, s2 = sort_perm(xs2)
#     i1 = 1
#     i2 = 1
#     rs = xs                     # should be undefined
#     for i in 1:D
#         if rs1[i1] < rs2[i2]
#             rs = setindex(rs, rs1[i1], i)
#             i1 += 1
#             if i1 == D
#                 for j in i+1:D
#                     rs = setindex(rs, rs2[i2], j)
#                     i2 += 1
#                 end
#                 break
#             end
#         else
#             rs = setindex(rs, rs2[i2], i)
#             i2 += 1
#             if i2 == D
#                 for j in i+1:D
#                     rs = setindex(rs, rs1[i1], j)
#                     i1 += 1
#                 end
#                 break
#             end
#         end
#     end
#
#     @assert i1 == D1
#     @assert i2 == D2
#     @assert i == D
#     @assert issorted(rs)
#     rs, s
# end

end


text:
using MolecularGraph: atomsymbol!, atomprop!, atom!

@testset "smarts.atom" begin

@testset "atomsymbol" begin
    state = SmartsParser("Cl", false)
    Chloride = atomsymbol!(state)
    @test Chloride == (:and => (:atomsymbol => :Cl, :isaromatic => false))
    @test state.pos == 3

    state = SmartsParser("Cr", false)
    Chromium = atomsymbol!(state)
    @test Chromium == (:and => (:atomsymbol => :C, :isaromatic => false))
    @test state.pos == 2

    state = SmartsParser("p", false)
    aromp = atomsymbol!(state)
    @test aromp == (:and => (:atomsymbol => :P, :isaromatic => true))
end

@testset "atomprop" begin
    state = SmartsParser("s", false)
    aroms = atomprop!(state)
    @test aroms == (:and => (:atomsymbol => :S, :isaromatic => true))

    state = SmartsParser("123", false)
    iso = atomprop!(state)
    @test iso == (:mass => 123)

    state = SmartsParser("H41", false)
    H4 = atomprop!(state)
    @test H4 == (:hcount => 4)
    @test state.pos == 3

    state = SmartsParser("X2", false)
    X2 = atomprop!(state)
    @test X2 == (:connectivity => 2)

    state = SmartsParser("Xe", false)
    Xe = atomprop!(state)
    @test Xe == (:atomsymbol => :Xe)
    @test state.pos == 3

    state = SmartsParser("Na", false)
    Na = atomprop!(state)
    @test Na == (:atomsymbol => :Na)
    @test state.pos == 3

    state = SmartsParser("Yv2", false)
    Yval = atomprop!(state)
    @test Yval == (:atomsymbol => :Y)

    state = SmartsParser("+23", false)
    chg1 = atomprop!(state)
    @test chg1 == (:charge => 2)
    @test state.pos == 3

    state = SmartsParser("----+", false)
    chg2 = atomprop!(state)
    @test chg2 == (:charge => -4)
    @test state.pos == 5

    state = SmartsParser("@+", false)
    stereo1 = atomprop!(state)
    @test stereo1 == (:stereo => :anticlockwise)
    @test state.pos == 2

    state = SmartsParser("@@?", false)
    stereo4 = atomprop!(state)
    @test stereo4 == (:not => (:stereo => :anticlockwise))
    @test state.pos == 4

    state = SmartsParser("\$([CH2]=*)", false)
    rec = atomprop!(state)
    @test rec == (:recursive => "[CH2]=*")
    @test state.pos == 11
end

@testset "atom" begin
    state = SmilesParser("Br", false)
    br = atom!(state)[1]
    @test br.symbol == :Br

    state = SmilesParser("[2H]", false)
    hyd = atom!(state)
    @test length(hyd) == 1
    @test hyd[1].mass == 2

    state = SmilesParser("[H2]", false)
    hmol = atom!(state)
    @test length(hmol) == 2

    state = SmilesParser("[H+]", false)
    proton = atom!(state)[1]
    @test proton.charge == 1

    state = SmilesParser("[14c@@H]", false)
    iso = atom!(state)
    @test iso[1].symbol == :C
    @test iso[1].mass == 14
    @test iso[1].isaromatic
    @test iso[1].stereo == :clockwise
    @test iso[2].symbol == :H
    @test length(iso) == 2
    @test state.pos == 9

    state = SmilesParser("[NH4+]", false)
    ammonium = atom!(state)
    @test length(ammonium) == 5
    @test ammonium[1].charge == 1
    @test ammonium[5].symbol == :H

    state = SmilesParser("[bH+]", false)
    aromB = atom!(state)
    @test aromB[1].symbol == :B
    @test aromB[1].isaromatic
    @test aromB[1].charge == 1
    @test aromB[2].symbol == :H

    state = SmilesParser("[as]", false)
    aromAs = atom!(state)[1]
    @test aromAs.symbol == :As
    @test aromAs.isaromatic

    state = SmilesParser("[Zn++]", false)
    zn = atom!(state)[1]
    @test zn.symbol == :Zn
    @test zn.charge == 2

    state = SmilesParser("[O-2]", false)
    ox = atom!(state)[1]
    @test ox.symbol == :O
    @test ox.charge == -2

end

@testset "smartsatom" begin
    state = SmartsParser("c", false)
    aromc = atom!(state)[1]
    @test aromc.query == (:and => (:atomsymbol => :C, :isaromatic => true))

    state = SmartsParser("[#16]", false)
    no16 = atom!(state)[1]
    @test no16.query == (:atomsymbol => :S)
end

end # smiles.atom


text:
const ParameterName=Union{Symbol,Expr}

"""
    RandomSearch(bounded=Distributions.Uniform,
                 positive_unbounded=Distributions.Gamma,
                 other=Distributions.Normal,
                 rng=Random.GLOBAL_RNG)

Instantiate a random search tuning strategy, for searching over
Cartesian hyperparameter domains, with customizable priors in each
dimension.

### Supported ranges

A single one-dimensional range or vector of one-dimensioinal ranges
can be specified. If not paired with a prior, then one is fitted,
according to fallback distribution types specified by the tuning
strategy hyperparameters. Specifically, in `RandomSearch`, the `range`
field of a `TunedModel` instance can be:

- a single one-dimensional range (`ParamRange` object) `r`

- a pair of the form `(r, d)`, with `r` as above and where `d` is:

    - a probability vector of the same length as `r.values` (`r` a
      `NominalRange`)

    - any `Distributions.UnivariateDistribution` *instance* (`r` a
      `NumericRange`)

    - one of the *subtypes* of `Distributions.UnivariateDistribution`
      listed in the table below, for automatic fitting using
      `Distributions.fit(d, r)`, a distribution whose support always
      lies between `r.lower` and `r.upper` (`r` a `NumericRange`)

- any pair of the form `(field, s)`, where `field` is the (possibly
  nested) name of a field of the model to be tuned, and `s` an
  arbitrary sampler object for that field. This means only that
  `rand(rng, s)` is defined and returns valid values for the field.

- any vector of objects of the above form

A range vector may contain multiple entries for the same model field,
as in `range = [(:lambda, s1), (:alpha, s), (:lambda, s2)]`. In that
case the entry used in each iteration is random.

distribution types  | for fitting to ranges of this type
--------------------|-----------------------------------
`Arcsine`, `Uniform`, `Biweight`, `Cosine`, `Epanechnikov`, `SymTriangularDist`, `Triweight` | bounded
`Gamma`, `InverseGaussian`, `Poisson` | positive (bounded or unbounded)
`Normal`, `Logistic`, `LogNormal`, `Cauchy`, `Gumbel`, `Laplace`  | any

`ParamRange` objects are constructed using the `range` method.

### Examples

    using Distributions

    range1 = range(model, :hyper1, lower=0, upper=1)

    range2 = [(range(model, :hyper1, lower=1, upper=10), Arcsine),
              range(model, :hyper2, lower=2, upper=Inf, unit=1, origin=3),
              (range(model, :hyper2, lower=2, upper=4), Normal(0, 3)),
              (range(model, :hyper3, values=[:ball, :tree]), [0.3, 0.7])]

    # uniform sampling of :(atom.λ) from [0, 1] without defining a NumericRange:
    struct MySampler end
    Base.rand(rng::Random.AbstractRNG, ::MySampler) = rand(rng)
    range3 = (:(atom.λ), MySampler())

### Algorithm

In each iteration, a model is generated for evaluation by mutating the
fields of a deep copy of `model`. The range vector is shuffled and the
fields sampled according to the new order (repeated fields being
mutated more than once). For a `range` entry of the form `(field, s)`
the algorithm calls `rand(rng, s)` and mutates the field `field` of
the model clone to have this value. For an entry of the form `(r, d)`,
`s` is substituted with `sampler(r, d)`. If no `d` is specified, then
sampling is uniform (with replacement) if `r` is a `NominalRange`, and
is otherwise given by the defaults specified by the tuning strategy
parameters `bounded`, `positive_unbounded`, and `other`, depending on
the field values of the `NumericRange` object `r`.

See also [`TunedModel`](@ref), [`range`](@ref), [`sampler`](@ref).

"""
mutable struct RandomSearch <: TuningStrategy
    bounded
    positive_unbounded
    other
    rng::Random.AbstractRNG
end

# Constructor with keywords
function RandomSearch(; bounded=Distributions.Uniform,
                      positive_unbounded=Distributions.Gamma,
                      other=Distributions.Normal,
                      rng=Random.GLOBAL_RNG)
    (bounded isa Type{<:Distributions.UnivariateDistribution} &&
        positive_unbounded isa Type{<:Distributions.UnivariateDistribution} &&
        other isa Type{<:Distributions.UnivariateDistribution}) ||
        error("`bounded`, `positive_unbounded` and `other` "*
              "must all be subtypes of "*
              "`Distributions.UnivariateDistribution`. ")

    _rng = rng isa Integer ? Random.MersenneTwister(rng) : rng
    return RandomSearch(bounded, positive_unbounded, other, _rng)
end

# `state` consists of a tuple of (field, sampler) pairs (that gets
# shuffled each iteration):
setup(tuning::RandomSearch, model, user_range, n, verbosity) =
    process_random_range(user_range,
                              tuning.bounded,
                              tuning.positive_unbounded,
                              tuning.other) |> collect

function MLJTuning.models(tuning::RandomSearch,
                          model,
                          history,
                          state, # tuple of (field, sampler) pairs
                          n_remaining,
                          verbosity)
    new_models = map(1:n_remaining) do _
        clone = deepcopy(model)
        Random.shuffle!(tuning.rng, state)
        for (fld, s) in state
            recursive_setproperty!(clone, fld, rand(tuning.rng, s))
        end
        clone
    end
    return new_models, state
end

function tuning_report(tuning::RandomSearch, history, field_sampler_pairs)

    fields = first.(field_sampler_pairs)
    parameter_scales = map(field_sampler_pairs) do (fld, s)
        scale(s)
    end

    plotting = plotting_report(fields, parameter_scales, history)

    return (plotting=plotting,)

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type tank_simplified
	tank_simplified()=begin
		new(
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Valve Constant",
				:Unit=>"m^2.5/h",
				:Default=>4
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Tank area",
				:Default=>2
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tank level"
			)),
			flow_vol(Dict{Symbol,Any}(
				:Brief=>"Input flow",
				:PosX=>0.3037,
				:PosY=>0
			)),
			flow_vol(Dict{Symbol,Any}(
				:Brief=>"Output flow",
				:PosX=>1,
				:PosY=>1
			)),
			[
				:(diff(A*Level) = Fin - Fout),
				:(Fout = k*sqrt(Level)),
			],
			[
				"Mass balance","Valve equation",
			],
			[:k,:A,],
			[:Level,:Fin,:Fout,]
		)
	end
	k::DanaReal
	A::area
	Level::length
	Fin::flow_vol
	Fout::flow_vol
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export tank_simplified
function setEquationFlow(in::tank_simplified)
	addEquation(1)
	addEquation(2)
end
function atributes(in::tank_simplified,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Tank"
	fields[:Brief]="Model of a simplified tank."
	fields[:Info]="== Specify ==
* the Inlet flow rate;

== Initial Conditions ==
* the tank initial level (Level);
"
	drive!(fields,_)
	return fields
end
tank_simplified(_::Dict{Symbol,Any})=begin
	newModel=tank_simplified()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
export GaussianMeanPrecision

"""
Description:

    A Gaussian with mean-precision parameterization:

    f(out,m,w) = 𝒩(out|m,w) = (2π)^{-D/2} |w|^{1/2} exp(-1/2 (out - m)' w (out - m))

Interfaces:

    1. out
    2. m (mean)
    3. w (precision)

Construction:

    GaussianMeanPrecision(out, m, w, id=:some_id)
"""
mutable struct GaussianMeanPrecision <: Gaussian
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function GaussianMeanPrecision(out, m, w; id=generateId(GaussianMeanPrecision))
        @ensureVariables(out, m, w)
        self = new(id, Array{Interface}(undef, 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:m] = self.interfaces[2] = associate!(Interface(self), m)
        self.i[:w] = self.interfaces[3] = associate!(Interface(self), w)

        return self
    end
end

slug(::Type{GaussianMeanPrecision}) = "𝒩"

format(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = "$(slug(GaussianMeanPrecision))(m=$(format(dist.params[:m])), w=$(format(dist.params[:w])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{GaussianMeanPrecision}; m=0.0, w=1.0) = ProbabilityDistribution{Univariate, GaussianMeanPrecision}(Dict(:m=>m, :w=>w))
ProbabilityDistribution(::Type{GaussianMeanPrecision}; m::Number=0.0, w::Number=1.0) = ProbabilityDistribution{Univariate, GaussianMeanPrecision}(Dict(:m=>m, :w=>w))
ProbabilityDistribution(::Type{Multivariate}, ::Type{GaussianMeanPrecision}; m=[0.0], w=transpose([1.0])) = ProbabilityDistribution{Multivariate, GaussianMeanPrecision}(Dict(:m=>m, :w=>w))

dims(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = length(dist.params[:m])

vague(::Type{GaussianMeanPrecision}) = ProbabilityDistribution(Univariate, GaussianMeanPrecision, m=0.0, w=tiny)
vague(::Type{GaussianMeanPrecision}, dims::Int64) = ProbabilityDistribution(Multivariate, GaussianMeanPrecision, m=zeros(dims), w=tiny*diageye(dims))
vague(::Type{GaussianMeanPrecision}, dims::Tuple{Int64}) = ProbabilityDistribution(Multivariate, GaussianMeanPrecision, m=zeros(dims), w=tiny*diageye(dims[1]))

unsafeMean(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = deepcopy(dist.params[:m]) # unsafe mean

unsafeMode(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = deepcopy(dist.params[:m]) # unsafe mode

unsafeVar(dist::ProbabilityDistribution{Univariate, GaussianMeanPrecision}) = 1.0/dist.params[:w] # unsafe variance
unsafeVar(dist::ProbabilityDistribution{Multivariate, GaussianMeanPrecision}) = diag(cholinv(dist.params[:w]))

unsafeCov(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = cholinv(dist.params[:w])

unsafeMeanCov(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = (deepcopy(dist.params[:m]), cholinv(dist.params[:w]))

unsafeWeightedMean(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = dist.params[:w]*dist.params[:m] # unsafe weighted mean

unsafePrecision(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = deepcopy(dist.params[:w]) # unsafe precision

unsafeWeightedMeanPrecision(dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType = (dist.params[:w]*dist.params[:m], deepcopy(dist.params[:w]))

logPdf(dist::ProbabilityDistribution{Univariate, GaussianMeanPrecision}, x) = -0.5*(log(2pi) - log(dist.params[:w]) + (x-dist.params[:m])^2*dist.params[:w])
logPdf(dist::ProbabilityDistribution{Multivariate, GaussianMeanPrecision}, x) = -0.5*(dims(dist)*log(2pi) - logdet(dist.params[:w]) + transpose(x-dist.params[:m])*dist.params[:w]*(x-dist.params[:m]))

isProper(dist::ProbabilityDistribution{Univariate, GaussianMeanPrecision}) = (floatmin(Float64) < dist.params[:w] < floatmax(Float64))
isProper(dist::ProbabilityDistribution{Multivariate, GaussianMeanPrecision}) = isRoundedPosDef(dist.params[:w])

function ==(t::ProbabilityDistribution{V, GaussianMeanPrecision}, u::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType
    (t === u) && return true
    isApproxEqual(t.params[:m], u.params[:m]) && isApproxEqual(t.params[:w], u.params[:w])
end

function ==(t::ProbabilityDistribution{V, GaussianMeanPrecision}, u::ProbabilityDistribution{V, F}) where {V<:VariateType, F<:Gaussian}
    (t === u) && return true
    isApproxEqual(t.params[:m], unsafeMean(u)) && isApproxEqual(t.params[:w], unsafePrecision(u))
end

# Average energy functional
function averageEnergy(::Type{GaussianMeanPrecision}, marg_out::ProbabilityDistribution{Univariate}, marg_mean::ProbabilityDistribution{Univariate}, marg_prec::ProbabilityDistribution{Univariate})
    (m_mean, v_mean) = unsafeMeanCov(marg_mean)
    (m_out, v_out) = unsafeMeanCov(marg_out)

    0.5*log(2*pi) -
    0.5*unsafeLogMean(marg_prec) +
    0.5*unsafeMean(marg_prec)*(v_out + v_mean + (m_out - m_mean)^2)
end

function averageEnergy(::Type{GaussianMeanPrecision}, marg_out::ProbabilityDistribution{Multivariate}, marg_mean::ProbabilityDistribution{Multivariate}, marg_prec::ProbabilityDistribution{MatrixVariate})
    (m_mean, v_mean) = unsafeMeanCov(marg_mean)
    (m_out, v_out) = unsafeMeanCov(marg_out)

    0.5*dims(marg_out)*log(2*pi) -
    0.5*unsafeDetLogMean(marg_prec) +
    0.5*tr( unsafeMean(marg_prec)*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)' ))
end

function averageEnergy(::Type{GaussianMeanPrecision}, marg_out_mean::ProbabilityDistribution{Multivariate, F}, marg_prec::ProbabilityDistribution{Univariate}) where F<:Gaussian
    (m, V) = unsafeMeanCov(marg_out_mean)

    0.5*log(2*pi) -
    0.5*unsafeLogMean(marg_prec) +
    0.5*unsafeMean(marg_prec)*( V[1,1] - V[1,2] - V[2,1] + V[2,2] + (m[1] - m[2])^2 )
end

function averageEnergy(::Type{GaussianMeanPrecision}, marg_out_mean::ProbabilityDistribution{Multivariate, F}, marg_prec::ProbabilityDistribution{MatrixVariate}) where F<:Gaussian
    (m, V) = unsafeMeanCov(marg_out_mean)
    d = Int64(dims(marg_out_mean)/2)

    0.5*d*log(2*pi) -
    0.5*unsafeDetLogMean(marg_prec) +
    0.5*tr( unsafeMean(marg_prec)*( V[1:d,1:d] - V[1:d,d+1:end] - V[d+1:end,1:d] + V[d+1:end,d+1:end] + (m[1:d] - m[d+1:end])*(m[1:d] - m[d+1:end])' ) )
end


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

JLLWrappers.@generate_main_file_header("Gettext")
JLLWrappers.@generate_main_file("Gettext", UUID("78b55507-aeef-58d4-861c-77aaff3498b1"))
end  # module Gettext_jll


text:
export DoubleIntR2
export init_traj_straightline, init_traj_geometricplan

mutable struct DoubleIntR2 <: DynamicsModel
  # state: r v p omega
  x_dim
  u_dim
  clearance

  # Parameters that can be updated
  f::Vector
  A
  B
end

function DoubleIntR2()
  x_dim = 4
  u_dim = 2
  clearance = 0.05
  DoubleIntR2(x_dim, u_dim, clearance, [], [], [])
end

function SCPParam(model::DoubleIntR2, fixed_final_time::Bool)
  convergence_threshold = 0.01

  SCPParam(fixed_final_time, convergence_threshold)
end

function SCPParam_GuSTO(model::DoubleIntR2)
  Δ0 = 1.
  ω0 = 1.
  ω_max = 1.0e10
  ε = 1.0e-6
  ρ0 = 0.01
  ρ1 = 0.3
  β_succ = 2.
  β_fail = 0.5
  γ_fail = 10.

  SCPParam_GuSTO(Δ0, ω0, ω_max, ε, ρ0, ρ1, β_succ, β_fail, γ_fail)
end

###############
# Gurobi stuff
###############
function cost_true(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2}) where T
  U, N = traj.U, SCPP.N
  dtp = traj_prev.dt
  Jm = 0
  for k in 1:N-1
    Jm += dtp*norm(U[:,k])^2
  end
  return Jm
end

function cost_true_convexified(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}) where {T,E}
  cost_true(traj, traj_prev, SCPP)
end

#############################
# Trajectory Initializations
#############################
function init_traj_straightline(TOP::TrajectoryOptimizationProblem{PointMassInSphere{T}, DoubleIntR2, E}) where {T,E}
  model, x_init, x_goal = TOP.PD.model, TOP.PD.x_init, TOP.PD.x_goal
  x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess
  N = TOP.N

  X = hcat(range(x_init, stop=x_goal, length=N)...)
  U = zeros(u_dim, N)
  Trajectory(X, U, tf_guess)
end

####################
# Constraint-adding
####################
function initialize_model_params!(SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, traj_prev::Trajectory) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  x_dim, u_dim = model.x_dim, model.u_dim
  Xp, Up = traj_prev.X, traj_prev.U

  model.f, model.A, model.B = [], A_dyn(Xp[:,1],robot,model), B_dyn(Xp[:,1],robot,model)
  for k = 1:N-1
    push!(model.f, f_dyn(Xp[:,k],Up[:,k],robot,model))
  end
end

function update_model_params!(SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, traj_prev::Trajectory) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  Xp, Up, f = traj_prev.X, traj_prev.U, model.f

  for k = 1:N-1
    update_f!(f[k], Xp[:,k], Up[:,k], robot, model)
  end
end

macro constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)
  quote
    X, U, Tf = $(esc(traj)).X, $(esc(traj)).U, $(esc(traj)).Tf
    Xp, Up, Tfp, dtp = $(esc(traj_prev)).X, $(esc(traj_prev)).U, $(esc(traj_prev)).Tf, $(esc(traj_prev)).dt
    robot, model, WS, x_init, x_goal = $(esc(SCPP)).PD.robot, $(esc(SCPP)).PD.model, $(esc(SCPP)).WS, $(esc(SCPP)).PD.x_init, $(esc(SCPP)).PD.x_goal
    x_dim, u_dim, N, dh = model.x_dim, model.u_dim, $(esc(SCPP)).N, $(esc(SCPP)).dh
    X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh
  end
end

## Dynamics constraints
function dynamics_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, k::Int, i::Int) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)

  return A_dyn_discrete(X[:,k],dtp,robot,model)*X[:,k] + B_dyn_discrete(X[:,k],dtp,robot,model)*U[:,k] - X[:,k+1]
end

# Get current dynamics structures for a time step
get_f(k::Int, model::DoubleIntR2) = model.f[k]
get_A(k::Int, model::DoubleIntR2) = model.A
get_B(k::Int, model::DoubleIntR2) = model.B

# Generate full dynamics structures for a time step
function f_dyn(x::Vector, u::Vector, robot::Robot, model::DoubleIntR2)
  x_dim = model.x_dim
  f = zeros(x_dim)
  update_f!(f, x, u, robot, model)
  return f
end

function update_f!(f, x::Vector, u::Vector, robot::Robot, model::DoubleIntR2)
  f[3:4] = 1/robot.mass*u
end

function A_dyn(x::Vector, robot::Robot, model::DoubleIntR2)
  kron([0 1; 0 0], Eye(2))
end

function B_dyn(x::Vector, robot::Robot, model::DoubleIntR2)
  B = zeros(4,2)
  B[3:4,1:2] = Eye(2)/robot.mass
  return B
end

# Generate full discrete update version of dynamics matrices for a time step
# TODO(ambyld): Rename these? Initialize these once?
function A_dyn_discrete(x, dt, robot::Robot, model::DoubleIntR2)
  kron([1 dt; 0 1], Eye(2))
end

function B_dyn_discrete(x, dt, robot::Robot, model::DoubleIntR2)
  [0.5*dt^2*Eye(2);
   dt*Eye(2)]/robot.mass
end

## Convex state inequality constraints
function csi_translational_velocity_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)
  return norm(X[3:4,k]) - robot.hard_limit_vel
end
## Convex control inequality constraints
function cci_translational_accel_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)
  return 1/robot.mass*norm(U[1:2,k]) - robot.hard_limit_accel
end

## Nonconvex state inequality constraints
function ncsi_obstacle_avoidance_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)

  rb_idx, env_idx = 1, i
  env_ = WS.btenvironment_keepout

  clearance = model.clearance

  r = get_workspace_location(traj, SCPP, k)
  dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)

  # See Eq. 12b in "Convex optimization for proximity maneuvering of a spacecraft with a robotic manipulator"
  return clearance - dist
end

## Nonconvex state inequality constraints (convexified)
function ncsi_body_obstacle_avoidance_constraints_convexified(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)

  env_ = WS.btenvironment_keepout
  env_idx = i
  clearance = model.clearance

  # Base
  rb_idx = 1
  r0 = get_workspace_location(traj_prev, SCPP, k)
  dist, xbody, xobs = BulletCollision.distance(env_, rb_idx, r0, env_idx)
  r = get_workspace_location(traj, SCPP, k)

  if dist < SCPP.param.obstacle_toggle_distance
    # See Eq. 12b in "Convex optimization for proximity maneuvering of a spacecraft with a robotic manipulator"
    nhat = dist > 0 ?
      (xbody-xobs)./norm(xbody-xobs) :
      (xobs-xbody)./norm(xobs-xbody)

    return clearance - (dist + nhat'*(r-r0))
  else
    return 0.
  end
end

## State trust region inequality constraints
function stri_state_trust_region(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)
  return norm(X[:,k]-Xp[:,k])^2
end

## Convex state inequality constraints
function ctri_control_trust_region(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)
  return norm(U[:,k]-Up[:,k])
end

function get_workspace_location(traj, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, k::Int, i::Int=0) where {T,E}
  return [traj.X[1:2,k]; 0]
end

## Constructing full list of constraint functions
function SCPConstraints(SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}) where {T,E}
  model = SCPP.PD.model
  x_dim, u_dim, N = model.x_dim, model.u_dim, SCPP.N
  WS = SCPP.WS

  SCPC = SCPConstraints()

  ## Dynamics constraints
  for k = 1:N-1
    push!(SCPC.dynamics, (dynamics_constraints, k, 0))
  end

  ## Convex state equality constraints
  # Init and goal (add init first for convenience of getting dual)
  for i = 1:x_dim
    push!(SCPC.convex_state_eq, (cse_init_constraints, 0, i))
  end
  for i = 1:x_dim
    push!(SCPC.convex_state_eq, (cse_goal_constraints, 0, i))
  end

  ## Convex state inequality constraints
  for k = 1:N
    push!(SCPC.convex_state_ineq, (csi_translational_velocity_bound, k, 0))
  end

  ## Nonconvex state equality constraints (convexified)
  nothing

  ## Nonconvex state inequality constraints (convexified)
  env_ = WS.btenvironment_keepout
  for k = 1:N, i = 1:length(env_.convex_env_components)
    push!(SCPC.nonconvex_state_convexified_ineq, (ncsi_body_obstacle_avoidance_constraints_convexified, k, i))
  end

  ## Convex control equality constraints
  nothing

  ## Convex control inequality constraints
  for k = 1:N-1
    push!(SCPC.convex_control_ineq, (cci_translational_accel_bound, k, 0))
  end

  ## State trust region ineqality constraints
  for k = 1:N
    push!(SCPC.state_trust_region_ineq, (stri_state_trust_region, k, 0))
  end

  ## Constrol trust region inequality constraints
  for k = 1:N-1
    push!(SCPC.control_trust_region_ineq, (ctri_control_trust_region, k, 0))
  end

  return SCPC
end

# TODO: Generalize this? Need to make A always a vector
function trust_region_ratio_gusto(traj, traj_prev::Trajectory, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_DoubleIntR2(traj, traj_prev, SCPP)
  fp, Ap = model.f, model.A
  num,den = 0, 0
  env_ = WS.btenvironment_keepout

  for k in 1:N-1
    linearized = fp[k] + Ap*(X[:,k]-Xp[:,k])
    num += norm(f_dyn(X[:,k],U[:,k],robot,model) - linearized)
    den += norm(linearized)
  end

  clearance = model.clearance

  # TODO(ambyld): Update for arm
  for k in 1:N
    r0 = get_workspace_location(traj, SCPP, k)
    r = get_workspace_location(traj_prev, SCPP, k)
    for (rb_idx, body_point) in enumerate(env_.convex_robot_components)
      for (env_idx,convex_env_component) in enumerate(env_.convex_env_components)
        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r0,env_idx)
        nhat = dist > 0 ?
          (xbody-xobs)./norm(xbody-xobs) :
          (xobs-xbody)./norm(xobs-xbody)
        linearized = clearance - (dist + nhat'*(r-r0))

        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)

        num += abs((clearance-dist) - linearized)
        den += abs(linearized)
      end
    end
  end
  return num/den
end

function get_dual_cvx(prob::Convex.Problem, SCPP::SCPProblem{PointMassInSphere{T}, DoubleIntR2, E}, solver) where {T,E}
  if solver == "Mosek"
    return -MathProgBase.getdual(prob.model)[1:SCPP.PD.model.x_dim]
  else
    return []
  end
end


text:
module MemTransportTests

using Thrift
using Test

function testmemtransport()
    println("\nTesting memory transport...")
    t = TMemoryTransport()
    p = TCompactProtocol(t)

    @test open(t) === nothing
    @test close(t) === nothing
    @test flush(t) === nothing
    @test isa(Thrift.rawio(t), IOBuffer)
    @test isopen(t)

    s1 = "Hello World"
    write(p, s1)
    s2 = read(p, typeof(s1))
    @test s2 == s1

    t = TMemoryTransport()
    p = TCompactProtocol(t)
    write(p, s1)
    t = TMemoryTransport(take!(t.buff))
    p = TCompactProtocol(t)
    s2 = read(p, typeof(s1))
    @test s2 == s1

    println("passed.")
end

testmemtransport()

end


text: const StaticTypes = DomainSets.StaticTypes


text:
using DelimitedFiles
import Base.size
import Base.getindex
import Base.setindex!
import Base.eachindex
import Base.CartesianIndices

function readtensor(filename)
    data = readdlm(filename, Float64)
    cartesianindices = []

    size_tensor = Tuple([Int(maximum(data[:,i])) for i in 1:length(data[1,:])-1])
    densetensor = Array{Union{Missing, Float64}}(missing, size_tensor)
    for i in 1:size(data)[1]
        I = CartesianIndex(Tuple(convert(Array{Int}, data[i,1:end-1])))
        push!(cartesianindices, I)
        densetensor[I] = data[i,end]
    end
    return densetensor, cartesianindices
end
function fill(T, inds, fill)
    if fill == "zeros"
        filled = zeros(size(T))
        for I in inds
            filled[I] = T[I]
        end
        return filled
    else
        println("fill not implemented")
    end
end

"""
struct SparseTensor{T,N} <: AbstractArray{T,N}
    # T is set to Float64. Should be ...
    elements_::Dict
    size_
    function SparseTensor(filename::String)
        elements = Dict()
        a = readdlm(filename, Float64)
        for i in 1:size(a)[1]
            index = CartesianIndex(Tuple(convert(Array{Int}, a[i,1:end-1])))
            elements[index] = a[i,end]
        end
        @assert length(keys(elements)) == length(unique(keys(elements)))
        s = Tuple(maximum(keys(elements)))
        new{Float64, length(a[1,:])-1}(elements, s)
    end
    function SparseTensor(elements::Dict)
        s = Tuple(maximum(keys(T.elements_)))
        new{Float64, length(collect(keys(elements))[1])}(elements, s)
    end
    function SparseTensor(elements::Dict, s::Tuple)
        new{Float64, length(collect(keys(elements))[1])}(elements, s)
    end
    function SparseTensor(s::Tuple)
        new{Float64, length(s)}(Dict(), s)
    end
end
function SparseTensor(old::SparseTensor, selection::Array)
    new = SparseTensor(size(old))
    for I in selection
        new[I] = old[I]
    end
    return new
end
function SparseTensor(a::Array)
    d = Dict()
    for I in CartesianIndices(a)
        d[I] = a[I]
    end
    return SparseTensor(d, size(a))
end
function save(T::SparseTensor, name)
    tosave = zeros(length(keys(T.elements_)), length(size(T))+1)
    for it in zip(1:length(keys(T.elements_)),keys(T.elements_))
        i, key  = it
        for j in 1:length(size(T))
            tosave[i,j] = key[j]
        end
        tosave[i,end] = T.elements_[key]
    end
    writedlm(name, tosave)
end
function size(T::SparseTensor)
    return T.size_
end
function getindex(T::SparseTensor, ind::CartesianIndex)
    for i in 1:length(ind)
        if ind[i]>T.size_[i]
            throw("index out of range")
        end
    end
    try
        return T.elements_[ind]
    catch e
        return missing
    end
end
function getindex(T::SparseTensor, ind...)
    ind = CartesianIndex(ind)
    for i in 1:length(ind)
        if ind[i]>T.size_[i]
            throw("index out of range")
        end
    end
    try
        return T.elements_[ind]
    catch e
        return missing
    end
end
function setindex!(T::SparseTensor, v, ind::CartesianIndex)
    for i in 1:length(ind)
        if ind[i]>T.size_[i]
            throw("index out of range")
        end
    end
    T.elements_[ind] = v
end
function setindex!(T::SparseTensor, v, ind...)
    ind = CartesianIndex(ind)
    for i in 1:length(ind)
        if ind[i]>T.size_[i]
            throw("index out of range")
        end
    end
    T.elements_[ind] = v
end
function eachindex(T::SparseTensor)
    return collect(keys(T.elements_))
end

function getindex(T::SparseTensor, inds::Array)
    res = Array{Float64}(undef, length(inds))
    for i in 1:length(inds)
        res[i] = T[inds[i]]
    end
    return res
end

function density(T::SparseTensor)
    return length(T.elements_) / prod(size(T))
end
function full(T::SparseTensor)
    if prod(size(T)) == length(T.elements_)
        return true
    elseif prod(size(T)) > length(T.elements_)
        return false
    else
        throw("should not be possible: more indices than size allows")
    end
end
function tensor(T::SparseTensor)
    @assert full(T)
    return tensor(T, "zeros")
end
function tensor(T::SparseTensor, howtofill::String)
    if howtofill=="zeros"
        tensor = zeros(size(T))
        for I in keys(T.elements_)
            tensor[I] = T.elements_[I]
        end
        return tensor
    else
        println("this fill method not implemented:  "* how)
    end
end
"""

max_stars_repo_name: padwulf/MyExample

max_issues_repo_name: padwulf/MyExample

max_forks_repo_name: padwulf/MyExample


text:
import HTTP

#************************ LOCALS TO CONFIGURE!!!! **************************
const email= "" #This is an enviroment variable that you need to setup
const search_term="(obesity[MeSH Major Topic]) AND (\"2010\"[Date - Publication] : \"2012\"[Date - Publication])"
const max_articles = 2
const overwrite_db=true
const verbose = false
const umls_user = get(ENV, "UMLS_USER", "")
const umls_pswd = get(ENV, "UMLS_PSSWD", "")
#************************ SQLite **************************
const db_path="./test_processes.sqlite"
#***************************************************************************

if isfile(db_path)
    rm(db_path)
end

const conn_sql = SQLite.DB(db_path)
PubMed.create_tables!(conn_sql)

@testset "Save and Search" begin

    println("-----------------------------------------")
    println("       Testing Search and Save")

    pubmed_search_and_save!(email, search_term, max_articles,
    conn_sql, verbose)
    #query the article table and make sure the count is correct
    all_pmids = BioMedQuery.PubMed.all_pmids(conn_sql)
    @test length(all_pmids) == max_articles

end

global credentials_set = get(ENV, "TRAVIS_SECURE_ENV_VARS", "true")=="true"

@testset "UMLS" begin

    if credentials_set
        println("-----------------------------------------")
        println("       Testing MESH2UMLS")

        append = false

        map_mesh_to_umls_async!(conn_sql, umls_user, umls_pswd; append_results=append, timeout=1)

        all_pairs_query = db_query(conn_sql, "SELECT mesh FROM mesh2umls;")
        all_pairs = all_pairs_query[1]
        @test length(all_pairs) > 0

        println("-----------------------------------------")
        println("       Testing Occurrences")
        umls_concept = "Disease or Syndrome"

        labels2ind, occur = umls_semantic_occurrences(conn_sql, umls_concept)

        @test length(keys(labels2ind)) > 0
        @test length(findall(x->x=="Obesity", collect(keys(labels2ind)))) ==1
    end

end

# remove temp files
if isfile(db_path)
    rm(db_path)
end
println("------------End Test Processes SQLite-----------")
println("------------------------------------------------")


text:
using FiniteDifferenceSchemes
using Test

@testset "LinearSchemes" begin
    upwind_3 = linear_scheme(-1, 1)
    vector_in = [1., 2., 3.]
    flux_negative_linear = flux_negative_function(flux_positive_linear,
                                                  upwind_3)
    @test flux_positive_linear(upwind_3, vector_in) ≈ 2.5
    @test flux_negative_linear(upwind_3, vector_in) ≈ 1.5
end

@testset "Weno5-JS" begin
    weno5 = weno5_js(1.e-8)
    vector_in = [1., 2., 3., 4., 5.]
    flux_negative_weno5 = flux_negative_function(flux_positive_weno5,
                                                 weno5)
    @test flux_positive_weno5(weno5, vector_in) ≈ 3.5
    @test flux_negative_weno5(weno5, vector_in) ≈ 2.5
end


text:
using DiffEqDiffTools, LinearAlgebra, SparseArrays, Test, LinearAlgebra, BlockBandedMatrices, ArrayInterface, BandedMatrices

fcalls = 0
function f(dx,x)
  global fcalls += 1
  for i in 2:length(x)-1
    dx[i] = x[i-1] - 2x[i] + x[i+1]
  end
  dx[1] = -2x[1] + x[2]
  dx[end] = x[end-1] - 2x[end]
  nothing
end
function oopf(x)
  global fcalls += 1
  vcat([-2x[1]+x[2]],[x[i-1]-2x[i]+x[i+1] for i in 2:length(x)-1],[x[end-1]-2x[end]])
end

function second_derivative_stencil(N)
  A = zeros(N,N)
  for i in 1:N, j in 1:N
      (j-i==-1 || j-i==1) && (A[i,j]=1)
      j-i==0 && (A[i,j]=-2)
  end
  A
end

J = DiffEqDiffTools.finite_difference_jacobian(oopf,rand(30))
@test J ≈ second_derivative_stencil(30)
_J = sparse(J)
@test fcalls == 31

_J2 = similar(_J)
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_J2,f,rand(30),colorvec=repeat(1:3,10))
@test fcalls == 4
@test _J2 ≈ _J

_J2 = similar(_J)
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_J2,f,rand(30),Val{:central},colorvec=repeat(1:3,10))
@test fcalls == 6
@test _J2 ≈ _J

_J2 = similar(_J)
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_J2,f,rand(30),Val{:complex},colorvec=repeat(1:3,10))
@test fcalls == 3
@test _J2 ≈ _J

_J2 = similar(_J)
_denseJ2 = collect(_J2)
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_denseJ2,f,rand(30),colorvec=repeat(1:3,10),sparsity=_J2)
@test fcalls == 4
@test sparse(_denseJ2) ≈ _J

_J2 = similar(_J)
_denseJ2 = collect(_J2)
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_denseJ2,f,rand(30),Val{:central},colorvec=repeat(1:3,10),sparsity=_J2)
@test fcalls == 6
@test sparse(_denseJ2) ≈ _J

_J2 = similar(_J)
_denseJ2 = collect(_J2)
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_denseJ2,f,rand(30),Val{:complex},colorvec=repeat(1:3,10),sparsity=_J2)
@test fcalls == 3
@test sparse(_denseJ2) ≈ _J

_J2 = similar(Tridiagonal(_J))
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_J2,f,rand(30),colorvec=repeat(1:3,10))
@test fcalls == 4
@test _J2 ≈ _J

_J2 = similar(_J2)
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_J2,f,rand(30),Val{:central},colorvec=repeat(1:3,10))
@test fcalls == 6
@test _J2 ≈ _J

_J2 = similar(_J2)
fcalls = 0
DiffEqDiffTools.finite_difference_jacobian!(_J2,f,rand(30),Val{:complex},colorvec=repeat(1:3,10))
@test fcalls == 3
@test _J2 ≈ _J

_Jb = BandedMatrices.BandedMatrix(similar(_J2),(1,1))
DiffEqDiffTools.finite_difference_jacobian!(_Jb, f, rand(30), colorvec=repeat(1:3,10))
@test _Jb ≈ _J

_Jtri = Tridiagonal(similar(_J2))
DiffEqDiffTools.finite_difference_jacobian!(_Jtri, f, rand(30), colorvec=repeat(1:3,10))
@test _Jtri ≈ _J

#https://github.com/JuliaDiffEq/DiffEqDiffTools.jl/issues/67#issuecomment-516871956
function f(out, x)
	x = reshape(x, 100, 100)
	out = reshape(out, 100, 100)
	for i in 1:100
		for j in 1:100
			out[i, j] = x[i, j] + x[max(i -1, 1), j] + x[min(i+1, size(x, 1)), j] +  x[i, max(j-1, 1)]  + x[i, min(j+1, size(x, 2))]
		end
	end
	return vec(out)
end
x = rand(10000)
Jbbb = BandedBlockBandedMatrix(Ones(10000, 10000), (fill(100, 100), fill(100, 100)), (1, 1), (1, 1))
Jsparse = sparse(Jbbb)
colorsbbb = ArrayInterface.matrix_colors(Jbbb)
DiffEqDiffTools.finite_difference_jacobian!(Jbbb, f, x, colorvec=colorsbbb)
DiffEqDiffTools.finite_difference_jacobian!(Jsparse, f, x, colorvec=colorsbbb)
@test Jbbb ≈ Jsparse
Jbb = BlockBandedMatrix(similar(Jsparse),(fill(100, 100), fill(100, 100)),(1,1));
colorsbb = ArrayInterface.matrix_colors(Jbb)
DiffEqDiffTools.finite_difference_jacobian!(Jbb, f, x, colorvec=colorsbb)
@test Jbb ≈ Jsparse


text:
function getctnf(netfile, outfile)
    a = readdlm(netfile, Int)
    m = Dict{Int,Int}()
    g = Graph()
    fout = open(outfile, "w")
    if size(a,2) < 3
        for i=1:size(a,1)
            x=copy(m)
            lpa_addedge!(g, m, a[i,1], a[i,2])
            println(fout, '#', i)
            println(fout, "+\t", a[i,1], '\t', a[i,2])
            commdiff(x, m, fout)
        end
    else
        for i=1:size(a,1)
            x=copy(m)
            println(fout, '#', i)
            if a[i,1] > 0
                lpa_addedge!(g, m, a[i,2], a[i,3])
                println(fout, "+\t", a[i,2], '\t', a[i,3])
            else
                lpa_deleteedge!(g, m, a[i,2], a[i,3])
                println(fout, "-\t", a[i,2], '\t', a[i,3])
            end
            commdiff(x, m, fout)
        end
    end
    close(fout)
    convertfile(outfile, "converted_$(outfile)")
end

function commdiff(m1, m2, fout)
    x = Dict{Int,Set{Int}}()
    y = Dict{Int,Set{Int}}()
    for (k,v) in m1
        if haskey(x,v)
            push!(x[v], k)
        else
            x[v]=Set(k)
        end
    end
    for (k,v) in m2
        if haskey(y,v)
            push!(y[v],k)
        else
            y[v]=Set(k)
        end
    end
    died_comms = setdiff(keys(x),keys(y))
    new_comms = setdiff(keys(y),keys(x))
    con_comms = intersect(keys(x),keys(y))
    for c in died_comms
        for u in x[c]
            println(fout, "-nc\t", u, '\t', c)
        end
        println(fout, "-c\t", c)
    end
    for c in new_comms
        println(fout, "+c\t", c)
        for u in y[c]
            println(fout, "+nc\t", u, '\t', c)
        end
    end
    for c in con_comms
        new_nodes = setdiff(y[c], x[c])
        for u in new_nodes
            println(fout, "+nc\t", u, '\t', c)
        end
        del_nodes = setdiff(x[c],y[c])
        for u in del_nodes
            println(fout, "-nc\t", u, '\t', c)
        end
    end
end

function convertfile(filein, fileout)
    fin = open(filein, "r")
    fout = open(fileout, "w")
    i = 0
    nmap = Dict{Int,Int}()
    for l in eachline(fin)
        items = split(chop(l))
        if items[1]=="+c"
            i+=1
            nmap[parse(Int,items[2])]=i
            println(fout, items[1], '\t', i)
        end
        if items[1]=="+nc" || items[1]=="-nc"
            println(fout, items[1],'\t', items[2],'\t', nmap[parse(Int,items[3])])
        end
        if items[1]=="-c"
            println(fout, items[1],'\t', nmap[parse(Int,items[2])])
        end
        if items[1]=="+" || items[1]=="-" || items[1][1]=='#'
            print(fout, l)
        end
    end
    close(fout)
    close(fin)
end


text:
#=  https://www.damtp.cam.ac.uk/research/dd/files/teaching/MB17Gog.pdf#subsubsection.1.3.4  =#
using Plots

δt = 0.001
T = 1
times = 0:δt:T
ntimes = length(times)

N = 10000

S = zeros(size(times))
S[1] = N - 1
I = zeros(size(times))
I[1] = 1
R = zeros(size(times))

# Rate of infection
β = 0.005
# Rate of recovery
ν = 4

# Randomness for σ > 0
σ = 1

for i ∈ 2:ntimes
    δS = - β * I[i-1] * S[i-1] * δt
	δR = ν * I[i-1] * δt
	δI = - δS - δR
    S[i] = max(S[i-1] + δS * exp(σ * randn()), 0)
    I[i] = max(I[i-1] + δI * exp(σ * randn()), 0)
	R[i] = max(R[i-1] + δR * exp(σ * randn()), 0)
	scale = (S[i] + I[i] + R[i])/N
	S[i] /= scale; I[i] /= scale; R[i] /= scale
end

lineplot = plot(times, [S I R], label = ["Susceptible" "Infectious" "Recovered"], color = [:gold :red :blue]);
areaplot_ = areaplot(times, [S I R], label = ["Susceptible" "Infectious" "Recovered"], color = [:gold :red :blue]);

display(lineplot)
display(areaplot_)


text:
module Eval

import LsoBase

# include ./eval/*
LsoBase.includedir(dirname(@__FILE__)*"/eval")

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
function add_satellite(obj::SP3File, arg0::JString)
    return jcall(obj, "addSatellite", void, (JString,), arg0)
end

function contains_satellite(obj::SP3File, arg0::JString)
    return jcall(obj, "containsSatellite", jboolean, (JString,), arg0)
end

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

function get_agency(obj::SP3File)
    return jcall(obj, "getAgency", JString, ())
end

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

function get_coordinate_system(obj::SP3File)
    return jcall(obj, "getCoordinateSystem", JString, ())
end

function get_data_used(obj::SP3File)
    return jcall(obj, "getDataUsed", JString, ())
end

function get_day_fraction(obj::SP3File)
    return jcall(obj, "getDayFraction", jdouble, ())
end

function get_epoch(obj::SP3File)
    return jcall(obj, "getEpoch", AbsoluteDate, ())
end

function get_epoch_interval(obj::SP3File)
    return jcall(obj, "getEpochInterval", jdouble, ())
end

function get_gps_week(obj::SP3File)
    return jcall(obj, "getGpsWeek", jint, ())
end

function get_julian_day(obj::SP3File)
    return jcall(obj, "getJulianDay", jint, ())
end

function get_number_of_epochs(obj::SP3File)
    return jcall(obj, "getNumberOfEpochs", jint, ())
end

function get_orbit_type(obj::SP3File)
    return jcall(obj, "getOrbitType", SP3File_SP3OrbitType, ())
end

function get_orbit_type_key(obj::SP3File)
    return jcall(obj, "getOrbitTypeKey", JString, ())
end

function get_satellite_count(obj::SP3File)
    return jcall(obj, "getSatelliteCount", jint, ())
end

function get_satellites(obj::SP3File)
    return jcall(obj, "getSatellites", Map, ())
end

function get_seconds_of_week(obj::SP3File)
    return jcall(obj, "getSecondsOfWeek", jdouble, ())
end

function get_time_system(obj::SP3File)
    return jcall(obj, "getTimeSystem", SP3File_TimeSystem, ())
end

function get_type(obj::SP3File)
    return jcall(obj, "getType", SP3File_SP3FileType, ())
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 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