text:
module TestShow

using DataFrames, Random, Test

function capture_stdout(f::Function)
    oldstdout = stdout
    rd, wr = redirect_stdout()
    f()
    redirect_stdout(oldstdout)
    size = displaysize(rd)
    close(wr)
    str = read(rd, String)
    close(rd)
    str, size
end

@testset "Basic show test with allrows and allcols" begin
    df = DataFrame(A = Int64[1:4;], B = ["x\"", "∀ε>0: x+ε>x", "z\$", "A\nC"],
                   C = Float32[1.0, 2.0, 3.0, 4.0])

    refstr = """
    4×3 DataFrame
    │ Row │ A     │ B           │ C       │
    │     │ Int64 │ String      │ Float32 │
    ├─────┼───────┼─────────────┼─────────┤
    │ 1   │ 1     │ x"          │ 1.0     │
    │ 2   │ 2     │ ∀ε>0: x+ε>x │ 2.0     │
    │ 3   │ 3     │ z\$          │ 3.0     │
    │ 4   │ 4     │ A\\nC        │ 4.0     │"""

    for allrows in [true, false], allcols in [true, false]
        io = IOBuffer()
        show(io, df, allcols=allcols, allrows=allrows)
        str = String(take!(io))
        @test str == refstr
        io = IOBuffer()
        show(io, MIME("text/plain"), df, allcols=allcols, allrows=allrows)
        str = String(take!(io))
        @test str == refstr
    end

    df = DataFrame(A = Vector{String}(undef, 3))
    @test sprint(show, df) == """
    3×1 DataFrame
    │ Row │ A      │
    │     │ String │
    ├─────┼────────┤
    │ 1   │ #undef │
    │ 2   │ #undef │
    │ 3   │ #undef │"""
end

@testset "displaysize test" begin
    Random.seed!(1)
    df_big = DataFrame(rand(Int64(10000000):Int64(20000000), 25,5))

    io = IOContext(IOBuffer(), :displaysize=>(11,40), :limit=>true)
    show(io, df_big)
    str = String(take!(io.io))
    @test str == """
    25×5 DataFrame. Omitted printing of 2 columns
    │ Row │ x1       │ x2       │ x3       │
    │     │ Int64    │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 11569985 │ 12178109 │ 17315979 │
    ⋮
    │ 24  │ 14660095 │ 13529407 │ 19204569 │
    │ 25  │ 16992761 │ 15379139 │ 13043102 │"""

    io = IOContext(IOBuffer(), :displaysize=>(11,40), :limit=>true)
    show(io, df_big, allcols=true)
    str = String(take!(io.io))
    @test str == """
    25×5 DataFrame
    │ Row │ x1       │ x2       │ x3       │
    │     │ Int64    │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 11569985 │ 12178109 │ 17315979 │
    ⋮
    │ 24  │ 14660095 │ 13529407 │ 19204569 │
    │ 25  │ 16992761 │ 15379139 │ 13043102 │

    │ Row │ x4       │ x5       │
    │     │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┤
    │ 1   │ 10701540 │ 17870314 │
    ⋮
    │ 24  │ 19137290 │ 16313933 │
    │ 25  │ 12314843 │ 17754964 │"""

    io = IOContext(IOBuffer(), :displaysize=>(11,40), :limit=>true)
    show(io, df_big, allrows=true, allcols=true)
    str = String(take!(io.io))
    @test str == """
    25×5 DataFrame
    │ Row │ x1       │ x2       │ x3       │
    │     │ Int64    │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 11569985 │ 12178109 │ 17315979 │
    │ 2   │ 19192686 │ 15027856 │ 16050701 │
    │ 3   │ 11140411 │ 13293718 │ 12677906 │
    │ 4   │ 11985459 │ 15660106 │ 13057018 │
    │ 5   │ 19201756 │ 11621062 │ 12785196 │
    │ 6   │ 18146532 │ 15031398 │ 17264273 │
    │ 7   │ 13797788 │ 12442579 │ 14445511 │
    │ 8   │ 10315853 │ 15868969 │ 17213493 │
    │ 9   │ 12916299 │ 18460443 │ 12300445 │
    │ 10  │ 13157897 │ 10639979 │ 10106305 │
    │ 11  │ 16871306 │ 17775325 │ 12544912 │
    │ 12  │ 15251225 │ 10795979 │ 17976274 │
    │ 13  │ 13028558 │ 12220954 │ 17347524 │
    │ 14  │ 13057734 │ 17805919 │ 18183984 │
    │ 15  │ 16120716 │ 12095491 │ 12250720 │
    │ 16  │ 17157319 │ 12697043 │ 11145594 │
    │ 17  │ 19264603 │ 16981136 │ 12984027 │
    │ 18  │ 11292293 │ 19779844 │ 18722894 │
    │ 19  │ 12360733 │ 13117985 │ 11836582 │
    │ 20  │ 13975865 │ 14529221 │ 18458601 │
    │ 21  │ 11416596 │ 17150526 │ 14731764 │
    │ 22  │ 19587351 │ 15978326 │ 17757430 │
    │ 23  │ 10802633 │ 17029758 │ 12105159 │
    │ 24  │ 14660095 │ 13529407 │ 19204569 │
    │ 25  │ 16992761 │ 15379139 │ 13043102 │

    │ Row │ x4       │ x5       │
    │     │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┤
    │ 1   │ 10701540 │ 17870314 │
    │ 2   │ 15507419 │ 16954480 │
    │ 3   │ 16317941 │ 10996749 │
    │ 4   │ 16740306 │ 18240586 │
    │ 5   │ 14628017 │ 14818074 │
    │ 6   │ 14239854 │ 12226254 │
    │ 7   │ 16972931 │ 17422692 │
    │ 8   │ 16291319 │ 18748371 │
    │ 9   │ 11568636 │ 19668632 │
    │ 10  │ 12073431 │ 11550526 │
    │ 11  │ 14037727 │ 19458202 │
    │ 12  │ 12911112 │ 17083732 │
    │ 13  │ 11111638 │ 18826082 │
    │ 14  │ 14534511 │ 18412780 │
    │ 15  │ 19921217 │ 12593752 │
    │ 16  │ 12449162 │ 19192439 │
    │ 17  │ 11091987 │ 17326521 │
    │ 18  │ 16884903 │ 19680264 │
    │ 19  │ 18348781 │ 16089361 │
    │ 20  │ 10189747 │ 13978329 │
    │ 21  │ 19285347 │ 11039190 │
    │ 22  │ 14664380 │ 17880065 │
    │ 23  │ 18398749 │ 12751655 │
    │ 24  │ 19137290 │ 16313933 │
    │ 25  │ 12314843 │ 17754964 │"""

    io = IOContext(IOBuffer(), :displaysize=>(11,40), :limit=>true)
    show(io, df_big, allrows=true, allcols=false)
    str = String(take!(io.io))
    @test str == """
    25×5 DataFrame. Omitted printing of 2 columns
    │ Row │ x1       │ x2       │ x3       │
    │     │ Int64    │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 11569985 │ 12178109 │ 17315979 │
    │ 2   │ 19192686 │ 15027856 │ 16050701 │
    │ 3   │ 11140411 │ 13293718 │ 12677906 │
    │ 4   │ 11985459 │ 15660106 │ 13057018 │
    │ 5   │ 19201756 │ 11621062 │ 12785196 │
    │ 6   │ 18146532 │ 15031398 │ 17264273 │
    │ 7   │ 13797788 │ 12442579 │ 14445511 │
    │ 8   │ 10315853 │ 15868969 │ 17213493 │
    │ 9   │ 12916299 │ 18460443 │ 12300445 │
    │ 10  │ 13157897 │ 10639979 │ 10106305 │
    │ 11  │ 16871306 │ 17775325 │ 12544912 │
    │ 12  │ 15251225 │ 10795979 │ 17976274 │
    │ 13  │ 13028558 │ 12220954 │ 17347524 │
    │ 14  │ 13057734 │ 17805919 │ 18183984 │
    │ 15  │ 16120716 │ 12095491 │ 12250720 │
    │ 16  │ 17157319 │ 12697043 │ 11145594 │
    │ 17  │ 19264603 │ 16981136 │ 12984027 │
    │ 18  │ 11292293 │ 19779844 │ 18722894 │
    │ 19  │ 12360733 │ 13117985 │ 11836582 │
    │ 20  │ 13975865 │ 14529221 │ 18458601 │
    │ 21  │ 11416596 │ 17150526 │ 14731764 │
    │ 22  │ 19587351 │ 15978326 │ 17757430 │
    │ 23  │ 10802633 │ 17029758 │ 12105159 │
    │ 24  │ 14660095 │ 13529407 │ 19204569 │
    │ 25  │ 16992761 │ 15379139 │ 13043102 │"""
end

@testset "IOContext parameters test" begin
    df = DataFrame(A = Int64[1:4;], B = ["x\"", "∀ε>0: x+ε>x", "z\$", "A\nC"],
                   C = Float32[1.0, 2.0, 3.0, 4.0])
    str1, size = capture_stdout() do
        show(df)
    end
    io = IOContext(IOBuffer(), :limit=>true, :displaysize=>size)
    show(io, df)
    str2 = String(take!(io.io))
    @test str1 == str2

    Random.seed!(1)
    df_big = DataFrame(rand(25,5))
    str1, size = capture_stdout() do
        show(df_big)
    end
    io = IOContext(IOBuffer(), :limit=>true, :displaysize=>size)
    show(io, df_big)
    str2 = String(take!(io.io))
    @test str1 == str2
end

@testset "SubDataFrame show test" begin
    df = DataFrame(A = Int64[1:4;], B = ["x\"", "∀ε>0: x+ε>x", "z\$", "A\nC"],
                   C = Float32[1.0, 2.0, 3.0, 4.0])
    subdf = view(df, [2, 3], :)
    io = IOBuffer()
    show(io, subdf, allrows=true, allcols=false)
    str = String(take!(io))
    @test str == """
    2×3 SubDataFrame
    │ Row │ A     │ B           │ C       │
    │     │ Int64 │ String      │ Float32 │
    ├─────┼───────┼─────────────┼─────────┤
    │ 1   │ 2     │ ∀ε>0: x+ε>x │ 2.0     │
    │ 2   │ 3     │ z\$          │ 3.0     │"""
    show(io, subdf, allrows=true)
    show(io, subdf, allcols=true)
    show(io, subdf, allcols=true, allrows=true)
end

@testset "Test showing StackedVector and RepeatedVector" begin
    A = DataFrames.StackedVector(Any[[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    @test sprint(show, A) == "[1, 2, 3, 4, 5, 6, 7, 8, 9]"
    A = DataFrames.RepeatedVector([1, 2, 3], 5, 1)
    @test sprint(show, A) == "[1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]"
    A = DataFrames.RepeatedVector([1, 2, 3], 1, 5)
    @test sprint(show, A) == "[1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]"
end

@testset "Test colors and non-standard values: missing and nothing" begin
    df = DataFrame(Fish = ["Suzy", "Amir"], Mass = [1.5, missing])
    @test sprint(show, df, context=:color=>true) == """
        2×2 DataFrame
        │ Row │ Fish   │ Mass     │
        │     │ \e[90mString\e[39m │ \e[90mFloat64⍰\e[39m │
        ├─────┼────────┼──────────┤
        │ 1   │ Suzy   │ 1.5      │
        │ 2   │ Amir   │ \e[90mmissing\e[39m  │"""

    df = DataFrame(A = [:Symbol, missing, :missing],
                   B = [missing, "String", "missing"],
                   C = [:missing, "missing", missing])
    @test sprint(show, df, context=:color=>true) == """
        3×3 DataFrame
        │ Row │ A       │ B       │ C       │
        │     │ \e[90mSymbol⍰\e[39m │ \e[90mString⍰\e[39m │ \e[90mAny\e[39m     │
        ├─────┼─────────┼─────────┼─────────┤
        │ 1   │ Symbol  │ \e[90mmissing\e[39m │ missing │
        │ 2   │ \e[90mmissing\e[39m │ String  │ missing │
        │ 3   │ missing │ missing │ \e[90mmissing\e[39m │"""

    df_nothing = DataFrame(A = [1.0, 2.0, 3.0], B = ["g", "g", nothing])
    @test sprint(show, df_nothing) == """
    3×2 DataFrame
    │ Row │ A       │ B      │
    │     │ Float64 │ Union… │
    ├─────┼─────────┼────────┤
    │ 1   │ 1.0     │ g      │
    │ 2   │ 2.0     │ g      │
    │ 3   │ 3.0     │        │"""
end

@testset "Test correct width computation" begin
    df = DataFrame([["a"]], [:x])
    @test sprint(show, df) == """
    1×1 DataFrame
    │ Row │ x      │
    │     │ String │
    ├─────┼────────┤
    │ 1   │ a      │"""
end

@testset "Test showing special types: strings with escapes, categorical and BigFloat" begin
    df = DataFrame(a = ["1\n1", "2\t2", "3\r3", "4\$4", "5\"5", "6\\6"])
    @test sprint(show, df) == """
    6×1 DataFrame
    │ Row │ a      │
    │     │ String │
    ├─────┼────────┤
    │ 1   │ 1\\n1   │
    │ 2   │ 2\\t2   │
    │ 3   │ 3\\r3   │
    │ 4   │ 4\$4    │
    │ 5   │ 5"5    │
    │ 6   │ 6\\\\6   │"""

    df = DataFrame(a = categorical([1,2,3]), b = categorical(["a", "b", missing]))
    @test sprint(show, df) == """
    3×2 DataFrame
    │ Row │ a            │ b             │
    │     │ Categorical… │ Categorical…⍰ │
    ├─────┼──────────────┼───────────────┤
    │ 1   │ 1            │ a             │
    │ 2   │ 2            │ b             │
    │ 3   │ 3            │ missing       │"""

    df = DataFrame(a = [big(1.0), missing])
    @test sprint(show, df) == """
    2×1 DataFrame
    │ Row │ a         │
    │     │ BigFloat⍰ │
    ├─────┼───────────┤
    │ 1   │ 1.0       │
    │ 2   │ missing   │"""
end

@testset "Test using :compact parameter of IOContext" begin
    df = DataFrame(x = [float(pi)])
    @test sprint(show, df) == """
        1×1 DataFrame
        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint(show, df, context=:compact=>false) == """
        1×1 DataFrame
        │ Row │ x                 │
        │     │ Float64           │
        ├─────┼───────────────────┤
        │ 1   │ 3.141592653589793 │"""
end

@testset "Test of DataFrameRows and DataFrameColumns" begin
    df = DataFrame(x = [float(pi)])
    @test sprint(show, eachrow(df)) == """
        1×1 DataFrameRows
        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint((io, x) -> show(io, x, summary=false), eachrow(df)) == """

        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint(show, eachcol(df)) == """
        1×1 DataFrameColumns (with names=false)
        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint((io, x) -> show(io, x, summary=false), eachcol(df)) == """

        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint(show, eachcol(df, true)) == """
        1×1 DataFrameColumns (with names=true)
        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint((io, x) -> show(io, x, summary=false), eachcol(df, true)) == """

        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""
end

@testset "Test empty data frame and DataFrameRow" begin
    df = DataFrame(x = [float(pi)])
    @test sprint(show, df[:, 2:1]) == "0×0 DataFrame\n"
    @test sprint(show, @view df[:, 2:1]) == "0×0 SubDataFrame\n"
    @test sprint(show, df[1, 2:1]) == "DataFrameRow"
end

@testset "consistency" begin
    df = DataFrame(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(df.c, 5)
    @test_throws AssertionError sprint(show, df)

    df = DataFrame(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(DataFrames._columns(df), df[:, :a])
    @test_throws AssertionError sprint(show, df)
end

end # module


text:
using Zygote, Flux, NNlib, Random
using Flux: onehot, chunk, batchseq, crossentropy
using Base.Iterators: partition
using StatsBase: wsample

# Optimizer and data loading utilities
include("utils.jl")
include("optimise.jl")

function train_model!(model, Xs, Ys, num_epochs = 100)
    opt = ADAM(0.001)
    avg_loss = 0.0

    # Run through our entire dataset a few times
    for epoch_idx in 1:num_epochs
        permutation = shuffle(1:length(Xs))
        for batch_idx in 1:length(permutation)
            x_batch, y_batch = Xs[permutation[batch_idx]], Ys[permutation[batch_idx]]
            # Calculate gradients upon the model for this batch of data,
            # summing crossentropy loss across time
            l, back = Zygote.forward(model) do model
				return sum(crossentropy.(model.(x_batch), y_batch))
			end
            grads = back(1f0)[1]

            # Update the model, then reset its internal LSTM states
            model = update!(opt, model, grads)
            Flux.reset!(model)

            avg_loss = avg_loss*0.98 + 0.02*l
            @info(epoch_idx, batch_idx, l, avg_loss)
        end
        @info("Done with epoch $(epoch_idx)!")
        println(sample(model, alphabet, 300))
    end

    return model
end

# Download dataset if it does not already exist
if !isfile("shakespeare_input.txt")
    download("https://cs.stanford.edu/people/karpathy/char-rnn/shakespeare_input.txt", "shakespeare_input.txt")
end

# Load data and alphabet
alphabet, Xs, Ys = load_data("shakespeare_input.txt")

# Define simple LSTM-based model to map from alphabet back on to alphabet,
# predicting the next letter in a corpus of Shakespeare text.
model = mapleaves(Flux.data, Chain(
    LSTM(length(alphabet), 512),
    LSTM(512, 512),
    Dense(512, length(alphabet)),
	softmax,
))

model = train_model!(model, Xs, Ys)

println("Behold, I speak:")
println(sample(model, alphabet, 1000))

ext: jl

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
using DataFrames
using CSV
using Random
using Plots
using PDMats
using Printf
include("transforms.jl")
include("statistics.jl")
include("validation.jl")
include("plotting.jl")

function gridsearch(f, a, b, numpts=45)
    h = (b-a)/(numpts-1)
    l, m = 0, 0
    for i = a:h:b
        val = f(i);
        if val>m
            m = val; l = i;
        end
    end
    return l, m
end

#macroeconomics time series data
df = DataFrame(CSV.File("data\\TB3MS.csv"))

z = convert(Matrix, df[:,2:2])
n = length(z)
s = collect(1:n) #data points equally spaced in time
portion_train = 0.3

rng = MersenneTwister(1234)
ind = shuffle(rng, s)
id_train = view(ind, 1:Int(ceil(n*portion_train)))
id_test = view(ind, (Int(ceil(n*portion_train))+1):n)

train_x = s[id_train]; test_x = s[id_test]
train_y = z[id_train]; test_y = z[id_test]

X = ones(length(train_x), 1); X0 = [1]
pλ = x -> 1; pθ = x -> 1
range_theta = [0.1 5]
range_lambda = [-5 5]
results = Array{Float64}(undef, length(test_x))#preallocate space for medians and interval widths
intwidths = Array{Float64}(undef, length(test_x))

if true #test and time pdfn and cdfn functions
pdf = model(X, X0, train_x, [test_x[6]], boxCox, boxCoxPrime, pθ, pλ, train_y, range_theta, range_lambda)
cdf = z0 ->  int1D(pdf, 0, z0, "0")
@time begin
un = cdf(30)
end
cdfn = x -> cdf(x)/un #normalized cdf
pdfn = x -> pdf(x)/un #normalized pdf
function testcdfn(x)
    @time begin
        r = cdfn(x)
    end
    return r
end

function testbisection(f, a, b)
    @time begin
        val = bisection(f, a, b)
    end
    return val
end

function testconfidence(f, med)
    @time begin
        val = confidence(f, med)
    end
    return val
end
plt2(pdfn, 0, 20)
end

if false #simply plot at multiple locations
    pltX = zeros(20, length(test_x))
    pltY = zeros(20, length(test_x))
    for i=1:length(test_x)#single-point predictions
        println(i)
        pdf = model(X, X0, train_x, [test_x[i]], boxCox, boxCoxPrime, pθ, pλ, train_y, range_theta, range_lambda)
        cdf = z0 ->  int1D(pdf, 0, z0, "1")
        un = cdf(30)
        cdfn = x -> cdf(x)/un #normalized cdf
        pdfn = x -> pdf(x)/un #normalized pdf
        x, y = plt(pdfn, 5, 18)
        pltX[:, i] = x; pltY[:, i] = y
    end
    display(plot(pltX, pltY, layout = length(test_x), ylim = (0, )))
end

if false #determine medians and confidence intervals
for i=1:length(test_x)#single-point predictions
    println(i)
    @time begin
    pdf = model(X, X0, train_x, [test_x[i]], boxCox, boxCoxPrime, pθ, pλ, train_y, range_theta, range_lambda)
    cdf = z0 ->  int1D(pdf, 0, z0, "0")
    @time begin
    un = cdf(25)
    end
    cdfn = x -> cdf(x)/un #normalized cdf
    pdfn = x -> pdf(x)/un #normalized pdf
    @time begin
    med = bisection(x -> cdfn(x)-0.5, 1e-3, 25, 1e-2, 9)
    #med, val = gridsearch(pdfn, 0, 30, 20) #MAP estimation
    @printf("ith pt med: %f \n", med)
    end
    @time begin
    intwidths[i] = confidence(pdfn, med) #compute width of credible interval corresponding to 95% confidence level
    @printf("ith pt intwidth: %f \n", intwidths[i])
    end
    results[i] = med
    end
end
display(plot(collect(1:n), z));
display(plot!(train_x, train_y, seriestype = scatter));
display(plot!(test_x, results, seriestype = :scatter));
#display(plot!(test_x, test_y, seriestype = :scatter));
display(plot!(repeat(test_x, 2), vcat(results-intwidths, results+intwidths), seriestype = :scatter, color=:blue))
#display(plot!(test_x, test_y))
savefig("treasury6.png");
end

if true #cross validation on training set
    _, Xs, Ys = cross_validate(X, train_x, boxCox, boxCoxPrime, pθ, pλ, train_y, range_theta, range_lambda, 500, 2, 24)
    display(plot(Xs, Ys, layout = length(train_y), legend=nothing, xtickfont = font(10, "Courier"), lw=0.5))
    display(plot!(train_y', layout = length(z), legend=nothing, seriestype = :vline, xtickfont = font(10, "Courier"), lw=0.5))
    savefig("results//treasury//treasury_delete_one_cross_validate2.jpg")
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Pkg
Pkg.add([
    "DFTK", "DoubleFloats", "FiniteDiff", "ForwardDiff",
    "GenericLinearAlgebra", "Infiltrator", "IntervalArithmetic",
    "KrylovKit", "LineSearches", "NLsolve", "Plots", "PyCall",
    "Unitful", "UnitfulAtomic"
])

using PyCall

if !isempty(PyCall.python)
    run(`$(PyCall.python) -m pip install pymatgen`)
end

using DFTK
DFTK.setup_threading() # print threading information


text:
using Base64

"""
    RawPoint

A point of data representing a single position within a time series. Each
point contains a read-only time and value attribute.

**Constructors**

```julia
RawPoint(time::Int64, value::Float64)
RawPoint(data::Dict{String,Any})
RawPoint(time, value)
```

**Arguments**

* time : the time (`Int64`) of a single value in the time series (in nanoseconds since the Unix epoch)
* value : the value (`Int64`) of a time series at a single point in time

**Notes**

The `RawPoint` constructors are used internally by calls such as `values`.

"""
struct RawPoint
    time::Int64
    value::Float64
end

"""
    StatPoint

An aggregated data point representing a summary or rollup of one or more
points of data within a single time series.

This aggregation point provides for the min, mean, max, count, and standard
deviation of all data values it spans. It is returned by windowing queries
such as windows or aligned_windows.

**Constructors**

```julia
StatPoint(time::Int64, min::Float64, mean::Float64, max::Float64,
          count::Int64, stddev::Float64)
StatPoint(data::Dict{String,Any})
StatPoint(time, min, mean, max, count, stddev)
```

**Arguments**

* time : a `Int64` for the time span represented by the aggregated values (in nanoseconds since the Unix epoch)
* min : a `Float64` representing the minimum value of points in this window
* mean : a `Float64` representing the average value of points in this window
* max : a `Float64` representing the maximum value of points in this window
* count : a `Int64` for the number of real values in this window
* stddev : a `Float64` representing the standard deviation of point values in this window

**Notes**

The `StatPoint` constructors are used internally by aggregation calls such as `windows` and `aligned_windows`.
"""
struct StatPoint
    time::Int64
    min::Float64
    mean::Float64
    max::Float64
    count::Int64
    stddev::Float64
end

"""
    Stream

An object that represents a specific time series stream in the BTrDB database.

**Constructors**

```julia
Stream(uuid::String, name::String, collection::String, tags::Dict{String,String},
       annotations::Dict{String,Any}, version::Int64, propertyVersion::Int64)
Stream(data::Dict{String,Any})
Stream(uuid, name, collection, tags, annotations, version, propertyVersion)
```

**Arguments**

* uuid : a `String` of the stream's unique identifier
* name : a `String` of the stream's friendly name
* collection : a `String` of the stream's collection (path) in the stream hierarchy
* tags : a `Dict{String,String}` of the stream's (internal use) metadata
* annotations : a `Dict{String,Any}` of the stream's public metadata
* version : a `Int64` that acts as a monotonically increasing version of the data
* propertyVersion : a `Int64` that acts as a monotonically increasing version of the metadata

**Notes**

The `Stream` constructors are used internally by any calls that will return Stream objects.
In general, this is unlikely to be used by end users but is still provided in the public API.

"""
mutable struct Stream
    uuid::String
    name::String
    collection::String
    tags::Dict{String,String}
    annotations::Dict{String,Any}
    version::Int64
    propertyVersion::Int64
end

###############################################################################
# Constructors
###############################################################################

function Stream(data::Dict{String,Any})
    tags = decodeTags(data["tags"])
    Stream(
        decodeUUID(data["uuid"]),
        tags["name"],
        data["collection"],
        tags,
        decodeAnnotations(data["annotations"]),
        0,
        0
    )
end

"""
    RawPoint

    A point of data representing a single position within a time series. Each point contains a read-only time and value attribute.
"""
function RawPoint(data::Dict{String,Any})
    return RawPoint(parse(Int64, data["time"]), data["value"])
end

"""
StatPoint

An aggregated data point representing a summary or rollup of one or more points of data within a single time series.

This aggregation point provides for the min, mean, max, count, and standard deviation of all data values it spans. It is returned by windowing queries such as windows or aligned_windows.
"""
function StatPoint(data::Dict{String,Any})
    if typeof(data["stddev"]) == String
        stddev = NaN
    else
        stddev = data["stddev"]
    end
    return StatPoint(
        parse(Int64, data["time"]),
        data["min"],
        data["mean"],
        data["max"],
        parse(Int64, data["count"]),
        stddev
    )
end

###############################################################################
# Module Only Utility Functions
###############################################################################

function dict2array(tags::Dict)
    data = []
    for (k, v) in tags
        push!(data, Dict{String, Any}("key" => k, "val" => Dict{String, String}("value" => v)))
    end
    return data
end

function decodeTags(data)
    tags = Dict{String, String}()
    for item in data
        tags[item["key"]] = item["val"]["value"]
    end
    return tags
end

function decodeAnnotations(data)
    annotations = Dict{String, Any}()
    for item in data
        annotations[item["key"]] = item["val"]["value"]
    end
    return annotations
end


text:
function master_formulation(rd::RoutingData, ::PathFormulation)
    m = _create_model(rd)
    rm = basic_routing_model_unitary(m, rd, PathFormulation()) # Includes convexity constraint.

    @variable(m, mu >= 0)
    @objective(m, Min, mu)

    if rd.model_robust_reformulation_traffic_matrices
        @variable(m, dual_alpha[e in edges(rd), k in demands(rd)] <= 0)
        @variable(m, dual_beta[e in edges(rd), e2 in edges(rd)] >= 0)
    else
        @variable(m, dual[e in edges(rd), e2 in edges(rd)] >= 0)
    end

    # Robust uncertainty set.
    if rd.model_robust_reformulation_traffic_matrices
        dual_constraints = Dict{Edge{Int}, Dict{Edge{Int}, ConstraintRef}}(
            e => Dict{Edge{Int}, ConstraintRef}() for e in edges(rd)
        ) # Edge -> demand -> constraint reference.
    end
    for e in edges(rd)
        # Uncertainty set for one edge.
        for d in demands(rd)
            rhs = AffExpr(0.0)
            for p in rd.demand_to_path_ids[d]
                if e in rd.paths_edges[p]
                    add_to_expression!(rhs, rm.routing[d, p])
                end
            end

            if rd.model_robust_reformulation_traffic_matrices
                dual_constraints[e][d] = @constraint(m, rhs + dual_alpha[e, d] <= 0)
                for p in rd.demand_to_path_ids[d]
                    @constraint(m, dual_alpha[e, d] + sum(dual_beta[e, e2] for e2 in rd.paths_edges[p]) >= 0)
                end
            else
                for p in rd.demand_to_path_ids[d]
                    @constraint(m, rhs <= sum(dual[e, e2] for e2 in rd.paths_edges[p]))
                end
            end
        end

        # Relate the main decision variables to the uncertainty sets.
        let dual_var = (rd.model_robust_reformulation_traffic_matrices ? dual_beta : dual)
            rhs = sum(dual_var[e, e2] * capacity(rd, e2) for e2 in edges(rd))
            rhs /= capacity(rd, e)
            @constraint(m, rhs <= mu)
        end
    end

    if rd.model_robust_reformulation_traffic_matrices
        return RoutingModel(rd, m, UnitaryFlows, rm.routing, mu=mu,
                            dual_alpha=dual_alpha, dual_beta=dual_beta,
                            constraints_uncertainty_set=dual_constraints,
                            constraints_convexity=rm.constraints_convexity)
    else
        return RoutingModel(rd, m, UnitaryFlows, rm.routing, mu=mu, dual=dual,
                            constraints_convexity=rm.constraints_convexity)
    end
end


text:
"""
    function test_2_breeding(base, qtl, nsib)
---
- 100 sires, 200 dams
- sibship size: 10, 20, 40 for both produciton and challenge group
- edit 20,000 fertilized eggs with a rate
  * also chech Hicky's editing method
- compare
  * no editing
  * editing
- order QTL on `2pqa²`.
- equal weights on both traits.
- 2021-March-06
"""
function test_2_breeding(base, qtl, nsib, weight)
    # It is better to setup simulation scenarios with a dictionary
    # Then pass it as a named tuple, so that we can conveniantly use `x.y`
    # using a dictionary also makes it easy to modify parameter structure.
    # the keys in this dictionry must exist to setup a simulation
    par = begin
        Parameters = Dict(:nSire=> 100,
                          :nDam => 200, # male:female = 1:2
                          :nSib => 2nsib,
                          :nC7e => nsib, # number of sib to be challenged
                          :nG8n => 10,   # number of generations
                          :nQTL => [1000, 1000],
                          :h²   => [.5, .5],
                          :p8e  => .5, # percentage of affected in the binary trait
                          :e19e => .8, # edit_successsful_rate
                          :w4t  => weight # weight on the binary trait EBV
                          )
        (; Parameters...)       # named tuple.  contents as above
    end
    prd, _ = breeding_program(base, par, qtl)
    serialize(joinpath(dat_dir, "run/breeding/w/$weight-$nsib.ser"), prd)
end

"""
    function sum_w_breeding()
---
Summarize genetic improvement using different weight on GEBV of binary trait.
The conclusion is `weight` should be some where between 2 and 3.
- 2021-Mar-11
"""
function sum_w_breeding()
    sib    = [10, 20, 40]
    weight = [3, 6, 9, 12]
    ave = Array{Float64, 4}(undef, 3, 4, 2, 10)
    for i in 1:3
        for j in 1:4
            rst = deserialize("dat/run/breeding/w/$(weight[j])-$(sib[i]).ser")
            k = 1
            for d in groupby(rst, :g8n)
                ave[i, j, 1, k] = mean(d.tbv)
                ave[i, j, 2, k] = mean(d.t2c)
                k += 1
            end
        end
    end
    for i in 1:3
        ave[i, 1, 1, :] .-= ave[i, 1, 1, 1]
        plot(ave[i, 1, 1, :], label="prd-$(sib[i])s-3w", dpi=300, legend=:topleft)
        for j in 2:4
            ave[i, j, 1, :] .-= ave[i, j, 1, 1]
            plot!(ave[i, j, 1, :], label="prd-$(sib[i])s-$(weight[j])w")
        end
        for j in 1:4
            ave[i, j, 2, :] .-= ave[i, j, 2, 1]
            plot!(ave[i, j, 2, :], label="tst-$(sib[i])s-$(weight[j])w", linestyle=:dash)
        end
        savefig(joinpath(dat_dir, "run/breeding/w/sib-$(sib[i]).pdf"))
    end
end

"""
    function sum_2_breeding()
---
Function `test_2_breeding` of v0.5.2 resulted in `breeding/{1,2,3}'.
This function is to summarize the results in 9 files.
Below was done Mar. 08, 2021
"""
function sum_2_breeding()
    for i in 1:3
        m = Float64[]
        for sib in [10, 20, 40]
            dat = deserialize(joinpath(dat_dir, "run/breeding/$i/$sib.ser"))
            for d in groupby(dat, :g8n)
                push!(m, mean(d.tbv))
                push!(m, mean(d.t2c))
            end
        end
        m = reshape(m, 6, :)
        t = reshape(m[1:20], 2, :)
        t .-= t[:, 1]           # set the first generation as 0.
        plot(t[1, :], label="Production-10-sibs", dpi=300, legend=:topleft)
        plot!(t[2, :], label="Challenge-10-sibs")
        t = reshape(m[21:40], 2, :)
        t .-= t[:, 1]           # set the first generation as 0.
        plot!(t[1, :], label="Production-20-sibs")
        plot!(t[2, :], label="Challenge-20-sibs")
        t = reshape(m[41:60], 2, :)
        t .-= t[:, 1]           # set the first generation as 0.
        plot!(t[1, :], label="Production-40-sibs")
        plot!(t[2, :], label="Challenge-40-sibs")
        savefig(joinpath(dat_dir, "run/breeding/fig/$i.pdf"))
    end
end

"""
    function test_breeding()
---
This is a independent function.  It tests selection on
TBV,
phenotypes,
and random variable (no selection).
The result is a pdf figure.
So far so good.
"""
function test_breeding()
    Parameters = Dict(:nSire=> 125,
                      :nDam => 250, # male:female = 1:2
                      :nSib => 100,
                      :nC7e => 40, # number of sib to be challenged
                      :nG8n => 40, # number of generations
                      :nQTL => [1000, 1000],
                      :h²   => [.2, .2],
                      :p8e  => .5, # percentage of affected in the binary trait
                      )
    par = (; Parameters...)     # named tuple.  contents as above

    # read back real data
    #base = deserialize(joinpath(dat_dir, "run/ns.ser")) # to save time
    @load "dat/run/base.jld2"

    # general parameters
    qtl = sim_QTL(base, par.nQTL...)
    @info "Create storage and generation one"
    snp₁, snp₂, prd, clg, f1, f2, f3, f4 = # f: from
        create_storage(base, par, qtl)
    # shorthands
    nSire, nDam, nSib, nChallenge, percentage, generation =
        par.nSire, par.nDam, par.nSib, par.nC7e, par.p8e, par.nG8n
    ip, ic, hp, hc = test_sets(nDam, nSib, nChallenge)
    l1, l2, l3, l4 = length.((ip, ic, hp, hc)) # lengths of generation 2+
    mp = zeros(generation)      # mean phenotype on generation
    average(ig) = begin
        tmp = select(prd, :g8n, :p7e)
        mean(tmp.p7e[tmp.g8n .== ig])
    end
    mp[1] = average(1)

    @info "!!!!! Select on TBV"
    a, b, c, d = f1, f2, f3, f4
    for ig in 2:generation
        @info "generation $ig"
        df = select(prd[(prd.g8n .== ig-1), :], :id, :sex, :tbv => :val)
        nuclear = select_nuclear(df, nSire, nDam)

        # simulate current generation
        ped = random_mate(nuclear, nSire, nDam, nSib)
        nxt = gdrop(snp₁, ped, base[:r])
        bv₁, p₁ = phenotype(nxt, qtl[1], par.h²[1])
        bv₂, p₂ = phenotype(nxt, qtl[2], par.h²[2], par.p8e)

        # update storage.  codes below are ugly.
        prd.tbv[f1+1:f1+l1] = bv₁[ip]
        prd.p7e[f1+1:f1+l1] =  p₁[ip]
        prd.sir[f1+1:f1+l1] = ped[ip, 1]
        prd.dam[f1+1:f1+l1] = ped[ip, 2]
        clg.tbv[f2+1:f2+l2] = bv₂[ic]
        clg.p7e[f2+1:f2+l2] =  p₂[ic]
        snp₁[:, f3+1:f3+l3] = nxt[:, hp]
        snp₂[:, f4+1:f4+l4] = nxt[:, hc]
        f1 += l1
        f2 += l2
        f3 += l3
        f4 += l4
        mp[ig] = average(ig)
    end
    plot(mp, label="TBV", legend=:topleft)

    @info "!!!!! Select on phenotypes"
    f1, f2, f3, f4 = a, b, c, d
    for ig in 2:generation
        @info "generation $ig"
        df = select(prd[(prd.g8n .== ig-1), :], :id, :sex, :p7e => :val)
        nuclear = select_nuclear(df, nSire, nDam)

        # simulate current generation
        ped = random_mate(nuclear, nSire, nDam, nSib)
        nxt = gdrop(snp₁, ped, base[:r])
        bv₁, p₁ = phenotype(nxt, qtl[1], par.h²[1])
        bv₂, p₂ = phenotype(nxt, qtl[2], par.h²[2], par.p8e)

        # update storage.  codes below are ugly.
        prd.tbv[f1+1:f1+l1] = bv₁[ip]
        prd.p7e[f1+1:f1+l1] =  p₁[ip]
        prd.sir[f1+1:f1+l1] = ped[ip, 1]
        prd.dam[f1+1:f1+l1] = ped[ip, 2]
        clg.tbv[f2+1:f2+l2] = bv₂[ic]
        clg.p7e[f2+1:f2+l2] =  p₂[ic]
        snp₁[:, f3+1:f3+l3] = nxt[:, hp]
        snp₂[:, f4+1:f4+l4] = nxt[:, hc]
        f1 += l1
        f2 += l2
        f3 += l3
        f4 += l4
        mp[ig] = average(ig)
    end
    plot!(mp, label="Phenotype")

    @info "!!!!! Random mating"
    f1, f2, f3, f4 = a, b, c, d
    for ig in 2:generation
        @info "generation $ig"
        # simulate current generation
        df = select(prd[(prd.g8n .== ig-1), :], :id, :sex)
        df.val = rand(length(df.id))
        nuclear = select_nuclear(df, nSire, nDam)
        ped = random_mate(nuclear, nSire, nDam, nSib)
        nxt = gdrop(snp₁, ped, base[:r])
        bv₁, p₁ = phenotype(nxt, qtl[1], par.h²[1])
        bv₂, p₂ = phenotype(nxt, qtl[2], par.h²[2], par.p8e)

        # update storage.  codes below are ugly.
        prd.tbv[f1+1:f1+l1] = bv₁[ip]
        prd.p7e[f1+1:f1+l1] =  p₁[ip]
        prd.sir[f1+1:f1+l1] = ped[ip, 1]
        prd.dam[f1+1:f1+l1] = ped[ip, 2]
        clg.tbv[f2+1:f2+l2] = bv₂[ic]
        clg.p7e[f2+1:f2+l2] =  p₂[ic]
        snp₁[:, f3+1:f3+l3] = nxt[:, hp]
        snp₂[:, f4+1:f4+l4] = nxt[:, hc]
        f1 += l1
        f2 += l2
        f3 += l3
        f4 += l4
        mp[ig] = average(ig)
    end
    plot!(mp, label="Random mating")
    savefig("doc/fig/test-breeding-v3.pdf")
end

"""
    function test_breeding_sim_base()
---
This is a copy of above function, except
1. the base was replace with a simulated one.
2. figures were drawn somewhere else
"""
function test_breeding_sim_base()
    Parameters = Dict(:nSire=> 125,
                      :nDam => 250, # male:female = 1:2
                      :nSib => 100,
                      :nC7e => 40, # number of sib to be challenged
                      :nG8n => 20, # number of generations
                      :nQTL => [1000, 1000],
                      :h²   => [.2, .2],
                      :p8e  => .5, # percentage of affected in the binary trait
                      )
    par = (; Parameters...)     # named tuple.  contents as above

    base = sim_base(1200, 29, 1e8, 1725)

    # general parameters
    qtl = sim_QTL(base, par.nQTL...)
    @info "Create storage and generation one"
    snp₁, snp₂, prd, clg, f1, f2, f3, f4 = # f: from
        create_storage(base, par, qtl)
    # shorthands
    nSire, nDam, nSib, nChallenge, percentage, generation =
        par.nSire, par.nDam, par.nSib, par.nC7e, par.p8e, par.nG8n
    ip, ic, hp, hc = test_sets(nDam, nSib, nChallenge)
    l1, l2, l3, l4 = length.((ip, ic, hp, hc)) # lengths of generation 2+
    mp = zeros(generation)      # mean phenotype on generation
    average(ig) = begin
        tmp = select(prd, :g8n, :p7e)
        mean(tmp.p7e[tmp.g8n .== ig])
    end
    mp[1] = average(1)

    @info "!!!!! Select on TBV"
    a, b, c, d = f1, f2, f3, f4
    for ig in 2:generation
        @info "generation $ig"
        df = select(prd[(prd.g8n .== ig-1), :], :id, :sex, :tbv => :val)
        nuclear = select_nuclear(df, nSire, nDam)

        # simulate current generation
        ped = random_mate(nuclear, nSire, nDam, nSib)
        nxt = gdrop(snp₁, ped, base[:r])
        bv₁, p₁ = phenotype(nxt, qtl[1], par.h²[1])
        bv₂, p₂ = phenotype(nxt, qtl[2], par.h²[2], par.p8e)

        # update storage.  codes below are ugly.
        prd.tbv[f1+1:f1+l1] = bv₁[ip]
        prd.p7e[f1+1:f1+l1] =  p₁[ip]
        prd.sir[f1+1:f1+l1] = ped[ip, 1]
        prd.dam[f1+1:f1+l1] = ped[ip, 2]
        clg.tbv[f2+1:f2+l2] = bv₂[ic]
        clg.p7e[f2+1:f2+l2] =  p₂[ic]
        snp₁[:, f3+1:f3+l3] = nxt[:, hp]
        snp₂[:, f4+1:f4+l4] = nxt[:, hc]
        f1 += l1
        f2 += l2
        f3 += l3
        f4 += l4
        mp[ig] = average(ig)
    end
    plot(mp, label="TBV", legend=:topleft)

    @info "!!!!! Select on phenotypes"
    f1, f2, f3, f4 = a, b, c, d
    for ig in 2:generation
        @info "generation $ig"
        df = select(prd[(prd.g8n .== ig-1), :], :id, :sex, :p7e => :val)
        nuclear = select_nuclear(df, nSire, nDam)

        # simulate current generation
        ped = random_mate(nuclear, nSire, nDam, nSib)
        nxt = gdrop(snp₁, ped, base[:r])
        bv₁, p₁ = phenotype(nxt, qtl[1], par.h²[1])
        bv₂, p₂ = phenotype(nxt, qtl[2], par.h²[2], par.p8e)

        # update storage.  codes below are ugly.
        prd.tbv[f1+1:f1+l1] = bv₁[ip]
        prd.p7e[f1+1:f1+l1] =  p₁[ip]
        prd.sir[f1+1:f1+l1] = ped[ip, 1]
        prd.dam[f1+1:f1+l1] = ped[ip, 2]
        clg.tbv[f2+1:f2+l2] = bv₂[ic]
        clg.p7e[f2+1:f2+l2] =  p₂[ic]
        snp₁[:, f3+1:f3+l3] = nxt[:, hp]
        snp₂[:, f4+1:f4+l4] = nxt[:, hc]
        f1 += l1
        f2 += l2
        f3 += l3
        f4 += l4
        mp[ig] = average(ig)
    end
    plot!(mp, label="Phenotype")

    @info "!!!!! Random mating"
    f1, f2, f3, f4 = a, b, c, d
    for ig in 2:generation
        @info "generation $ig"
        # simulate current generation
        df = select(prd[(prd.g8n .== ig-1), :], :id, :sex)
        df.val = rand(length(df.id))
        nuclear = select_nuclear(df, nSire, nDam)
        ped = random_mate(nuclear, nSire, nDam, nSib)
        nxt = gdrop(snp₁, ped, base[:r])
        bv₁, p₁ = phenotype(nxt, qtl[1], par.h²[1])
        bv₂, p₂ = phenotype(nxt, qtl[2], par.h²[2], par.p8e)

        # update storage.  codes below are ugly.
        prd.tbv[f1+1:f1+l1] = bv₁[ip]
        prd.p7e[f1+1:f1+l1] =  p₁[ip]
        prd.sir[f1+1:f1+l1] = ped[ip, 1]
        prd.dam[f1+1:f1+l1] = ped[ip, 2]
        clg.tbv[f2+1:f2+l2] = bv₂[ic]
        clg.p7e[f2+1:f2+l2] =  p₂[ic]
        snp₁[:, f3+1:f3+l3] = nxt[:, hp]
        snp₂[:, f4+1:f4+l4] = nxt[:, hc]
        f1 += l1
        f2 += l2
        f3 += l3
        f4 += l4
        mp[ig] = average(ig)
    end
    plot!(mp, label="Random mating")
    savefig("doc/fig/test-breeding-v4.pdf")
end


text:
@inline triangle_area(n1::AbstractVec,n2::AbstractVec,n3::AbstractVec) =  0.5*norm((n2-n1) × (n3-n1))

@inline triangle_center(n1::AbstractVec,n2::AbstractVec,n3::AbstractVec) = (n1+n2+n3)/3

@inline triangle_perimeter(n1::AbstractVec,n2::AbstractVec,n3::AbstractVec) = norm(n2-n1)+norm(n3-n2)+norm(n1-n3)

@inline normal_to_2Dline(n1::Vec2D,n2::Vec2D) = Vec2D(n2.y - n1.y, n1.x - n2.x)
@inline function normal_to_2Dline(I::NTuple{2,Integer},nodes::AbstractArray{<:Vec2D,1}) ; @inbounds r = normal_to_2Dline(nodes[I[1]], nodes[I[2]]); r; end

@inline line_center(n1::AbstractVec,n2::AbstractVec) = 0.5*(n1+n2)
@inline function line_center(I::NTuple{2,Integer},nodes::AbstractArray{<:Vec2D,1}) ; @inbounds r = line_center(nodes[I[1]], nodes[I[2]]); return r; end


text:
function to_rna(dna)
    mapping = Dict('G' => 'C', 'C' => 'G', 'T' => 'A', 'A' => 'U')
    !issubset(Set(collect(dna)), keys(mapping)) && throw(ErrorException("Invalid DNA"))
    map(c->mapping[c], dna)
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
# utility functions for PressureDrop package

"""
`read_survey(<named arguments>)`

Reads in a wellbore deviation survey from a delimited file and returns a Wellbore object for use with pressure drop calculations.

Assumes the column order for the file is (MD, INC, TVD, <optional ID>), in U.S. field units, where:
    MD = measured depth, ft
    Inc = inclination, degrees from vertical
    TVD = true vertical depth, ft
    ID = pipe hydraulic diameter, inches

# Arguments

- `path::String`: absolute or relative path to survey file
- `delim::Char = ','`: file delimiter
- `skiplines::Int64 = 1`: number of lines to skip before survey data starts; assumes a 1-line header by default
- `maxdepth::Union{Bool, Real} = false`: If set to a real number, drop survey data past a certain measured depth. If false, keep the entire survey.
- `id_included::Bool = false`: whether the survey segment ID is stored as a fourth column. This is the easiest option to include tapered strings.
- `id::Real = 2.441`: the diameter to assume for the entire wellbore length, if the ID is not included in the survey file.
- `valves::Union{GasliftValves, Nothing} = nothing`: set of gas lift valves to add corresponding depths to the final Wellbore object via interpolation
- `allow_negatives::Bool = false`: allow negative depths on the survey inputs
"""
function read_survey(;path::String, delim::Char = ',', skiplines::Int64 = 1, maxdepth::Union{Bool, Real} = false, id_included::Bool = false, id::Real = 2.441, valves::Union{GasliftValves, Nothing} = nothing, allow_negatives::Bool = false)

    nlines = countlines(path) - skiplines
    ncols = id_included ? 4 : 3
    output = Array{Float64, 2}(undef, nlines, ncols)
    filestream = open(path, "r")

    try
        for skip in 1:skiplines
            readline(filestream)
        end

        for (index, line) in enumerate(eachline(filestream))
            parsedline = parse.(Float64, split(line, delim, keepempty = false))
            output[index, :] = parsedline
        end
    finally
        close(filestream)
    end

    if maxdepth != false
        output = output[output[:,1] .<= maxdepth, :]
    end

    if id_included
        return Wellbore(output[:,1], output[:,2], output[:,3], output[:,4], valves, allow_negatives)
    else
        return Wellbore(output[:,1], output[:,2], output[:,3], id, valves, allow_negatives)
    end
end

"""
`read_valves(;path::String, delim::Char = ',', skiplines::Int64 = 1)`

Expects a CSV with columns for [measured depth (ft)], [test rack opening pressure (psig)], [R-value (dimensionless)], [port size (diameter in 64ths inch)].

Indicate orifice valves with an R-value and PTRO of 0.

# Arguments
- `path::String`: absolute or relative path to survey file
- `delim::Char = ','`: file delimiter
- `skiplines::Int64 = 1`: number of lines to skip before survey data starts; assumes a 1-line header by default
"""
function read_valves(;path::String, delim::Char = ',', skiplines::Int64 = 1)

    nlines = countlines(path) - skiplines
    output = Array{Float64, 2}(undef, nlines, 4)
    filestream = open(path, "r")

    try
        for skip in 1:skiplines
            readline(filestream)
        end

        for (index, line) in enumerate(eachline(filestream))
            parsedline = parse.(Float64, split(line, delim, keepempty = false))
            output[index, :] = parsedline
        end
    finally
        close(filestream)
    end

    return GasliftValves(output[:,1], output[:,2], output[:,3], output[:,4]) #constructor will parse appropriately
end

"""
`interpolate(ref_array, property::Array{Real,1}, point)`

Interpolate between points without any bounds checking.
"""
function interpolate(ref_array::Array{T,1} where T <: Real, property::Array{T,1} where T <: Real, point)

    index_above = searchsortedlast(ref_array, point)
    index_below = index_above + 1
    x1 = ref_array[index_above]

    if x1 == point
        interpolated_value = property[index_above]
    else
        x2 = ref_array[index_below]
        y1, y2 = property[index_above], property[index_below]
        interpolated_value = y1 + (y2 - y1)/(x2 - x1) * (point - x1)
    end

    return interpolated_value
end

"""
`interpolate_all(well::Wellbore, properties::Array{Array{Real,1},1}, points::Array{Real,1})`

Interpolate multiple points for multiple properties with bounds checking.
"""
function interpolate_all(well::Wellbore, properties::Array{Array{T,1},1} where T <: Real, points::Array{T,1} where T <: Real)

    if !(all(length.(properties) .== length(well.md)))
        throw(DimensionMismatch("Property array lengths and number of wellbore survey points must match."))
    end

    if !(all(points .<= well.md[end]) && all(points .>= well.md[1]))
        throw(DimensionMismatch("Interpolation points cannot be outside wellbore measured depth."))
    end

    results = Array{Float64, 2}(undef, length(points), length(properties))
    for (index, property) in enumerate(properties)
        results[:, index] = map(pt -> interpolate(well.md, property, pt), points)
    end

    return results
end


text:
# # Linear regression

# We estimate simple linear regression model with a half-T prior.
# First, we load the packages we use.

using TransformVariables, LogDensityProblems, DynamicHMC, DynamicHMC.Diagnostics
using MCMCDiagnostics
using Parameters, Statistics, Random, Distributions
import ForwardDiff              # use for AD

# Then define a structure to hold the data: observables, covariates, and the degrees of
# freedom for the prior.

"""
Linear regression model ``y ∼ Xβ + ϵ``, where ``ϵ ∼ N(0, σ²)`` IID.

Flat prior for `β`, half-T for `σ`.
"""
struct LinearRegressionProblem{TY <: AbstractVector, TX <: AbstractMatrix,
                               Tν <: Real}
    "Observations."
    y::TY
    "Covariates"
    X::TX
    "Degrees of freedom for prior."
    ν::Tν
end

# Then make the type callable with the parameters *as a single argument*.

function (problem::LinearRegressionProblem)(θ)
    @unpack y, X, ν = problem   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    loglikelihood(Normal(0, σ), y .- X*β) + logpdf(TDist(ν), σ)
end

# We should test this, also, this would be a good place to benchmark and
# optimize more complicated problems.

N = 100
X = hcat(ones(N), randn(N, 2));
β = [1.0, 2.0, -1.0]
σ = 0.5
y = X*β .+ randn(N) .* σ;
p = LinearRegressionProblem(y, X, 1.0);
p((β = β, σ = σ))

# For this problem, we write a function to return the transformation (as it varies with the
# number of covariates).

function problem_transformation(p::LinearRegressionProblem)
    as((β = as(Array, size(p.X, 2)), σ = asℝ₊))
end

# Wrap the problem with a transformation, then use ForwardDiff for the gradient.

t = problem_transformation(p)
P = TransformedLogDensity(t, p)
∇P = ADgradient(:ForwardDiff, P);

# Finally, we sample from the posterior. `chain` holds the chain (positions and
# diagnostic information), while the second returned value is the tuned sampler
# which would allow continuation of sampling.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000);

# We use the transformation to obtain the posterior from the chain.

posterior = transform.(t, results.chain);

# Extract the parameter posterior means: `β`,

posterior_β = mean(first, posterior)

# then `σ`:

posterior_σ = mean(last, posterior)

# Effective sample sizes (of untransformed draws)

ess = vec(mapslices(effective_sample_size,
                    DynamicHMC.position_matrix(results.chain);
                    dims = 2))

# NUTS-specific statistics

summarize_tree_statistics(results.tree_statistics)


text:
################################################################################
#  Copyright 2021, Marta Vanin                                                 #
################################################################################
module SE_framework_paper_results

import DataFrames, CSV, Distributions, Query
import PowerModelsDistribution
import PowerModelsDistributionStateEstimation
import Plots: savefig, ylims!, xlims!
import Plots
import Random, Statistics

# Define Pkg constants
const _DF = DataFrames
const _DST = Distributions
const _PMD = PowerModelsDistribution
const _PMDSE = PowerModelsDistributionStateEstimation
const _RAN = Random
const _SEF = SE_framework_paper_results
const _STT = Statistics

include("plotting/all_plots.jl")

include("scripts/case_study_A.jl")
include("scripts/case_study_B.jl")
include("scripts/case_study_C.jl")
include("scripts/case_study_D.jl")
include("scripts/case_study_E.jl")

include("utils/linear_ivr_functions.jl")
include("utils/other_utils.jl")

export run_case_study_A, run_case_study_B, run_case_study_C, run_case_study_D, run_case_study_E
export plot_result_caseA,plot_result_caseB,plot_result_caseC,plot_result_caseD
export savefig # so users do not have to install plots to save their figures

end


text:
struct Nano_Day <: Dates.AbstractTime
    days::Day
    nanosecs::Nanosecond
end

function NanoDay(days::Day, nanosecs::Nanosecond)
    dy, ns = Dates.value(days), Dates.value(nanosecs)
    ddy, ns = fldmod(ns, NanosecondsPerDay)
    dy += ddy
    Day_Nano(Day(dy), Nanosecond(ns))
end

NanoDay(datetime::DateTime) =
    Nano_Day(Date(datetime).instant.periods, Time(datetime).instant)

NanoDay(date::Date) =
    Nano_Day(Date(datetime).instant.periods, Nanosecond0)

NanoDay(time::Time) =
    Nano_Day(Day0, time.instant)

NanoDay(nd::NanoDate) =
    Nano_Day(Date(nd).instant.periods, Time(nd).instant)

Base.:(+)(nd::Nano_Day, wk::Week) =
    Nano_Day(Day(Dates.value(nd.days) + Dates.value(wk) * DaysPerWeek), nd.nanosecs)

Base.:(+)(nd::Nano_Day, dy::Day) =
    Nano_Day(nd.days + dy, nd.nanosecs)

Base.:(+)(nd::Nano_Day, hr::Hour) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) + Dates.value(hr) * NanosecondsPerHour))

Base.:(+)(nd::Nano_Day, mi::Minute) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) + Dates.value(mi) * NanosecondsPerMinute))

Base.:(+)(nd::Nano_Day, sc::Second) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) + Dates.value(sc) * NanosecondsPerSecond))

Base.:(+)(nd::Nano_Day, ms::Millisecond) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) + Dates.value(ms) * NanosecondsPerMillisecond))

Base.:(+)(nd::Nano_Day, cs::Microsecond) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) + Dates.value(ms) * NanosecondsPerMicrosecond))

Base.:(+)(nd::Nano_Day, ns::Nanosecond) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) + Dates.value(ns)))

Base.:(-)(nd::Nano_Day, wk::Week) =
    Nano_Day(Day(Dates.value(nd.days) - Dates.value(wk) * DaysPerWeek), nd.nanosecs)

Base.:(-)(nd::Nano_Day, dy::Day) =
    Nano_Day(nd.days - dy, nd.nanosecs)

Base.:(-)(nd::Nano_Day, hr::Hour) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) - Dates.value(hr) * NanosecondsPerHour))

Base.:(-)(nd::Nano_Day, mi::Minute) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) - Dates.value(mi) * NanosecondsPerMinute))

Base.:(-)(nd::Nano_Day, sc::Second) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) - Dates.value(sc) * NanosecondsPerSecond))

Base.:(-)(nd::Nano_Day, ms::Millisecond) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) - Dates.value(ms) * NanosecondsPerMillisecond))

Base.:(-)(nd::Nano_Day, cs::Microsecond) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) - Dates.value(ms) * NanosecondsPerMicrosecond))

Base.:(-)(nd::Nano_Day, ns::Nanosecond) =
    NanoDay(nd.days, Nanosecond(Dates.value(nd.nanosecs) - Dates.value(ns)))


text:
"""
	solar(cmd0::String="", arg1=nothing; kwargs...)

Calculate and plot the day-night terminator and the civil, nautical and astronomical twilights.

Full option list at [`solar`]($(GMTdoc)solar.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_B)
- **C** | **formated** :: [Type => Bool]

    ($(GMTdoc)solar.html#c)
- **G** | **fill** :: [Type => Str | Number]

    ($(GMTdoc)solar.html#g)
- **I** | **sun** :: [Type => Bool | Tuple | NamedTuple]

    ($(GMTdoc)solar.html#i)
- $(GMT.opt_P)
- **M** | **dump** :: [Type => Bool]

    ($(GMTdoc)solar.html#m)
- $(GMT.opt_R)
- **N** | **invert** :: [Type => Bool]

    ($(GMTdoc)solar.html#n)
- **T** | **terminators** :: [Type => Bool | Tuple | NamedTuple]

    ($(GMTdoc)solar.html#t)
- **W** | **pen** :: [Type => Str | Tuple]

    ($(GMTdoc)solar.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bo)
- $(GMT.opt_h)
- $(GMT.opt_o)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function solar(cmd0::String="", arg1=nothing; first=true, kwargs...)

	arg2 = nothing		# May be needed if GMTcpt type is sent in via C
	N_args = (arg1 === nothing) ? 0 : 1

	gmt_proggy = (IamModern[1]) ? "solar "  : "pssolar "
	length(kwargs) == 0 && N_args == 0 && return monolitic(gmt_proggy, cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12cd/0d")
	cmd, = parse_common_opts(d, cmd, [:bo :c :h :o :p :t :UVXY :params], first)
	cmd  = parse_these_opts(cmd, d, [[:C :formated], [:M :dump], [:N :invert]])

	cmd  = add_opt_fill(cmd, d, [:G :fill], 'G')
	cmd  = add_opt(d, cmd, 'I', [:I :sun], (pos="",date="+d",TZ="+z"))
	cmd  = add_opt(d, cmd, 'T', [:T :terminators], (term="",date="+d",TZ="+z"))
	cmd *= opt_pen(d, 'W', [:W :pen])

	opt_extra = "";		do_finish = true
	if (occursin("-I", cmd) || occursin("-I", cmd0))
		opt_extra = "-I";		do_finish = false;	cmd = replace(cmd, opt_J => "")
	end
	cmd, K = finish_PS_nested(d, gmt_proggy * cmd, K, O, [:coast :logo :text])
	return finish_PS_module(d, cmd, opt_extra, K, O, do_finish, arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
solar!(cmd0::String="", arg1=nothing; kw...) = solar(cmd0, arg1; first=false, kw...)

const pssolar  = solar				# Alias
const pssolar! = solar!				# Alias


text:
struct HH3DPlaneWave{T,P} <: Functional
    direction::P
    wavenumber::T
    amplitude::T
end

function HH3DPlaneWave(direction, wavenumber; amplitude=1)
    w, a = promote(wavenumber, amplitude)
    HH3DPlaneWave(direction, w, a)
end

function (f::HH3DPlaneWave)(r)
    d = f.direction
    k = f.wavenumber
    a = f.amplitude
    a * exp(-im*k*dot(d,r))
end

struct NormalDerivative{F} <: Functional
    field::F
end

const ∂n = Val{:normalderivative}
(::Type{Val{:normalderivative}})(f) = NormalDerivative(f)

function (f::NormalDerivative{T})(manipoint) where T<:HH3DPlaneWave
    d = f.field.direction
    k = f.field.wavenumber
    a = f.field.amplitude
    n = normal(manipoint)
    r = cartesian(manipoint)
    -im*k*a * dot(d,n) * exp(-im*k*dot(d,r))
end

integrand(::NormalDerivative, test_vals, field_vals) = dot(test_vals[1], field_vals)


text:
# Copyright (c) 2021, Eric Sabo
# All rights reserved.
#
# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

include("linearcode.jl")

abstract type AbstractReedMullerCode <: AbstractLinearCode end

mutable struct ReedMullerCode <: AbstractReedMullerCode
    F::Union{FqNmodFiniteField, Nemo.GaloisField, AbstractAlgebra.GFField{Int64}}
    n::Integer
    k::Integer
    d::Union{Integer, Missing}
    r::Integer
    m::Integer
    G::Union{gfp_mat, fq_nmod_mat}
    Gorig::Union{gfp_mat, fq_nmod_mat, Missing}
    H::Union{gfp_mat, fq_nmod_mat}
    Horig::Union{gfp_mat, fq_nmod_mat, Missing}
    Gstand::Union{gfp_mat, fq_nmod_mat}
    Hstand::Union{gfp_mat, fq_nmod_mat}
    weightenum::Union{WeightEnumerator, Missing}
end

r(C::ReedMullerCode) = C.r
m(C::ReedMullerCode) = C.m

function show(io::IO, C::AbstractReedMullerCode)
    if get(io, :compact, false)
        println(io, "[$(length(C)), $(dimension(C)), $(minimumdistance(C))]_$(order(field(C))) Reed Muller code RM($(r(C)), $(m(C))).")
    else
        println(io, "[$(length(C)), $(dimension(C)), $(minimumdistance(C))]_$(order(field(C))) Reed Muller code RM($(r(C)), $(m(C))).")
        println(io, "Generator matrix: $(dimension(C)) × $(length(C))")
        for i in 1:dimension(C)
            print(io, "\t")
            for j in 1:length(C)
                if j != length(C)
                    print(io, "$(C.G[i, j]) ")
                elseif j == length(C) && i != dimension(C)
                    println(io, "$(C.G[i, j])")
                else
                    print(io, "$(C.G[i, j])")
                end
            end
        end
    end
end

function ReedMullergeneratormatrix(q::Integer, r::Integer, m::Integer)
    (0 ≤ r ≤ m) || error("Reed Muller codes require 0 ≤ r ≤ m, received r = $r and m = $m.")

    if q == 2
        F, _ = FiniteField(2, 1, "α")
        if r == m
            M = MatrixSpace(F, 2^m, 2^m)
            return M(1)
        elseif r == 0
            M = MatrixSpace(F, 1, 2^m)
            return M([1 for i in 1:2^m])
        else
            Grm1 = ReedMullergeneratormatrix(q, r, m - 1)
            Gr1m1 = ReedMullergeneratormatrix(q, r - 1, m - 1)
            M = MatrixSpace(F, size(Gr1m1, 1), size(Gr1m1, 2))
            return vcat(hcat(Grm1, Grm1), hcat(M(0), Gr1m1))
        end
    else
        error("Nonbinary Reed Muller codes have not yet been implemented.")
    end
end

function ReedMullerCode(q::Integer, r::Integer, m::Integer, verify::Bool=true)
    (1 ≤ r < m) || error("Reed Muller codes require 1 ≤ r < m, received r = $r and m = $m.")
    m < 64 || error("This Reed Muller code requires the implmentation of BigInts. Change if necessary.")
    q == 2 || error("Nonbinary Reed Muller codes have not yet been implemented.")

    G = ReedMullergeneratormatrix(q, r, m)
    H = ReedMullergeneratormatrix(q, m - r - 1, m)
    Gstand, Hstand = _standardform(G)

    if verify
        size(G, 2) == 2^m || error("Generator matrix computed in ReedMuller has the wrong number of columns; received: $(size(G, 2)), expected: $(BigInt(2)^m).")
        k = sum([binomial(m, i) for i in 0:r])
        size(G, 1) == k || error("Generator matrix computed in ReedMuller has the wrong number of rows; received: $(size(G, 1)), expected: $k.")
        if size(H) == (2^m - k, k)
            H = deepcopy(H')
        end
        for r in 1:size(Gstand, 1)
            iszero(Gstand[r, :] * H') || error("Column swap appeared in _standardform.")
        end
    end

    return ReedMullerCode(base_ring(G), size(G, 2), size(G, 1), 2^(m - r), r, m, G, missing,
        H, missing, Gstand, Hstand, missing)
end

function dual(C::ReedMullerCode)
    # return ReedMullerCode(C.q, C.m - C.r - 1, C.m)
    return ReedMullerCode(field(C), length(C), length(C) - dimension(C), 2^(r(C) + 1),
        m(C) - r(C) - 1, m(C), paritycheckmatrix(C), originalparitycheckmatrix(C), generatormatrix(C),
        originalgeneratormatrix(C), paritycheckmatrix(C, true), generatormatrix(C, true), missing)
end

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "hidapi"
version = v"0.10.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/libusb/hidapi.git", "f6d0073fcddbdda24549199445e844971d3c9cef")
]

# Script template (structure is the same for all platforms, directory for make differs)
# The additional eudev-dev package is to allow the ./configure step. For linux there are
# two backends: linux and libusb. This script builds only the 'libusb' backend
# which doesn't needs libudev...
script = raw"""
cd $WORKSPACE/srcdir/hidapi/

if [[ ${target} == *-linux-* ]]; then
    # We need this only to trick `configure` into thinking that udev is available,
    # but we aren't going to actually use it
    apk add eudev-dev
fi

./bootstrap
./configure --prefix=${prefix} --host=${target} --build=${MACHTYPE}
# Install all license files
install_license LICENSE*.txt

if [[ ${target} == *-mingw* ]]; then
    cd windows
elif [[ ${target} == *-apple-* ]]; then
    cd mac
else
    cd libusb
fi

make -j${nproc}
make install
"""

products = [
    LibraryProduct(["libhidapi","libhidapi-libusb"], :hidapi)
]

dependencies = [
    Dependency("libusb_jll")
]

platforms = [p for p in supported_platforms() if !Sys.isfreebsd(p)]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module CatBoost

using PyCall
using Conda
using DataFrames
using OrderedCollections
using Tables

#####
##### Exports
#####

export catboost

export CatBoostRegressor, CatBoostClassifier
export fit!, cv, predict, predict_proba
export Pool
export pandas_to_df

# Datasets API.
export load_dataset

#####
##### _init_
#####

const catboost = PyNULL()
const catboost_datasets = PyNULL()
const pandas = PyNULL()

function load_python_deps!()
    copy!(catboost, pyimport("catboost"))
    copy!(catboost_datasets, pyimport("catboost.datasets"))
    copy!(pandas, pyimport("pandas"))
    return nothing
end

function __init__()
    try
        load_python_deps!()
    catch ee
        if PyCall.conda
            Conda.pip_interop(true)
            Conda.pip("install", ["catboost", "pandas"])
            load_python_deps!()
        else
            typeof(ee) <: PyCall.PyError || rethrow(ee)
            @warn("""
                 Python Dependencies not installed
                 Please either:
                 - Rebuild PyCall to use Conda, by running in the julia REPL:
                 - `ENV["PYTHON"]=""; Pkg.build("PyCall"); Pkg.build("CatBoost")`
                 - Or install the depencences, eg by running pip
                 - `pip install catboost pandas`
                 """)
        end
    end

    @doc """
        cv(pool::PyObject; kwargs...) -> DataFrame

    Accepts a [`CatBoost.Pool`](@ref) positional argument to specify the training data,
    and keyword arguments to configure the settings. See the python documentation below
    for what keyword arguments are accepted.

    ---

    ## Python documentation for `catboost.cv`

    $(@doc catboost.cv)
    """
    cv

    return nothing
end

#####
##### Data
#####

"""
    Pool(data; label=nothing, cat_features=nothing, text_features=nothing,
         pairs=nothing, delimiter='\t', has_header=false, weight=nothing,
         group_id = nothing, group_weight=nothing, subgroup_id=nothing,
         pairs_weight=nothing, baseline=nothing, features_names=nothing,
         thread_count = -1) -> PyObject

Creates a `Pool` object holding training data and labels. `data` may also be passed
as a keyword argument.

"""
function Pool(data; kwargs...)
    return catboost.Pool(to_catboost(data);
                         (k => to_catboost(v) for (k, v) in pairs(kwargs))...)
end
Pool(; kwargs...) = catboost.Pool(; (k => to_catboost(v) for (k, v) in pairs(kwargs))...)

#####
##### Models
#####

CatBoostRegressor(args...; kwargs...) = catboost.CatBoostRegressor(args...; kwargs...)
CatBoostClassifier(args...; kwargs...) = catboost.CatBoostClassifier(args...; kwargs...)

fit!(cbm, args...; kwargs...) = cbm.fit(all_to_catboost(args)...; kwargs...)
predict(cbm, args...; kwargs...) = cbm.predict(all_to_catboost(args)...; kwargs...)
function predict_proba(cbm, args...; kwargs...)
    return cbm.predict_proba(all_to_catboost(args)...; kwargs...)
end

#####
##### Cross validation
#####

cv(pool::PyObject; kwargs...) = pandas_to_df(catboost.cv(pool; kwargs...))

#####
##### Conversion utilities
#####

"""
    to_catboost(arg)

`to_catboost` is called on each argument passed to [`fit`](@ref), [`predict`](@ref), [`predict_proba`](@ref), and [`cv`](@ref)
to allow customization of the conversion of Julia types to python types. If `to_catboost` emits a Julia type, then
PyCall will try to convert it appropriately (automatically).

By default, `to_catboost` simply checks if the argument satisfies `Tables.istable(arg)`, and if so, it outputs
a corresponding pandas table, and otherwise passes it on.

To customize the conversion for custom types, provide a method for this function.
"""
to_catboost(arg) = Tables.istable(arg) ? to_pandas(arg) : arg

# utility for calling `to_catboost` on each argument of a function
all_to_catboost(args) = (to_catboost(arg) for arg in args)

# the Julia-side code does not copy the columns, but the `pandas.DataFrame`
# constructor seems to make a copy here. Maybe that can be avoided?
function to_pandas(tbl)
    # ensure we have a column table
    col_table = Tables.columns(tbl)
    # write it in a way that pandas will understand (after PyCall conversion)
    dict_table = OrderedDict(col => Tables.getcolumn(col_table, col)
                             for col in Tables.columnnames(col_table))
    return pandas.DataFrame(; data=dict_table)
end

function pandas_to_df(pandas_df::PyObject)
    colnames = map(pandas_df.columns) do c
        ret = c isa PyObject ? PyAny(c) : c
        return ret isa Int ? ret + 1 : ret
    end
    df = DataFrame(Any[Array(getproperty(pandas_df, c).values) for c in colnames],
                   map(Symbol, colnames))
    return df
end

#####
##### Datasets API
#####

function load_dataset(dataset_name::Symbol)
    train, test = getproperty(catboost_datasets, dataset_name)()
    return train, test
end

end # module


text:
function does_data_file_exist(path_to_data_file::String)::Bool

    if (ispath(path_to_data_file) == true)
        return true
    end

    return false
end

max_stars_repo_name: Pooksoft/BinomialEquityAssetPriceModel

max_issues_repo_name: Pooksoft/BinomialEquityAssetPriceModel

max_forks_repo_name: Pooksoft/BinomialEquityAssetPriceModel


text:
using SparseArrays

# Dummy functions to allow not specifying these characteristics
function F_dummy(ẋ, xc, xd, parms, t)
	fill!(ẋ, 0)
	nothing
end

# Dummy flow to be used in rejection algorithm
function Phi_dummy(out, xc, xd, parms, t)
	# vector field used for the continuous variable
	# trivial dynamics
	out[1,:] .= xc
	out[2,:] .= xc
	nothing
end

mutable struct PDMPJumpTime{Tc <: Real, Td}
	tstop_extended::Tc
	lastjumptime::Tc
	njumps::Td

	# fields required for the rejection method
	lambda_star::Tc					# bound on the total rate
	ppf::Vector{Tc}
	reject::Bool					# boolean to know whether to reject or not the step
	fictitous_jumps::Td
end

struct PDMPCaracteristics{TF, TR, TJ, vecc, vecd, vecrate, Tparms}
	F::TF						# vector field for ODE between jumps
	R::TR			    		# rate function for jumps
	pdmpjump::TJ

	xc::vecc					# current continuous variable
	xd::vecd					# current discrete variable

	xc0::vecc					# initial continuous variable
	xd0::vecd					# initial discrete variable
	ratecache::vecrate			# to hold the rate vector for inplace computations. Also used to initialise rate as this can be an issue for StaticArrays.jl
	parms::Tparms				# container to hold parameters to be passed to F, R, Delta

	function PDMPCaracteristics(F, R, Delta, nu::Tnu, xc0::vecc, xd0::vecd, parms::Tparms) where {Tc, Td, Tparms, Tnu <: AbstractMatrix{Td},
						vecc <: AbstractVector{Tc},
						vecd <: AbstractVector{Td}}
		jump = Jump(nu, Delta)
		rate = dualcache(get_rate_prototype(jump, Tc))
		ratefunction = VariableRate(R)
		return new{typeof(F), typeof(ratefunction), typeof(jump), vecc, vecd, typeof(rate), Tparms}(F, ratefunction, jump, copy(xc0), copy(xd0), copy(xc0), copy(xd0), rate, parms)
	end

	function PDMPCaracteristics(F, R::TR, Delta, nu::Tnu, xc0::vecc, xd0::vecd, parms::Tparms) where {Tc, Td, Tparms, Tnu <: AbstractMatrix{Td},
						vecc <: AbstractVector{Tc},
						vecd <: AbstractVector{Td},
						TR <: AbstractRate}
		jump = Jump(nu, Delta)
		rate = dualcache(get_rate_prototype(jump, Tc))
		return new{typeof(F), typeof(R), typeof(jump), vecc, vecd, typeof(rate), Tparms}(F, R, jump, copy(xc0), copy(xd0), copy(xc0), copy(xd0), rate, parms)
	end
end

function PDMPCaracteristics(F, R, nu::Tnu, xc0::vecc, xd0::vecd, parms::Tparms) where {Tc, Td, Tparms, Tnu <: AbstractMatrix{Td}, vecc <: AbstractVector{Tc}, vecd <: AbstractVector{Td}}
	return PDMPCaracteristics(F, R, Delta_dummy, nu, xc0, xd0, parms)
end

function init!(pb::PDMPCaracteristics)
	pb.xc .= pb.xc0
	pb.xd .= pb.xd0
	init!(pb.R)
end

"""
	PDMPProblem(F, R, Delta, nu, xc0, xd0, p, tspan)

Create a PDMP problem to be simulated. To define a PDMP Problem, you first need to give the function `F` and the initial condition `xc0` which define an ODE: dxc/dt = F(xc(t),xd(t),p,t). Jumps are defined as a Jump process which changes states at some rate `R`. Note, that in between jumps, `xd` is constant but `xc` is allowed to evolve.

## Arguments
- `F`: inplace function `F(du, u, p, t)` representing the vector field
- `R`: the function to compute the transition rates. It should be specified in-place as `R(rate::AbstractVector, xc, xd, p, t, issum::Bool)` where it mutates `rate`. Note that a boolean `issum` is provided and the behavior of `R` should be as follows
    1. if `issum == true`, we only require `R` to return the total rate, *e.g.* `sum(rate)`. We use this formalism because sometimes you can compute the `sum` without mutating `rate`.
    2. if `issum == false`, `R` must populate `rate` with the updated rates

We then need to provide the way the jumps affect the state variable. There are two possible ways here:
    1. either give a transition matrix `nu`: it will only affect the discrete component `xd` and leave `xc` unaffected.
    2. give a function to implement jumps `Delta(xc, xd, parms, t, ind_reaction::Int64)` where you can mutate `xc,xd` or `parms`. The argument `ind_reaction` is the index of the reaction at which the jump occurs. See `examples/pdmp_example_eva.jl` for an example.
- `Delta` [Optional]: the function to effect the jumps
- `nu` [Optional]: the transition matrix
- `xc0`: the initial condition of the continuous part
- `xd0`: the initial condition of the discrete part
- `p`: the parameters to be provided to the functions `F, R, Delta`
- `tspan`: The timespan for the problem.

## Constructors

- `PDMPProblem(F,R,Delta,nu,xc0,xd0,p,tspan)`
- `PDMPProblem(F,R,nu,xc0,xd0,p,tspan)` when ones does not want to provide the function `Delta`
- `PDMPProblem(F,R,Delta,reaction_number::Int64,xc0,xd0,p,tspan)` when ones does not want to provide the transition matrix `nu`. The length `reaction_number` of the rate vector must then be provided.

We also provide a wrapper to [DiffEqJump.jl](https://github.com/JuliaDiffEq/DiffEqJump.jl). This is quite similar to how a `JumpProblem` would be created.

- `PDMPProblem(prob, jumps...)` where `prob` can be an `ODEProblem`. For an example, please consider `example/examplediffeqjumpwrapper.jl`.
"""
struct PDMPProblem{Tc, Td, vectype_xc <: AbstractVector{Tc},
						vectype_xd <: AbstractVector{Td},
						Tcar}
	tspan::Vector{Tc}				    			# final simulation time interval, we use an array to be able to mutate it
	simjptimes::PDMPJumpTime{Tc, Td}				# space to save result
	time::Vector{Tc}
	Xc::VectorOfArray{Tc, 2, Array{vectype_xc, 1}}	# continuous variable history
	Xd::VectorOfArray{Td, 2, Array{vectype_xd, 1}}	# discrete variable history
	# variables for debugging
	rate_hist::Vector{Tc}							# to save the rates for debugging purposes
	caract::Tcar									# struct for characteristics of the PDMP
end

pushTime!(pb::PDMPProblem, t) = push!(pb.time, t)
pushXc!(pb::PDMPProblem, xc) = push!(pb.Xc, xc)
pushXd!(pb::PDMPProblem, xd) = push!(pb.Xd, xd)

function init!(pb::PDMPProblem)
	init!(pb.caract)
	pb.simjptimes.tstop_extended = -log(rand())
	pb.simjptimes.lastjumptime = pb.tspan[1]
	pb.simjptimes.njumps = 0
	pb.simjptimes.fictitous_jumps = 0
	resize!(pb.time, 1)
	resize!(pb.rate_hist, 1)
	resize!(pb.Xc.u, 1)
	resize!(pb.Xd.u, 1)
end

# callable struct used in the iterator interface
function (prob::PDMPProblem)(u, t, integrator)
	t == prob.simjptimes.tstop_extended
end

function PDMPProblem(F::TF, R::TR, DX::TD, nu::Tnu,
				xc0::vecc, xd0::vecd, parms::Tp,
				tspan) where {Tc, Td, Tnu <: AbstractMatrix{Td}, Tp, TF ,TR ,TD, vecc <: AbstractVector{Tc}, vecd <:  AbstractVector{Td}}
	ti, tf = tspan
	caract = PDMPCaracteristics(F, R, DX, nu, xc0, xd0, parms)
	return PDMPProblem{Tc, Td, vecc, vecd, typeof(caract)}(
			[ti, tf],
			PDMPJumpTime{Tc, Td}(Tc(0), ti, 0, Tc(0), Vector{Tc}([0, 0]), false, 0),
			[ti],
			VectorOfArray([copy(xc0)]), VectorOfArray([copy(xd0)]),
			Tc[],
			caract)
end

function PDMPProblem(F, R, nu::Tnu, xc0::vecc, xd0::vecd, parms,
				tspan; kwargs...) where {Tc, Td, Tnu <: AbstractMatrix{Td}, vecc <: AbstractVector{Tc}, vecd <:  AbstractVector{Td}}
	return PDMPProblem(F, R, Delta_dummy, nu, xc0, xd0, parms, tspan; kwargs...)
end

function PDMPProblem(F, R, Delta, reaction_number::Int64, xc0::vecc, xd0::vecd, parms,
				tspan; kwargs...) where {Tc, Td, vecc <: AbstractVector{Tc}, vecd <:  AbstractVector{Td}}
	return PDMPProblem(F, R, Delta, spzeros(Int64, reaction_number, length(xd0)), xc0, xd0, parms, tspan; kwargs...)
end


text:
using LinearAlgebra: I, norm

@testset "I" begin
    @testset "ident" begin
        @test SPICE._ident() == Array{Float64}(I, 3, 3)
    end
    @testset "ilumin" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            et = str2et("2007 FEB 3 00:00:00.000")
            trgepc, obspos, trmpts = edterm("UMBRAL", "SUN", "MOON", et, "IAU_MOON", "LT+S", "EARTH", 3)
            expected_trgepc = 223732863.86351672
            expected_obspos = [394721.1024056578753516078,
                               27265.11780063395417528227,
                               -19069.08478859506431035697]
            expected_trmpts = [-1.53978381936825627463e+02,
                               -1.73056331949840728157e+03,
                               1.22893325627419600088e-01,
                               87.37506200891714058798,
                               864.40670594653545322217,
                               1504.56817899807947469526,
                               42.213243378688254,
                               868.21134651980412,
                               -1504.3223922609538]
            @test trgepc ≈ expected_trgepc
            @test obspos ≈ expected_obspos
            for i in eachindex(trmpts, expected_trmpts)
                @test trmpts[i] ≈ expected_trmpts[i]
            end
            iluet0, srfvec0, phase0, solar0, emissn0 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,1], abcorr="LT+S")
            @test rad2deg(solar0) ≈ 90.269765819
            iluet1, srfvec1, phase1, solar1, emissn1 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,2], abcorr="LT+S")
            @test rad2deg(solar1) ≈ 90.269765706
            iluet2, srfvec2, phase2, solar2, emissn2 = ilumin("Ellipsoid", "MOON", et, "IAU_MOON",
                                                              "EARTH", trmpts[:,3], abcorr="LT+S")
            @test rad2deg(solar2) ≈ 90.269765730
        finally
            kclear()
        end
    end
    @testset "illumf" begin
        try
            furnsh(path(CORE, :lsk),
                   path(CORE, :pck),
                   path(CORE, :spk),
                   path(CASSINI, :pck),
                   path(CASSINI, :sat_spk),
                   path(CASSINI, :tour_spk),
                   path(CASSINI, :fk),
                   path(CASSINI, :ck),
                   path(CASSINI, :sclk),
                   path(CASSINI, :ik))
            et = str2et("2013 FEB 25 11:50:00 UTC")
            camid = bodn2c("CASSINI_ISS_NAC")
            shape, obsref, bsight, bounds = getfov(camid, 4)
            # run sincpt on boresight vector
            spoint, etemit, srfvec = sincpt("ELLIPSOID", "ENCELADUS", et, "IAU_ENCELADUS", "CN+S",
                                            "CASSINI", obsref, bsight)
            trgepc2, srfvec2, phase, incid, emissn, visible, lit = illumf("ELLIPSOID", "ENCELADUS",
                                                                          "SUN", et, "IAU_ENCELADUS",
                                                                          "CN+S", "CASSINI", spoint)
            @test rad2deg(phase) ≈ 161.82854377660345
            @test rad2deg(incid) ≈ 134.92108561449996
            @test rad2deg(emissn) ≈ 63.23618556218115
            @test !lit # Incidence angle is greater than 90deg
            @test visible # Emission angle is less than 90deg
            @test_throws SpiceError illumf("ELLIPSOID", "norbert", "SUN", et, "IAU_ENCELADUS", "CN+S",
                                           "CASSINI", spoint)
            @test_throws ArgumentError illumf("ELLIPSOID", "ENCELADUS", "SUN", et, "IAU_ENCELADUS", "CN+S",
                                              "CASSINI", spoint[1:2])
        finally
            kclear()
        end
    end
    @testset "illumg" begin
        try
            furnsh(path(CORE, :lsk),
                   path(CORE, :pck),
                   path(CORE, :spk),
                   path(CASSINI, :pck),
                   path(CASSINI, :sat_spk),
                   path(CASSINI, :tour_spk),
                   path(CASSINI, :fk),
                   path(CASSINI, :ck),
                   path(CASSINI, :sclk),
                   path(CASSINI, :ik))
            et = str2et("2013 FEB 25 11:50:00 UTC")
            spoint, trgepc, srfvec = subpnt("NEAR POINT/ELLIPSOID", "ENCELADUS", et, "IAU_ENCELADUS",
                                            "CN+S", "EARTH")
            trgepc2, srfvec2, phase, incid, emissn = illumg("ELLIPSOID", "ENCELADUS", "SUN", et,
                                                            "IAU_ENCELADUS", "CASSINI", spoint, abcorr="CN+S")

            @test rad2deg(phase) ≈ 161.859925246638
            @test rad2deg(incid) ≈ 18.47670084384343
            @test rad2deg(emissn) ≈ 143.6546170649875
        finally
            kclear()
        end
    end
    @testset "inedpl" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            time = "Oct 31 2002, 12:55:00 PST"
            et = str2et(time)
            state, ltime = spkezr("EARTH", et, "J2000", "LT+S", "SUN")
            pos = state[1:3]
            radii = bodvrd("EARTH", "RADII", 3)
            pos = [pos[1] / radii[1]^2.0,
                   pos[2] / radii[2]^2.0,
                   pos[3] / radii[3]^2.0]
            plane = nvc2pl(pos, 1.0)
            term = inedpl(radii[1], radii[2], radii[3], plane)
            expected_center = [0.21512031, 0.15544527, 0.067391641]
            expected_smajor = [-3.73561164720596843836e+03, 5.16970328302375583007e+03,
                               1.35988201424391742850e-11]
            expected_sminor = [-1276.33357469839393161237, -922.27470443423590040766,
                               6159.97371233560443215538]
            @test center(term) ≈ expected_center atol=1e-6
            @test semi_major(term) ≈ expected_smajor
            @test semi_minor(term) ≈ expected_sminor
            @test norm(semi_major(term)) ≈ 6378.1365 rtol=1e-4
            @test norm(semi_minor(term)) ≈ 6358.0558 rtol=1e-4
            @test_throws SpiceError inedpl(-radii[1], radii[2], radii[3], plane)
        finally
            kclear()
        end
    end
    @testset "inelpl" begin
        try
            furnsh(path(CORE, :pck))
            radii = bodvrd("SATURN", "RADII", 3)
            vertex = [100.0 * radii[1], 0.0, radii[1] * 100.0]
            limb = edlimb(radii..., vertex)
            normal = [0.0, 0.0, 1.0]
            point = [0.0, 0.0, 0.0]
            plane = nvp2pl(normal, point)
            nxpts, xpt1, xpt2 = inelpl(limb, plane)
            expected_xpt1 = [602.68000, 60264.9865, 0.0]
            expected_xpt2 = [602.68000, -60264.9865, 0.0]
            @test nxpts == 2.0
            @test expected_xpt1 ≈ xpt1
            @test expected_xpt2 ≈ xpt2
        finally
            kclear()
        end
    end
    @testset "inrypl" begin
        try
            furnsh(path(CORE, :pck))
            radii = bodvrd("SATURN", "RADII", 3)
            vertex = [3.0 * radii[1], 0.0, radii[3] * 0.5]
            dire = [0.0, cos(deg2rad(30.0)), -1.0 * sin(deg2rad(30.0))]
            normal = [0.0, 0.0, 1.0]
            point = [0.0, 0.0, 0.0]
            plane = nvp2pl(normal, point)
            nxpts, xpt = inrypl(vertex, dire, plane)
            expected_xpt = [180804.0, 47080.6050513, 0.0]
            @test nxpts == 1
            @test xpt ≈ expected_xpt
        finally
            kclear()
        end
    end
    @testset "insrtc" begin
        cell = SpiceCharCell(10, 10)
        clist = ["aaa", "bbb", "ccc", "bbb"]
        for c in clist
            insrtc!(cell, c)
        end
        @test collect(cell) == ["aaa", "bbb", "ccc"]
    end
    @testset "insrtd" begin
        cell = SpiceDoubleCell(8)
        dlist = [0.5, 2.0, 30.0, 0.01, 30.0]
        for d in dlist
            insrtd!(cell, d)
        end
        @test collect(cell) == [0.01, 0.5, 2.0, 30.0]
    end
    @testset "insrti" begin
        cell = SpiceIntCell(8)
        ilist = [1, 2, 30, 1, 30]
        for i in ilist
            insrti!(cell, i)
        end
        @test collect(cell) == [1, 2, 30]
    end
    @testset "inter" begin
        i1 = SpiceIntCell(8)
        i2 = SpiceIntCell(8)
        insrti!(i1, 1)
        insrti!(i1, 2)
        insrti!(i2, 1)
        insrti!(i2, 3)
        out = inter(i1, i2)
        @test out == [1]

        d1 = SpiceDoubleCell(8)
        d2 = SpiceDoubleCell(8)
        insrtd!(d1, 1.0)
        insrtd!(d1, 2.0)
        insrtd!(d2, 1.0)
        insrtd!(d2, 3.0)
        out = inter(d1, d2)
        @test out == [1.0]

        c1 = SpiceCharCell(8, 8)
        c2 = SpiceCharCell(8, 8)
        insrtc!(c1, "1")
        insrtc!(c1, "2")
        insrtc!(c2, "1")
        insrtc!(c2, "3")
        out = inter(c1, c2)
        @test out == ["1"]
    end
	@testset "intmax" begin
		@test SPICE._intmax() == typemax(SPICE.SpiceInt)
	end
	@testset "intmin" begin
		@test SPICE._intmin() == typemin(SPICE.SpiceInt)
	end
    @testset "invert" begin
        m1 = [0.0 -1.0 0.0; 0.5 0.0 0.0; 0.0 0.0 1.0]
        expected = [0.0 2.0 0.0; -1.0 0.0 0.0; 0.0 0.0 1.0]
        mout = SPICE._invert(m1)
        @test expected ≈ mout
        @test inv(m1) ≈ mout
    end
    @testset "invort" begin
        m1 = eul2m(3, 2, 1, 3, 2, 1)
        @test inv(m1) ≈ SPICE._invort(m1)
    end
    @testset "isordv" begin
        @test isperm([1, 2]) == SPICE._isordv([1, 2])
        @test isperm([1, 2, 3]) == SPICE._isordv([1, 2, 3])
        @test isperm([1, 2, 3, 4]) == SPICE._isordv([1, 2, 3, 4])
        @test isperm([1, 1, 1]) == SPICE._isordv([1, 1, 1])
    end
    @testset "isrchc" begin
        array = ["1", "0", "4", "2"]
        @test SPICE._isrchc("4", array) == findfirst(array .== "4")
        @test SPICE._isrchc("2", array) == findfirst(array .== "2")
    end
    @testset "isrchd" begin
        array = [1.0, 0.0, 4.0, 2.0]
        @test SPICE._isrchd(4.0, array) == findfirst(array .== 4.0)
        @test SPICE._isrchd(2.0, array) == findfirst(array .== 2.0)
    end
    @testset "isrchi" begin
        array = [1, 0, 4, 2]
        @test SPICE._isrchi(4, array) == findfirst(array .== 4)
        @test SPICE._isrchi(2, array) == findfirst(array .== 2)
    end
    @testset "isrot" begin
        m = eul2m(3, 2, 1, 3, 2, 1)
        @test isrot(m, 0.0001, 0.0001)
        @test_throws ArgumentError isrot(m[1:2, 1:2], 0.0001, 0.0001)
        @test_throws SpiceError isrot(m, -0.0001, 0.0001)
    end
    @testset "iswhsp" begin
        @test isempty(strip("   ")) == SPICE._iswhsp("   ")
        @test isempty(strip("spice")) == SPICE._iswhsp("spice")
    end
end


text:
# code from https://github.com/MikeInnes/BrainForth/blob/master/src/brainfuck.jl
using Sockets

mutable struct Tape
  count::Int
  pos::Int
  tape::Vector{UInt8}
end

Tape() = Tape(0, 1, [0])

Base.getindex(t::Tape) = t.tape[t.pos]
Base.setindex!(t::Tape, v) = t.tape[t.pos] = v

function Base.show(io::IO, t::Tape)
  print(io, "[$(t.count)] ")
  for i = 1:length(t.tape)
    print(io, t.tape[i], i == t.pos ? "* " : " ")
  end
end

Base.:(==)(a::Tape, b::Tape) = a.tape == b.tape

function left!(t::Tape)
  t.pos == length(t.tape) && t.tape[end] == 0 && pop!(t.tape)
  t.pos == 1 ? pushfirst!(t.tape, 0) : (t.pos -= 1)
  return
end

function right!(t::Tape)
  t.pos == 1 && t.tape[1] == 0 && (popfirst!(t.tape); t.pos -= 1)
  t.pos == length(t.tape) && push!(t.tape, 0)
  t.pos += 1
  return
end

clip(n) = n > 255 ? n - 256 : n < 0 ? n + 256 : n
inc!(t::Tape) = (t[] = clip(t[] + 1))
dec!(t::Tape) = (t[] = clip(t[] - 1))

function read!(t::Tape, io::IO)
  t[] = read(io, UInt8)
end

function write!(t::Tape, io::IO)
  write(io, t[])
end

# Gets ~370 MHz

function interpret(t::Tape, bf; input::IO = stdin, output::IO = stdout)
  loops = Int[]
  scan = 0
  ip = 1
  @inbounds while ip <= length(bf)
    t.count += 1
    op = bf[ip]
    if op == '['
      scan > 0 || t[] == 0 ? (scan += 1) :
      push!(loops, ip)
    elseif op == ']'
      scan > 0 ? (scan -= 1) :
      t[] == 0 ? pop!(loops) :
      (ip = loops[end])
    elseif scan == 0
      op == '+' ? inc!(t) :
      op == '-' ? dec!(t) :
      op == '<' ? left!(t) :
      op == '>' ? right!(t) :
      op == ',' ? read!(t, input) :
      op == '.' ? write!(t, output) :
      op == '#' ? println(t) :
        nothing
    end
    ip += 1
  end
  return t
end

interpret(t::Tape, bf::String; kws...) = interpret(t, collect(bf); kws...)

interpret(bf; kws...) = interpret(Tape(), bf; kws...)

function main(text)
  interpret(text)
end

println("JIT warming up")
main(">++[<+++++++++++++>-]<[[>+>+<<-]>[<+>-]++++++++[>++++++++<-]>[-]<<>++++++++++[>++++++++++[>++++++++++[>++++++++++[>++++++++++[>++++++++++[>++++++++++[-]<-]<-]<-]<-]<-]<-]<-]++++++++++")

println("bench")

text = open(ARGS[1]) do file
  read(file, String)
end

function notify(msg)
  try
    socket = connect("localhost", 9001)
    write(socket, msg)
    close(socket)
  catch
    # standalone usage
  end
end

notify("Julia\t$(getpid())")

x = @timed main(text)
println("Elapsed: $(x[2]), Allocated: $(x[3]), GC Time: $(x[4])")

notify("stop")

max_forks_repo_name: Kreyren/benchmarks


text:
Base.sizeof(t::RandomProjectionTree) = sizeof(t.root)

Base.sizeof(n::RandomProjectionTreeNode) = sizeof(n.indices) + sizeof(n.isleaf) +
                                           sizeof(n.hyperplane) + sizeof(n.offset) +
                                           sizeof(n.leftchild) + sizeof(n.rightchild)

num_nodes(t::RandomProjectionTree) = num_nodes(t.root)

function num_nodes(n::RandomProjectionTreeNode)
    n.isleaf || true + num_nodes(n.leftchild) + num_nodes(n.rightchild)
end

num_leaves(t::RandomProjectionTree) = num_leaves(t.root)

function num_leaves(n::RandomProjectionTreeNode)
    n.isleaf || num_leaves(n.leftchild) + num_leaves(n.rightchild)
end

leaves(t::RandomProjectionTree) = leaves(t.root)

function leaves(n::RandomProjectionTreeNode)
    n.isleaf && return [n]
    return append!(leaves(n.leftchild), leaves(n.rightchild))
end

max_nnz(t::RandomProjectionTree) = max_nnz(t.root)

function max_nnz(n::RandomProjectionTreeNode)
    n.isleaf && return 0
    return max(nnz(n.hyperplane), max_nnz(n.leftchild), max_nnz(n.rightchild))
end


text:
using Manifolds, Manopt, Test

@testset "Subgradient Plan" begin
    M = Euclidean(2)
    x = [4.0, 2.0]
    x0 = [5.0, 2.0]
    o = SubGradientMethodOptions(M, x0, StopAfterIteration(200), ConstantStepsize(1.0))
    o.∂ = [1.0, 0.0]
    f = y -> distance(M, y, x)
    ∂f = y -> if distance(M, x, y) == 0
        zero_tangent_vector(M, y)
    else
        -2 * log(M, y, x) / max(10 * eps(Float64), distance(M, x, y))
    end
    p = SubGradientProblem(M, f, ∂f)
    oR = solve(p, o)
    xHat = get_solver_result(oR)
    @test get_initial_stepsize(p, o) == 1.0
    @test get_stepsize(p, o, 1) == 1.0
    @test get_last_stepsize(p, o, 1) == 1.0
    # Check Fallbacks of Problen
    @test get_cost(p, x) == 0.0
    @test norm(M, x, get_subgradient(p, x)) == 0
    @test_throws ErrorException get_gradient(p, o.x)
    @test_throws ErrorException get_proximal_map(p, 1.0, o.x, 1)
    o2 = subgradient_method(M, f, ∂f, copy(x0); return_options=true)
    xhat2 = get_solver_result(o2)
    @test f(xhat2) <= f(x0)
end


text:
@testset "Test Boltzmann solvers" begin
    include(joinpath("BoltzmannSolver", "test_damping_1v.jl"))
end


text:
module Amphitheatre

using MeshCat
using GeometryTypes
using Colors
using PlotUtils
using Sockets: @ip_str, IPAddr, IPv4, IPv6
#
using Mongoc
using GraffSDK
using JSON2
using Base64
# using Caesar
using RoME

#
using DocStringExtensions
# using Rotations
using TransformUtils
using CoordinateTransformations, Rotations
using StaticArrays
# using FileIO
# using Base64

const CTs = CoordinateTransformations
const TUs = TransformUtils

export
	BasicFactorGraphPose,
	AbstractAmphitheatre,
	visualize,
	Visualizer,
	startMeshCatVisualizer,
	stopAmphiVis!,
	#Common visse
	visPose!,
	visPoint!,
	#camera
	CameraModel,
	#DepthImages
	cloudFromDepthImageClampZ,
	cloudFromDepthImage,
	#Amphis
	TagkTl,
	AbstractAmphitheatre,
	plDrawProp,
	BasicFactorGraphPose,
	BasicGraffPose,
	TagOnPose,
	PCloudFactorGraphPose,
	GraffTagOnPose,
	GraffCloudOnPose,
	#re-export
	RGBA,
	ColorGradient, cgrad,
	@ip_str, IPAddr, IPv4, IPv6

include("../Common/CameraModel.jl")
include("../Common/DepthImages.jl")

include("common.jl")
include("amphis.jl")
include("pointclouds.jl")
include("pointcloudAmphis.jl")
include("reprojectAmphis.jl")

global runAmphi = true

function stopAmphiVis!()
	global runAmphi
  	runAmphi = false
  	nothing
end

"""
    $(SIGNATURES)
Initialize empty visualizer window with home axis.  New browser window will be opened based on `start_browser=true`.
"""
function startMeshCatVisualizer(;host=ip"127.0.0.1",
								 start_browser::Bool=true,
                                 draworigin::Bool=true,
                                 originscale::Float64=1.0)

    viz = MeshCat.Visualizer()
    if draworigin
      setobject!( viz[:origin], Triad(originscale) )
    end

    # open a new browser tab if required
	open(viz, host=host, start_browser=start_browser)

    return viz
end

"""
   $(SIGNATURES)
High level interface to launch webserver process that draws a factor_graph_vis_type <: AbstractAmphitheatre, using Three.js and MeshCat.jl.
User factor_graph_vis_type should provide a visualize!(vis::Visualizer, factor_graph_vis_variables::T<:AbstractAmphitheatre ) function.
"""
function visualize(visdatasets::Vector{AbstractAmphitheatre};
					trans=Translation(0.0,0.0,0.0),
					quat::Rotations.Quat=Quat(1.0,0.0,0.0,0.0),
				   	host=ip"127.0.0.1",
					start_browser::Bool=true,
				    draworigin::Bool=true,
				    originscale::Float64=1.0)
    #
    global runAmphi

    runAmphi = true

    # the visualizer object itself
    vis = startMeshCatVisualizer(host=host,
								 start_browser=start_browser,
								 draworigin=draworigin,
								 originscale=originscale)
	setGlobalDrawTransform!(vis, trans=trans, quat=quat)

    # run the visualization loop #TODO add is task done to avoid multiple tasks getting lost
    runAphiTask = @async begin
		while runAmphi
			# iterate through all datasets #vir wat staan rose_fgl?
			for rose_fgl in visdatasets
				# each dataset should provide an visualize function
				visualize!(vis, rose_fgl)
			end
			# take a break and repeat
			sleep(1)
		end
		@info "visualize is finalizing."
	end

    return vis, runAphiTask
end

function restartVisualize(vis, visdatasets)
	global runAmphi
	runAmphi = true
	runAphiTask = @async begin
		while runAmphi
			# iterate through all datasets #vir wat staan rose_fgl?
			for rose_fgl in visdatasets
				# each dataset should provide an visualize function
				visualize!(vis, rose_fgl)
			end
			# take a break and repeat
			sleep(1)
		end
		@info "visualize is finalizing."
	end
	return runAphiTask
end

end


text:
module test_h3_lib_inspection

using Test
using H3.Lib

Lib.h3GetResolution
Lib.h3GetBaseCell
Lib.stringToH3
Lib.h3ToString
Lib.h3IsValid
Lib.h3IsResClassIII
Lib.h3IsPentagon

@test Lib.h3GetResolution(0x085283473fffffff) == 5

@test Lib.h3GetBaseCell(0x085283473fffffff) == 20

@test Lib.stringToH3("85283473fffffff") == 0x085283473fffffff

# https://medium.com/stuart-engineering/ruby-bindings-and-extensions-91c794eb9acd
bufSz = 17 # 16 hexadecimal characters plus a null terminator
buf = Base.unsafe_convert(Cstring, "")
Lib.h3ToString(0x085283473fffffff, buf, bufSz)
@test Base.unsafe_string(buf) == "85283473fffffff"

@test Lib.h3IsValid(0x85283473fffffff) == 1
@test Lib.h3IsValid(0x5004295803a8800) == 0

@test Lib.h3IsResClassIII(0x085283473fffffff) == 1

@test Lib.h3IsPentagon(0x085283473fffffff) == 0

end # module test_h3_lib_inspection


text:
using Requires

function __init__()
    @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" begin
        # FIXME: Should be restricted to dimensionless quantities because the
        # argument of a kernel function is in units of the interpolation grid
        # step size.
        convert_coordinate(T::Type{<:AbstractFloat}, c::Unitful.Quantity) =
            convert(T, c.val)
    end
end


text:
module JeszenszkiBasis

export
    AbstractSzbasis,
    Szbasis,
    RestrictedSzbasis,

    @sz_str,

    num_vectors,
    serial_num,
    site_max,
    sub_serial_num,
    to_str

include("basis.jl")
include("indexing.jl")
include("iteration.jl")
include("string.jl")
include("utilities.jl")

end