text:
using IntcodeVM, Combinatorics

# Let's go back to day 2 part 1 for some tests.

const day2_tests = ["1,9,10,3,2,3,11,0,99,30,40,50", "1,0,0,0,99", "2,3,0,3,99", "2,4,4,5,99,0", "1,1,1,4,99,5,6,0,99"];

for test in day2_tests
    println(IntcodeVM.run(parse.(Int,split(test,",")))[1]);
end

# And day 5

const day5_test = ["1002,4,3,4,33"];
const day5_test_result = IntcodeVM.run(parse.(Int,split(day5_test[1],",")))[1];
println("Day 5 Test: $(day5_test_result) ($(day5_test_result == [1002, 4, 3, 4, 99] ? "passed" : "failed"))");

const day5_part2_tests = ["3,9,8,9,10,9,4,9,99,-1,8", "3,9,7,9,10,9,4,9,99,-1,8",
    "3,3,1108,-1,8,3,4,3,99", "3,3,1107,-1,8,3,4,3,99"]
for test in day5_part2_tests
    (result,outputs) = IntcodeVM.run(parse.(Int,split(test,",")), inputs=[1],
            in=devnull, out=devnull);
    println("Day 5 Test: $(result) produces outputs $(outputs)")
end

# Day 7

const day7_examples = ["3,15,3,16,1002,16,10,16,1,16,15,15,4,15,99,0,0",
"3,23,3,24,1002,24,10,24,1002,23,-1,23,101,5,23,23,1,24,23,23,4,23,99,0,0",
"3,31,3,32,1002,32,10,32,1001,31,-2,31,1007,31,0,33,1002,33,7,33,1,33,31,31,1,32,31,31,4,31,99,0,0,0"]

function chainAmplifiers(code, sequence)
    inputSignal = 0
    for phaseSetting in sequence
        c = copy(code)
        inputSignal = popfirst!(IntcodeVM.run(c, inputs=[phaseSetting, inputSignal], in=devnull, out=devnull)[2])
    end
    return inputSignal
end

function bruteForce(code, phaseSettings::UnitRange{Int})
    maxSignal = (-1, undef)
    for order in permutations(phaseSettings)
        signal = chainAmplifiers(code, order)
        if (signal > maxSignal[1])
            maxSignal = (signal, order)
        end
    end
    return maxSignal
end

for i in 1:length(day7_examples)
    print("Day 7 Example $(i): ");
    println(bruteForce(parse.(Int,split(day7_examples[i],",")), 0:4))
end

# Day 9!

const day9_examples = ["109,1,204,-1,1001,100,1,100,1008,100,16,101,1006,101,0,99",
"1102,34915192,34915192,7,4,7,99,0",
"104,1125899906842624,99"];

for test in day9_examples
    print("Day 9 Test: ");
    println(IntcodeVM.run(parse.(Int,split(test,",")))[2]);
end

#=
--------------------;
109,1,              ; increase relative offset by 1
204,-1,             ; output vm[-1] using the relative offset
1001,100,1,100,     ; increment vm[100] by 1
1008,100,16,101,    ; write result of (vm[100] == 16) to vm[101]
1006,101,0,         ; if vm[101] == 0 goto 0
99                  ; exit
--------------------;
=#


text:
using LinearAlgebra

function cheb4c(N)

    # The function x,D4 = cheb4c(N) computes the fourth derivative
    # matrix D4 on Chebyshev interior points, incorporating the
    # clamped boundary conditions.
    # It is adapted from the function with the same name in dmsuite
    # developed for Matlab by JAC Weideman and SC Reddy.
    #
    # Adapted by M. Beneitez at KTH Mechanics.
    # beneitez@mech.kth.se
    #
    # Input:
    # N-2:  Order of the differentiation matrix.
    #
    # Output:
    # yF: Interior Chebyshev points
    # D4: Size [N-2,N-2] Fourth derivative matrix
    #

    eye = Matrix{Float64}(I, N-2, N-2);

    n1 = Int(floor(N/2-1)); n2  = Int(ceil(N/2-1));     # Indices used for flipping trick.

    k = collect(1:N-2);                        # Compute theta vector.
    th = k.*pi./(N-1)

    x = sin.(pi.*collect(N-3:-2:3-N)./(2*(N-1))); # Compute Chebyshev points.

    s = [sin.(th[1:n1]);reverse(sin.(th[1:n2]), dims = 1)]

    alpha = s.^4;
    beta1 = -4*s.^2 .*x./alpha;
    beta2 = 4*(3*x.^2 .-1)./alpha;
    beta3 = 24*x./alpha;
    beta4 = 24 ./alpha;
    B = [beta1'; beta2'; beta3'; beta4'];

    T = repeat(th/2,1,N-2);
    DX = 2*sin.(T'+T).*sin.(T'-T);          # Trigonometric identity.
    DX = [DX[1:n1,:]; -reverse(reverse(DX[1:n2,:], dims = 2), dims=1)];   # Flipping trick.
    DX[eye.==true] = ones(N-2,1);                      # Put 1's on the main diagonal of DX.

    ss = s.^2 .*(-1.0).^k;
    S  = repeat(ss, 1, N-2);
    C  = S./S';

    Z  = 1 ./DX;
    Z[eye.==true] = zeros(N-2,1);

    X  = Z';
    X  = X[X .!= 0.0];
    X  = reshape(X,N-3,N-2);

    Y  = ones(N-3,N-2);
    D  = Matrix{Float64}(I, N-2, N-2);
    DM = zeros(N-2,N-2,4);

    for ell = 1:4
        Y = cumsum([B[ell,:]'; ell.*Y[1:N-3,:].*X], dims = 1);
        D = ell.*Z.*(C.*repeat(diag(D), 1, N-2)-D);
        D[eye.==true] = Y[N-2,:];
        DM[:,:,ell] = D;
    end

    D4 = DM[:,:,4];

    return x,D4

end

ext: jl


text:
@testset "1604.alert-using-same-key-card-three-or-more-times-in-a-one-hour-period.jl" begin
    @test alert_names(
        ["daniel", "daniel", "daniel", "luis", "luis", "luis", "luis"],
        ["10:00", "10:40", "11:00", "09:00", "11:00", "13:00", "15:00"],
    ) == ["daniel"]
    @test alert_names(
        ["leslie", "leslie", "leslie", "clare", "clare", "clare", "clare"],
        ["13:00", "13:20", "14:00", "18:00", "18:51", "19:30", "19:49"],
    ) == ["clare", "leslie"]
    @test alert_names(["john", "john", "john"], ["23:58", "23:59", "00:01"]) == String[]
    @test alert_names(
        ["alice", "alice", "alice", "bob", "bob", "bob", "bob"],
        ["12:01", "12:00", "18:00", "21:00", "21:20", "21:30", "23:00"],
    ) == ["bob"]
end


text:
immutable Biweight <: ContinuousUnivariateDistribution
    μ::Float64
    σ::Float64

    Biweight(μ::Real, σ::Real) = (@check_args(Biweight, σ > zero(σ)); new(μ, σ))
    Biweight(μ::Real) = new(μ, 1.0)
    Biweight() = new(0.0, 1.0)
end

@distr_support Biweight d.μ - d.σ d.μ + d.σ

## Parameters
params(d::Biweight) = (d.μ, d.σ)

## Properties
mean(d::Biweight) = d.μ
median(d::Biweight) = d.μ
mode(d::Biweight) = d.μ

var(d::Biweight) = d.σ^2 / 7.0
skewness(d::Biweight) = 0.0
kurtosis(d::Biweight) = -2.9523809523809526  # = 1/21-3

## Functions
function pdf(d::Biweight, x::Float64)
    u = abs(x - d.μ) / d.σ
    u >= 1.0 ? 0.0 : 0.9375 * (1 - u^2)^2 / d.σ
end

function cdf(d::Biweight, x::Float64)
    u = (x - d.μ) / d.σ
    u <= -1.0 ? 0.0 :
    u >= 1.0 ? 1.0 :
    0.0625 * (u + 1.0)^3 * @horner(u,8.0,-9.0,3.0)
end

function ccdf(d::Biweight, x::Float64)
    u = (d.μ - x) / d.σ
    u <= -1.0 ? 1.0 :
    u >= 1.0 ? 0.0 :
    0.0625 * (u + 1.0)^3 * @horner(u,8.0,-9.0,3.0)
end

@quantile_newton Biweight

function mgf(d::Biweight, t::Float64)
    a = d.σ*t
    a2 = a^2
    a == 0 ? 1.0 :
    15.0 * exp(d.μ * t) * (-3.0 * cosh(a) + (a + 3.0/a) * sinh(a)) / (a2^2)
end

function cf(d::Biweight, t::Float64)
    a = d.σ * t
    a2 = a^2
    a == 0 ? 1.0+0.0im :
    -15.0 * cis(d.μ * t) * (3.0 * cos(a) + (a - 3.0/a) * sin(a)) / (a2^2)
end


text:
module MutID

using Compat
using DataFrames

function preprocess(csv_file, idx_mut)
    out_file= string(csv_file[1:end-4], "_new.tsv")
    out = open(out_file, "w")
    open(csv_file) do file
        while !eof(file)
            line = readline(file)
            #@show line
            row  = convert(Array{(@compat String),1}, split(line, ","))
            new_row = row[1:7], row[]
            #@show row
            mutid = row[idx_mut]
            mutids = convert(Array{(@compat String),1}, split(split(mutid,"\\")[2][4:end], ","))
            num_mutid = length(mutids)
            for i in 1:num_mutid
                out_row = vcat(row, mutids[i])
                @show length(out_row)
                write(out, join(out_row,"\t"), '\n')
            end
        end
    end
    close(out)
end

#preprocess("data/re.csv", 8)
#preprocess("data/stomach.csv", 6)
cosmic_path = "/haplox/data/cosmic/GRch37/v74/CosmicCompleteExport.tsv"
cosmic = readtable(cosmic_path, header=true) #x17
re = readtable("data/re.csv", header=false) # end
sto = readtable("data/stomach.csv", header=false) #end

function preprocess_dataframe(df, idx_mut)
    mutids = map(x->split(x,"\\")[2][4:end], df[idx_mut]) |>s-> map(x->split(x,","), s)
    ret_df = DataFrame()
    rowidx_df = 0
    for mutid in mutids
        rowidx_df += 1
        num_id = length(mutid)
        for i in 1:num_id
            ret_df = vcat(ret_df, df[rowidx_df,:])
        end
    end
    new_mutids = reduce(vcat, mutids) |> x->convert(Array{UTF8String,1}, x)

    hcat(ret_df, new_mutids)
end

re = preprocess_dataframe(re, 8)
sto = preprocess_dataframe(sto, 6)
#@show head(re,1)
#@show head(cosmic,1)
rename!(re, :x1_1, :Mutation_ID)
rename!(sto, :x1_1, :Mutation_ID)
ret_re = join(re, cosmic, on=:Mutation_ID, kind=:left)
ret_sto = join(sto, cosmic, on=:Mutation_ID, kind=:left)

writetable("ret_re.csv", ret_re)
writetable("ret_stomach.csv", ret_sto)

end


text:
#=
  Bayesian Data Analysis of Wine dataset

  I'm continue to testingMvNormal models to see how well it does on
  various datasets.

  From http://archive.ics.uci.edu/ml/datasets/Wine
  """
    Original Owners:

    Forina, M. et al, PARVUS -
    An Extendible Package for Data Exploration, Classification and Correlation.
    Institute of Pharmaceutical and Food Analysis and Technologies, Via Brigata Salerno,
    16147 Genoa, Italy.

    Donor:

    Stefan Aeberhard, email: stefan '@' coral.cs.jcu.edu.au

    Data Set Information:

    These data are the results of a chemical analysis of wines grown in the same region in Italy
    but derived from three different cultivars. The analysis determined the quantities of 13
    constituents found in each of the three types of wines.

    I think that the initial data set had around 30 variables, but for some reason I only have the
    13 dimensional version. I had a list of what the 30 or so variables were,
    but a.) I lost it, and b.), I would not know which 13 variables are
    included in the set.

    The attributes are (dontated by Riccardo Leardi, riclea '@' anchem.unige.it )
    1) Alcohol
    2) Malic acid
    3) Ash
    4) Alcalinity of ash
    5) Magnesium
    6) Total phenols
    7) Flavanoids
    8) Nonflavanoid phenols
    9) Proanthocyanins
    10)Color intensity
    11)Hue
    12)OD280/OD315 of diluted wines
    13)Proline

    In a classification context, this is a well posed problem with "well behaved" class structures.
    A good data set for first testing of a new classifier, but not very challenging.
  """

  The dataset is in the file wine.csv. The first entry is the class: 1, 2, 3 (Green,Red,Yellow)

  This version use different training and test sets.

  * training pct 0.75:
    Num correct: 40.0
    Num incorrect: 4.0
    Accuracy: 0.9090909090909091
    Confusion matrix:
    3×3 Matrix{Float64}:
    14.0   3.0  0.0
     0.0  17.0  1.0
     0.0   0.0  9.0
    num bad: 4

    Distribution of the categories (training dataset):
    2.0 : 53 (39.55223880597015%)
    3.0 : 39 (29.104477611940297%)
    1.0 : 42 (31.34328358208955%)

=#
using Turing, StatsPlots, Distributions
using DataFrames, StatsBase
using CSV
include("jl_utils.jl")

#
# Read dataset and preprocessing
#
df = CSV.read("wine.csv",DataFrame,header=false)

display(describe(df))
num_rows,num_cols = size(df)
println("Num rows (full dataset): $num_rows")

data = Matrix(df)
x = data[:,2:num_cols]
y = data[:,1]

# training_pct = 0.90
training_pct = 0.75
# training_pct = 0.66
# training_pct = 0.33
# training_pct = 0.11
# training_pct = 0.11
# training_pct = 3/num_rows

x_train, x_test, y_train, y_test = split_train_test(x,y,training_pct)

num_train_rows = length(y_train)
num_test_rows = length(y_test)

println("Num train rows: $num_train_rows num_test_rows: $num_test_rows")

categories = unique(y) |> sort
println("Mean value for each species: ", categories)
display(categories .|> s -> (s,mean(x[y.==s,:], dims=1)))
println("\n")

means_train = mean(x_train,dims=1)
stds_train = std(x_train,dims=1)
x_train = (x_train .- means_train) ./ stds_train

means_test = mean(x_test,dims=1)
stds_test = std(x_test,dims=1)
x_test = (x_test .- means_test) ./ stds_test

println("train: mean(x_train): ", mean(x_train,dims=1))
println("train: std(x_train): ", std(x_train,dims=1))

model = linear_regression_model(x_train,y_train)

# println("Prior:")
# chns = sample(model, Prior(), 10_000)
# display(chns)

println("\nModel:")
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, HMC(0.1,6), 1_000)
chns = sample(model, NUTS(), 1_000)

display(chns)

# display(plot(chns))
println("\nPredictions:")

mpred = linear_regression_model(x_test,missing)
pred = predict(mpred,chns)
display(pred)

cats = Dict(1=>1, 2=>2, 3=>3)
confusion, bads = confusion_matrix(y_test,pred,cats,false)
println("Confusion matrix:")
display(confusion)
println("num bad: $(length(bads))")

println("\nDistribution of the categories (training dataset):")
for (k,v) in make_hash(y_train)
  println("$k : $v ($(100*v/num_train_rows)%)")
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using SafeTestsets

@safetestset "Data Series                                                                     " begin include("dataseries_tests.jl") end
@safetestset "Time Series                                                                     " begin include("timeseries_tests.jl") end
@safetestset "Atomic Solutions                                                                " begin include("atomic_solutions_tests.jl") end
@safetestset "Abstract Solution                                                               " begin include("abstract_solution_tests.jl") end
@safetestset "Deterministic Solutions                                                         " begin include("deterministic_solutions_tests.jl") end
@safetestset "Stochastic Solutions                                                            " begin include("stochastic_solutions_tests.jl") end


text:
module ScientificTypes

# Dependencies
using Reexport
@reexport using ScientificTypesBase
export scitype, elscitype, scitype_union
using Tables
using CategoricalArrays
using ColorTypes
using PrettyTables
using Dates
import Distributions

# exports
export coerce, coerce!, autotype, schema, levels, levels!

# --------------------------------------------------------------------------------------
# Abbreviations

const ST  = ScientificTypesBase
const Arr = AbstractArray
const CArr = CategoricalArray
const Cat = CategoricalValue
const COLS_SPECIALIZATION_THRESHOLD = 30
const ROWS_SPECIALIZATION_THRESHOLD = 10000
const SCHEMA_SPECIALIZATION_THRESHOLD = Tables.SCHEMA_SPECIALIZATION_THRESHOLD

#---------------------------------------------------------------------------------------
# Define convention
struct DefaultConvention <: Convention end
const CONV = DefaultConvention()
# -------------------------------------------------------------
# vtrait function, returns either `Val{:table}()` or `Val{:other}()`
vtrait(X) = Val{ifelse(Tables.istable(X), :table, :other)}()

# -------------------------------------------------------------
# Includes

include("coerce.jl")
include("schema.jl")
include("autotype.jl")

include("convention/scitype.jl")
include("convention/coerce.jl")

end # module


text:
"""
This function calls rstan.
It can be used as follows:

    rstan(model_code, data, iter = 1000; pars = missing, nchains = 4, adapt_delta = 0.8, max_treedepth = 10)

    Feel free to edit this function and submit a pull request to add more functionality.
    The function can be edited by calling

    @edit rstan("", Dict())

    and a text editor should open up, pointing to this function.
    You may want to set your editor first, eg:
    ENV["EDITOR"] = "code"
"""
function rstan(model_code, data, iter = 1000; pars = missing, nchains = 4, adapt_delta = 0.8, max_treedepth = 10)
    @rput model_code
    @rput data
    @rput nchains
    @rput iter
    @rput pars
    R"res <- stan(model_code = model_code, data = data, iter = iter, pars = pars, chains = nchains, control = list(adapt_delta = $adapt_delta, max_treedepth = $max_treedepth))"
    R"samples <- rstan::extract(res)"
    println(R"res")
    @rget samples
end


text:
using RetroSignalModel
using Optim
using CSV

optim_params(targetratio=5)

dfparams = optim_params_threads(120; targetratio=5)

CSV.write("solution_rtgMTK_optim_5x.csv", dfparams)


text:
"""
    F = JopEnvelope(spc[, power=0.5, damping=0.0])

where `F` is the envelope operator with doman and range given by `spc::JetSpace`.
The Envelope is taken along the fastest dimension of the space.  For example, if
`spc=JetSpace{Float64,10,11}` and `A=rand(spc)`, then the envelope would be along each
column of `A`.

The envelope of `d` is computed as: `(d^2 + (Hd)^2 + damping)^(power/2)` where `Hd` is the Hilbert
transform of `d`.  The evelope is computed when `power=1`.  If `power=2`, then the square
of the envelope is computed, and so on.

# Notes
The passed in power is multiplied by 1/2, and is the power for simple envelope,
applied to the sum of the squares (d^2 + (Hd)^2). For example:

    power = 2.0 : (d^2 + (H d)^2)^1

If power < 2, and damping is not > 0, you may get Indian Bread (Nan) when envelope value is zero
Default damping factor is eps(T). If you know your traces are not zero, set damping=0 in constructor to avoid over damping.
"""
function JopEnvelope(spc::JetSpace{T,N}, power::Real=1.0; damping::Real=eps(T)) where {T,N}
    JopNl(dom = spc, rng = spc, f! = JopEnvelope_f!, df! = JopEnvelope_df!, df′! = JopEnvelope_df′!,
        upstate! = JopEnvelope_upstate!, s = (power=T(power/2), damping=T(damping), eu=zeros(spc), hu=zeros(spc)))
end
export JopEnvelope

function JopEnvelope_hu_and_eu!(m::AbstractArray, power, damping, hu, eu)
    hu .= imag.(hilbert(m))
    eu .= m.^2 .+ hu.^2 .+ damping
    #if (any(eu .== 0))
    #    @error "Zero eu value. Possibly envelope is zero and damping is also set to zero or something negative! \n
    #                Set damping to a small positive value"
    #end
end

JopEnvelope_upstate!(m::AbstractArray, s::NamedTuple) = JopEnvelope_hu_and_eu!(m, s.power, s.damping, s.hu, s.eu)

function JopEnvelope_f!(d::AbstractArray, m::AbstractArray; power, damping, hu, eu)
    JopEnvelope_hu_and_eu!(m, power, damping, hu, eu)
    d .= eu.^power
end

# (1/e) [diag(u) + diag(Hu) H] δ
# (1/e) [u ∘ δ + Hu ∘ H δ]
function JopEnvelope_df!(d::AbstractArray, m::AbstractArray; mₒ, power, hu, eu, kwargs...)
    d .= (power * 2) .* (mₒ .* m .+ hu .* imag.(hilbert(m))) .* (eu.^(power - 1))
end

# [diag(u) - H diag(Hu)] (1/e) δ
# (u ∘ δ)/e - H(Hu ∘ δ/e)
function JopEnvelope_df′!(m::AbstractArray, d::AbstractArray; mₒ, hu, eu, power, kwargs...)
    Ed = (power * 2) .* d .* (eu.^(power - 1))
    m .= mₒ .* Ed .- imag(hilbert(hu .* Ed))
end


text:
using TDD
using Test

@testset "TDD.jl" begin
    # Write your tests here.

    #1.

    undirected_graph = [
        [2, 3], #nodes reacheable from node 1
        [1],    #nodes reacheable from node 2
        [1],    #nodes reacheable from node 3
        [5],    #nodes reacheable from node 4
        [4]     #nodes reacheable from node 5
    ]

    A_undirected_graph = Bool[               #A short for "Adjacency Matrix Format"
                        1 1 1 0 0;
                        1 1 0 0 0;
                        1 0 1 0 0;
                        0 0 0 1 1;
                        0 0 0 1 1;
                        ]

    #based on Figure 1 from homework, node order 2 --> 1 --> 3; 4 --> 5
    directed_graph = [
        [3],    #nodes reacheable from node 1
        [1],    #nodes reacheable from node 2
        [],     #nodes reacheable from node 3
        [5],    #nodes reacheable from node 4
        []      #nodes reacheable from node 5
    ]

    A_directed_graph = Bool[
                        1 0 1 0 0;
                        1 1 0 0 0;
                        0 0 1 0 0;
                        0 0 0 1 1;
                        0 0 0 0 1;
                        ]

    #new test crazy case - a unidirectional train directed graph
    # 8 --> 7 --> 6 --> 5 --> 4 --> 3 --> 2 --> 1
    train_graph = [
        [],  #nodes reacheable from node 1
        [1], #nodes reacheable from node 2
        [2],
        [3],
        [4],
        [5],
        [6],
        [7]  #nodes reacheable from node 8
    ]

    A_train_graph = Bool[
                    1 0 0 0 0 0 0 0;
                    1 1 0 0 0 0 0 0;
                    0 1 1 0 0 0 0 0;
                    0 0 1 1 0 0 0 0;
                    0 0 0 1 1 0 0 0;
                    0 0 0 0 1 1 0 0;
                    0 0 0 0 0 1 1 0;
                    0 0 0 0 0 0 1 1;
                    ]

    @test DCneighbors(undirected_graph, 1) == Set([3, 2, 1]) #Sets are unordered
    @test DCneighbors(undirected_graph, 2) == Set([2, 1])
    @test DCneighbors(undirected_graph, 3) == Set([3, 1])
    @test DCneighbors(undirected_graph, 4) == Set([4, 5])
    @test DCneighbors(undirected_graph, 5) == Set([5, 4])

    @test DCneighbors(directed_graph, 1) == Set([1, 3])
    @test DCneighbors(directed_graph, 2) == Set([2, 1])
    @test DCneighbors(directed_graph, 3) == Set([3])
    @test DCneighbors(directed_graph, 4) == Set([4, 5])
    @test DCneighbors(directed_graph, 5) == Set([5])

    @test DCneighbors(train_graph, 1) == Set([1])
    @test DCneighbors(train_graph, 2) == Set([2,1])
    @test DCneighbors(train_graph, 3) == Set([3,2])
    @test DCneighbors(train_graph, 4) == Set([4,3])
    @test DCneighbors(train_graph, 5) == Set([5,4])
    @test DCneighbors(train_graph, 6) == Set([6,5])
    @test DCneighbors(train_graph, 7) == Set([7,6])
    @test DCneighbors(train_graph, 8) == Set([8,7])

    @test_throws ErrorException("node does not exist") DCneighbors(undirected_graph, 6)

    #2

    @test ICneighbors(undirected_graph, 1) == Set([1, 2, 3])
    @test ICneighbors(undirected_graph, 2) == Set([1, 2, 3])
    @test ICneighbors(undirected_graph, 3) == Set([1, 2, 3])
    @test ICneighbors(undirected_graph, 4) == Set([4, 5])
    @test ICneighbors(undirected_graph, 5) == Set([4, 5])

    @test ICneighbors(directed_graph, 1) == Set([1, 2, 3])
    @test ICneighbors(directed_graph, 2) == Set([1, 2, 3])
    @test ICneighbors(directed_graph, 3) == Set([1, 2, 3])
    @test ICneighbors(directed_graph, 4) == Set([4, 5])
    @test ICneighbors(directed_graph, 5) == Set([4, 5])

    @test_throws ErrorException("node does not exist") ICneighbors(undirected_graph, 6)

    #new tests
    @test ICneighbors(train_graph, 1) == Set([1, 2, 3,4,5,6,7,8])
    @test ICneighbors(train_graph, 2) == Set([1, 2, 3,4,5,6,7,8])
    @test ICneighbors(train_graph, 3) == Set([1, 2, 3,4,5,6,7,8])
    @test ICneighbors(train_graph, 4) == Set([1, 2, 3,4,5,6,7,8])
    @test ICneighbors(train_graph, 5) == Set([1, 2, 3,4,5,6,7,8])
    @test ICneighbors(train_graph, 6) == Set([1, 2, 3,4,5,6,7,8])
    @test ICneighbors(train_graph, 7) == Set([1, 2, 3,4,5,6,7,8])
    @test ICneighbors(train_graph, 8) == Set([1, 2, 3,4,5,6,7,8])

    #3

    @test AllComponents(directed_graph) == Set( [ Set([1, 2, 3]), Set([4, 5]) ] ) #for some reason need to wrap Set argument in []
    @test AllComponents(undirected_graph) == Set( [ Set([1, 2, 3]), Set([4, 5]) ] )
    @test AllComponents(train_graph) == Set( [ Set([1, 2, 3,4,5,6,7,8]) ] )

    #4 tests for adjacency matrix format -- see matrices above, tests below (copied, changed matrices input)

        @test DCneighbors(A_undirected_graph, 1) == Set([3, 2, 1])                              #Sets are unordered
        @test DCneighbors(A_undirected_graph, 2) == Set([2, 1]) #passes erreously
        @test DCneighbors(A_undirected_graph, 3) == Set([3, 1]) #passes erreously
        @test DCneighbors(A_undirected_graph, 4) == Set([4, 5])
        @test DCneighbors(A_undirected_graph, 5) == Set([5, 4])

        @test DCneighbors(A_directed_graph, 1) == Set([1, 3])
        @test DCneighbors(A_directed_graph, 2) == Set([2, 1]) #passes erreously
        @test DCneighbors(A_directed_graph, 3) == Set([3])
        @test DCneighbors(A_directed_graph, 4) == Set([4, 5])
        @test DCneighbors(A_directed_graph, 5) == Set([5])

        @test DCneighbors(A_train_graph, 1) == Set([1]) #passes erreously
        @test DCneighbors(A_train_graph, 2) == Set([2,1]) #passes erreously
        @test DCneighbors(A_train_graph, 3) == Set([3,2])
        @test DCneighbors(A_train_graph, 4) == Set([4,3])
        @test DCneighbors(A_train_graph, 5) == Set([5,4])
        @test DCneighbors(A_train_graph, 6) == Set([6,5])
        @test DCneighbors(A_train_graph, 7) == Set([7,6])
        @test DCneighbors(A_train_graph, 8) == Set([8,7])

        @test_throws ErrorException("node does not exist") DCneighbors(undirected_graph, 6) #passes erreously (though expected? this is technicallythetruth)

        #4-2

        @test ICneighbors(A_undirected_graph, 1) == Set([1, 2, 3])
        @test ICneighbors(A_undirected_graph, 2) == Set([1, 2, 3])
        @test ICneighbors(A_undirected_graph, 3) == Set([1, 2, 3])
        @test ICneighbors(A_undirected_graph, 4) == Set([4, 5])
        @test ICneighbors(A_undirected_graph, 5) == Set([4, 5])

        @test ICneighbors(A_directed_graph, 1) == Set([1, 2, 3])
        @test ICneighbors(A_directed_graph, 2) == Set([1, 2, 3])
        @test ICneighbors(A_directed_graph, 3) == Set([1, 2, 3])
        @test ICneighbors(A_directed_graph, 4) == Set([4, 5])
        @test ICneighbors(A_directed_graph, 5) == Set([4, 5])

        @test_throws ErrorException("node does not exist") ICneighbors(undirected_graph, 6) #passes erreously (though expected? this is technicallythetruth)

        #new tests
        @test ICneighbors(A_train_graph, 1) == Set([1, 2, 3,4,5,6,7,8])
        @test ICneighbors(A_train_graph, 2) == Set([1, 2, 3,4,5,6,7,8])
        @test ICneighbors(A_train_graph, 3) == Set([1, 2, 3,4,5,6,7,8])
        @test ICneighbors(A_train_graph, 4) == Set([1, 2, 3,4,5,6,7,8])
        @test ICneighbors(A_train_graph, 5) == Set([1, 2, 3,4,5,6,7,8])
        @test ICneighbors(A_train_graph, 6) == Set([1, 2, 3,4,5,6,7,8])
        @test ICneighbors(A_train_graph, 7) == Set([1, 2, 3,4,5,6,7,8])
        @test ICneighbors(A_train_graph, 8) == Set([1, 2, 3,4,5,6,7,8])

        #4-3

        @test AllComponents(A_directed_graph) == Set( [ Set([1, 2, 3]), Set([4, 5]) ] ) #for some reason need to wrap Set argument in []
        @test AllComponents(A_undirected_graph) == Set( [ Set([1, 2, 3]), Set([4, 5]) ] )
        @test AllComponents(A_train_graph) == Set( [ Set([1, 2, 3,4,5,6,7,8]) ] )

        #note at this point, 7 tests for Adjacency Matrix input pass, even though I have not implemented this...
        #need to go back figure out which tests erreously passed
        #I assume the error exceptions. But also which others?
        #Yep erroneous passes labelled in comments see above. Now gotta figure out why
        #relevant bugs are in DCneighbors
        #Ok so ICneighbors if input some Bool_matrix and a node, it does Bool_matrix[node] which just returns true/false, so makes sense all those tests fail right now
        #AllComponents depends on ICneighbors so also makes sense those tests fail
        #But DCneighbors does something different
        #If I run DCneighbors(A_undirected_graph, somenode) for example it returns some Set{Int64} which may or may not pass the test
        #Because it is just doing `Set(vcat(node, graph_input[node]))` which returns a Set{Int64}([node, 0 or 1 depending on if indexed value from Bool_matrix]) which may or may not be the right answer
        #not gonna fix DCneighbors right now so that it correctly fails the tests, afraid gonna mess something up
        #just keep this in mind. When I implement the functions with multiple dispatch specifying method for input Bool matrix this will end up being moot

end


text:
module DataSets

using ..GPs
using StaticArrays

"""
TrainingSet{NX, NY, T}

A data structure to store inputs (x) and outputs (y) for use in training a
Gaussian process. Wherever possible data is not copied but data structures are
reinterpreted for speed. However, this could lead to inadvertant data
modification by the user if sufficient care is not taken.
"""
struct TrainingSet{NX, NY, T}
    x::Vector{SVector{NX, T}}
    y::Vector{SVector{NY, T}}
    xm::Matrix{T}
    ym::Matrix{T}

    function TrainingSet(x::Vector{SVector{NX, T}}, y::Vector{SVector{NY, T}}) where {NX, NY, T <: Number}
        if length(x) != length(y)
            throw(ArgumentError("x and y must have the same length"))
        end
        new{NX, NY, T}(x,
                       y,
                       reinterpret(T, x, (NX, length(x))),
                       reinterpret(T, y, (NY, length(y))))
    end
    function TrainingSet(x::Matrix{T}, y::Matrix{T}) where {T <: Number}
        if size(x, 2) != size(y, 2)
            throw(ArgumentError("x and y must have the same number of columns"))
        end
        new{size(x, 1), size(y, 1), T}(reinterpret(SVector{size(x, 1), T}, x, (size(x, 2),)),
                                       reinterpret(SVector{size(y, 1), T}, y, (size(y, 2),)),
                                       x,
                                       y)
    end
end

TrainingSet(x, y) = TrainingSet(toSVector(x), toSVector(y))

export TrainingSet

samplecount(data::TrainingSet) = length(data.x)

export samplecount

end


text:
export Predicate
export all
export exists

"""
    Represents a predicate.

A proposition-valued function of one or more variables.

fields: value and m.
"""
struct Predicate <: Logic
    value::Logic
    m::Number
    Predicate(p::Logic, m::Number) = new(p, m)
end

"""
    exists(x, P)

Construct a proposition with the given Predicate `P` and a single variable `x`.
"""
exists(x, P::Predicate) = !all(x, !P)

"""
    all(x, P)

Construct a proposition with the given Predicate `P` and a single variable `x`.
"""
all(x, P::Predicate) = value(P)

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
export RCM

#------------------------------------------------------------------------------
# definition & construction

struct RCM{A<:Real,𝒻₁,𝒻₂,𝒻₃,𝒻₄,C<:AbstractNumericalCore}
    Pₑ::Vector{A}
    Tₑ::Vector{A} #fossil initial condition
    Pᵣ::Vector{A}
    P::Vector{A}
    T::Vector{A}
    g::A #gravity [m/s^2]
    𝒻μ::𝒻₁ #mean molar mass function (T,P)->μ [kg/mole]
    𝒻S::𝒻₂ #stellar insolation at TOA ν->I [W/m^2/cm^-1]
    𝒻a::𝒻₃ #surface albedo ν->a [unitless]
    𝒻cₚ::𝒻₄ #heat capacity function (T,P)->cₚ [J/kg/K]
    cₛ::A #surface heat capacity [J/m^2/K]
    𝒜::AcceleratedAbsorber
    ν::Vector{Float64}
    core::C
    F::FluxPack{A}
    np::Int64 #number of atmospheric cells (including surface)
    nν::Int64
    nᵤ::Vector{Int64} #update counters for rows, just curious
    R::Vector{A} #net radiative fluxes at cell edges
    H::Vector{A} #cell heating rate [K/s]
    J::Matrix{A} #Jacobian
end

function Base.show(io::IO, ℛ::RCM{A}) where {A}
    @unpack Pₑ, ν, T, core, g, np, nν = ℛ
    print(io, "RCM{$A,$(typeof(core))}\n")
    print(io, "  gravity is $g m/s^2\n")
    Pmin, Pmax = round(minimum(Pₑ), sigdigits=4), round(maximum(Pₑ), sigdigits=4)
    print(io, "  $np pressure layers ∈ [$Pmin, $Pmax] Pa\n")
    νmin, νmax = round(minimum(ν), sigdigits=4), round(maximum(ν), sigdigits=4)
    print(io, "  $nν wavenumber coordinates ∈ [$νmin, $νmax] cm^-1\n")
    Tmin, Tmax = round(minimum(T), sigdigits=4), round(maximum(T), sigdigits=4)
    print(io, "  temperature currently ∈ [$Tmin, $Tmax] K")
end

function RCM(Pₑ::AbstractVector{<:Real},
             Tₑ::AbstractVector{<:Real},
             g::Real,
             𝒻μ,
             𝒻S,
             𝒻a,
             𝒻cₚ,
             cₛ::Real,
             absorbers...;
             core::AbstractNumericalCore=Discretized(),
             radmul::Int=2)
    #demand pressure coordiantes are ascending
    W = typeof(float(Tₑ[1]))
    idx = sortperm(Pₑ)
    Pₑ = collect(W, Pₑ[idx])
    Tₑ = collect(W, Tₑ[idx])
    #check the numbers
    np = length(Pₑ)
    @assert length(Tₑ) == np "must have same number of initial temperature and pressure values"
    #values at midpoints between edges (cell centers) and surface
    P = similar(Pₑ)
    T = similar(Tₑ)
    for i ∈ 1:np-1
        P[i] = (Pₑ[i] + Pₑ[i+1])/2
        T[i] = (Tₑ[i] + Tₑ[i+1])/2
    end
    P[end] = Pₑ[end]
    T[end] = Tₑ[end]
    #build up extra radiative nodes with weighted averaging
    @assert ((radmul % 2 == 0) | (radmul == 1)) "radmul must be an even integer or 1"
    nrad = radmul*(np - 1) + 1
    Pᵣ = similar(Pₑ, nrad)
    P₁ = @view Pₑ[1:end-1]
    Pᵣ[1:np-1] .= P₁
    P₂ = @view Pₑ[2:end]
    i = np
    for j ∈ 2:radmul
        w₁ = j - 1
        w₂ = radmul - w₁
        @. Pᵣ[i:i+np-2] = (w₁*P₁ + w₂*P₂)/radmul
        i += np - 1
    end
    Pᵣ[end] = Pₑ[end]
    sort!(Pᵣ)
    #parse absorbers
    𝒰, _, nν = unifyabsorbers(absorbers)
    #prepare accelerated absorption calculations
    𝒜 = AcceleratedAbsorber(Tₑ, Pₑ, 𝒰)
    #handle types
    cₛ = convert(W, cₛ)
    g = convert(W, g)
    #make room for radiative calculations
    F = FluxPack(nrad, nν, W)
    #cross-section update tracker
    nᵤ = zeros(Int64, np)
    #make room for heating calculations
    R = zeros(W, np)
    H = zeros(W, np)
    J = zeros(W, np, np)
    #construct
    RCM(Pₑ, Tₑ, Pᵣ, P, T, g, 𝒻μ, 𝒻S, 𝒻a, 𝒻cₚ, cₛ, 𝒜, 𝒜.ν, core, F, np, nν, nᵤ, R, H, J)
end

#------------------------------------------------------------------------------
# time derivatives and such

export heating!
function heating!(ℛ::RCM, δᵤ::Real=1)::Nothing
    @unpack Pₑ, Tₑ, Pᵣ, P, T, g, 𝒻μ, 𝒻S, 𝒻a, 𝒻cₚ, cₛ, 𝒜, core, F, np, nν, nᵤ, R, H = ℛ

    𝒻T = AtmosphericProfile(P, T)
    radiate!(F, core, Pᵣ, g, 𝒻T, 𝒻μ, 𝒻S, 𝒻a, 𝒜)

    #====================================================================
    Coordinates are a little tricky here. The "upward" flux Fup is going
    out of the atmosphere, away from the surface. However, the positive
    direction is the other way because pressure is increasing toward the
    surface and the model discretization is along pressure coordinates.
    So, although the net flux is typically F⁺ - F⁻, it must be reversed
    here to respect the true positive direction.
    ====================================================================#
    𝒻Fnet = AtmosphericProfile(Pᵣ, F.Fnet)
    @. R = -𝒻Fnet(Pₑ) #F.F⁻ - F.F⁺

    #==========================
    heating rates
    ==========================#
    for i ∈ 1:np-1
        #heat capacity [J/kg/K]
        cₚ = 𝒻cₚ(T[i],P[i])
        #pressure layer thickness [Pa]
        ΔP = Pₑ[i+1] - Pₑ[i]
        #radiation out of the cell minus radiation into it
        ΔR = R[i] - R[i+1]
        #cell heating rate [K/s]
        H[i] = (g/cₚ)*ΔR/ΔP
    end
    #surface heating
    H[end] = R[end]/cₛ
    #H[end] = (F.F⁺[1] - F.F⁻[end])/cₛ

    return nothing
end

export step!
function step!(ℛ::RCM, Δt::Real)::Nothing
    heating!(ℛ)
    ℛ.T .+= Δt*ℛ.H
    return nothing
end

export jacobian!
function jacobian!(ℛ::RCM, ϵ::Real=1)::Nothing
    heating!(ℛ)
    H = ℛ.H[:]
    T = ℛ.T
    for i ∈ 1:ℛ.np
        #perturb
        T[i] += ϵ
        #evaluate
        heating!(ℛ)
        #finite diff
        for j ∈ 1:length(H)
            ℛ.J[j,i] = (ℛ.H[j] - H[j])/ϵ
        end
        #restore
        T[i] -= ϵ
    end
    return nothing
end


text:
# # Prior model implementation in BAT.jl
#
# In this notebook, the implementation of the prior model is demonstrated,
# starting with a simple two-component model and scaling up to the full 9 component model.

using Distributions, StatsBase, LinearAlgebra
using Plots, SpecialFunctions, Printf, Random, ValueShapes
using BAT, DensityInterface
const sf = SpecialFunctions;

# ## Simple two-component model: gluons
#
# First we start with a simpler two-component model, and show all the
# steps explicity for clarity.
#
# ### Forward model
#
# The gluon distributions are parametrised by
# ```math
# x g(x) = A_{g1} x^{\lambda_{g1}}(1-x)^5 + A_{g2} x^{\lambda_g2}(1-x)^{K_g}.
# ```
#
# We also want to impose
#
# ```math
# \int_0^1 x g(x) dx = A_{g1} B(\lambda_{g1}+1, 1+5) + A_{g2} B(\lambda_{g2}+1, K_g+1)= 1,
# ```
#
# with B(.,.) the Beta function.
#
# Start by defining some useful functions:

function xg1x(x, λ_g1, θ_2)
    A_g1 = θ_2 / sf.beta(λ_g1 + 1, 5 + 1)
    return A_g1 * x^λ_g1 * (1 - x)^5
end

function xg2x(x, λ_g2, K_g, θ_1)
    A_g2 = θ_1 / sf.beta(λ_g2 + 1, K_g + 1)
    return A_g2 * x^λ_g2 * (1 - x)^K_g
end

function xgx(x, λ_g1, λ_g2, K_g, θ)
    xg1 = xg1x(x, λ_g1, θ[2])
    xg2 = xg2x(x, λ_g2, K_g, θ[1])
    return xg1 + xg2
end

# Choose true values for the high-level parameters and show what the resulting model looks like.

θ = [0.5, 0.5]
λ_g1 = -0.7 # rand(Uniform(-1, 0))
λ_g2 = 0.5 # rand(Uniform(0, 1))
K_g = 3 # rand(Uniform(2, 10))
truths = (θ = θ, λ_g1 = λ_g1, λ_g2 = λ_g2, K_g = K_g);

A_g2 = θ[1] / sf.beta(λ_g2+1, K_g+1)
A_g1 = θ[2] / sf.beta(λ_g1+1, 5+1);

# Check integral = 1
total = A_g1 * sf.beta(λ_g1+1, 5+1) + A_g2 * sf.beta(λ_g2+1, K_g+1)
print("Integral = ", total)

# Plot true model
x_grid = range(0, stop=1, length=50)

xg2 = A_g2 * x_grid.^λ_g2 .* (1 .- x_grid).^K_g
xg1 = A_g1 * x_grid.^λ_g1 .* (1 .- x_grid).^5

plot(x_grid, [xg1x(x, λ_g1, θ[2]) for x in x_grid],
    alpha=0.7, label="x g1(x)", lw=3, color="blue")
plot!(x_grid, [xg2x(x, λ_g2, K_g, θ[1]) for x in x_grid],
    alpha=0.7, label="x g2(x)", lw=3, color="green")
plot!(x_grid, [xgx(x, λ_g1, λ_g2, K_g, θ) for x in x_grid],
    alpha=0.7, label="x g1(x) + x g2(x)", lw=3, color="red")
plot!(xlabel="x")

# Now, for the purposes of testing the prior implementation,
# sample some data from this distribution assuming that the data
# are produced by integrating over the function in different bins,
# and multiplying by some factor.
# Then, plot the model and data to compare.

bins = 0.0:0.05:1.0
bin_widths = bins[2:end] - bins[1:end-1]
bin_centers = (bins[1:end-1] + bins[2:end]) / 2

N = 1000
nbins = size(bin_centers)[1]

expected_counts = zeros(nbins)
observed_counts = zeros(Integer, nbins)
for i in 1:nbins
    xg = xgx(bin_centers[i], λ_g1, λ_g2, K_g, θ) * N
    expected_counts[i] = bin_widths[i] * xg
    observed_counts[i] = rand(Poisson(expected_counts[i]))
end

plot(bin_centers, [xgx(x, λ_g1, λ_g2, K_g, θ) for x in bin_centers] .* bin_widths * N,
    alpha=0.7, label="Expected", lw=3, color="red")
scatter!(bin_centers, observed_counts, lw=3, label="Observed", color="black")

# Store the data in a simple dict to pass to the likelihood later.

data = Dict()
data["N"] = N
data["bin_centers"] = bin_centers;
data["observed_counts"] = observed_counts;
data["bin_widths"] = bin_widths;

# ### Fit
#
# To fit this example data, we choose a prior over our hyperparameters `θ`, `λ_g1`, `λ_g2` and `K_g`.
#
# We decide to choose a sensible Dirichlet prior, and have a look at
# some samples to help understand what this means.

dirichlet = Dirichlet([1, 1])
test = rand(dirichlet, 1000)
plot(append!(Histogram(0:0.1:1), test[1, :]))
plot!(append!(Histogram(0:0.1:1), test[2, :]))

# So, we see that this means all weights are equally likely for
# both components.
#
# Now we can define the prior:
#
# ```math
# \mathrm{Prior} - P(\theta) P(\lambda_{g1}) P(\lambda_{g2}) P(K_g).
# ```

prior = NamedTupleDist(
    θ = Dirichlet([1, 1]),
    λ_g1 = Uniform(-1, 0), #Truncated(Normal(-0.7, 0.1), -1, 0),
    λ_g2 = Uniform(0, 1), #Truncated(Normal(0.7, 0.1), 0, 1)
    K_g =  Uniform(2, 10), # Truncated(Normal(1, 2), 0, 10),
);

# Define a simple poisson likelihood that describes the test data that we generated above.

likelihood = let d = data, f = xgx

    observed_counts = d["observed_counts"]
    bin_centers = d["bin_centers"]
    bin_widths = d["bin_widths"]
    N = data["N"]

    logfuncdensity(function (params)
        function bin_log_likelihood(i)
            xg = f(
                bin_centers[i], params.λ_g1, params.λ_g2, params.K_g, params.θ
            )
            expected_counts = bin_widths[i] * xg * N
            logpdf(Poisson(expected_counts), observed_counts[i])
        end

        idxs = eachindex(observed_counts)
        ll_value = bin_log_likelihood(idxs[1])
        for i in idxs[2:end]
            ll_value += bin_log_likelihood(i)
        end

        return ll_value
    end)
end

# The prior and likelihood can be passed to `BAT.jl` via the
# `PosteriorDensity`. We can then sample this density using `bat_sample()`.

posterior = PosteriorDensity(likelihood, prior);
samples = bat_sample(
    posterior,
    MCMCSampling(mcalg=MetropolisHastings(), nsteps=10^4, nchains=4)
).result;

# The `SampledDensity` gives a quick overview of the results.

SampledDensity(posterior, samples)

#  ### Visualise results
# We can (roughly) check how well the fit reconstructs the truth with a simple comparison.

x_grid = range(0, stop=1, length=50)

plot()
for i in 1:200
    random_ind = rand(1:7000)
    s = samples[random_ind].v
    xg = [xgx(x, s.λ_g1, s.λ_g2, s.K_g, s.θ) for x in bin_centers]
    plot!(bin_centers, xg .* bin_widths * N, alpha=0.1, lw=3,
        color="darkorange", label="",)
end

xg = [xgx(x, λ_g1, λ_g2, K_g, θ) for x in bin_centers]
plot!(bin_centers, xg .* bin_widths * N, alpha=0.7,
    label="Expected", lw=3, color="red")

scatter!(bin_centers, observed_counts, lw=3, label="Observed", color="black")

plot!(xlabel="x")

# In the above plot, the red line represents the truth, and the set
# of fainter lines represent samples from the posterior.
#
# We can also look at marginal distributions for different parameters...
plot(
    samples, :(λ_g2),
    mean = true, std = true,
    nbins = 50,
)

# ## Full model with all components
#
# We can now extend this approach to the full 9 components,
# there is not so much documentation here, as it follows the above case.
#
#
# ### Forward model
# Here the parametrisation of each component in decreasing order of importance.
# We also use the helper functions provided by `PartonDensity` to keep things tidy.

using PartonDensity

#

pdf_params = PDFParameters(λ_u=0.7, K_u=4.0, λ_d=0.5, K_d=6.0,
    λ_g1=-0.4, λ_g2=0.7, K_g=6.0, λ_q=0.5, weights=[1, 0.5, 0.3, 0.2, 0.1, 0.1, 0.1])

# Sanity check

int_xtotx(pdf_params)

# Plot true model

plot_input_pdfs(pdf_params)

# Generate example data

bins = 0.0:0.05:1.0
bin_widths = bins[2:end] - bins[1:end-1]
bin_centers = (bins[1:end-1] + bins[2:end]) / 2

N = 1000
nbins = size(bin_centers)[1]

expected_counts = zeros(nbins)
observed_counts = zeros(Integer, nbins)
for i in 1:nbins
    xt = xtotx(bin_centers[i], pdf_params) * N
    expected_counts[i] = bin_widths[i] * xt
    observed_counts[i] = rand(Poisson(expected_counts[i]))
end

# Plot data and expectation

plot(bin_centers, [xtotx(x, pdf_params) for x in bin_centers] .* bin_widths * N,
    alpha=0.7, label="Expected", lw=3, color="red")
scatter!(bin_centers, observed_counts, lw=3, label="Observed", color="black")

# Store the data

data = Dict()
data["N"] = N
data["bin_centers"] = bin_centers;
data["observed_counts"] = observed_counts;
data["bin_widths"] = bin_widths;

# ### Fit
#
# Prior

prior = NamedTupleDist(
    θ = Dirichlet([1, 1, 1, 1, 1, 1, 1]),
    λ_u = Truncated(Normal(pdf_params.λ_u, 0.1), 0, 1), #  Uniform(0, 1),
    K_u = Uniform(2, 6),
    λ_d = Truncated(Normal(pdf_params.λ_d, 0.1), 0, 1), # Uniform(0, 1),
    K_d = Uniform(2, 8),
    λ_g1 = Uniform(-1, 0),
    λ_g2 = Uniform(0, 1),
    K_g =  Uniform(2, 8),
    λ_q = Truncated(Normal(pdf_params.λ_q, 0.1), -1, 0),
);

# Likelihood

likelihood = let d = data, f = xtotx

    observed_counts = d["observed_counts"]
    bin_centers = d["bin_centers"]
    bin_widths = d["bin_widths"]
    N = data["N"]

    logfuncdensity(function (params)
        function bin_log_likelihood(i)
            xt = f(bin_centers[i], params.λ_u, params.K_u, params.λ_d, params.K_d,
                    params.λ_g1, params.λ_g2, params.K_g, params.λ_q, Vector(params.θ))
            expected_counts = bin_widths[i] * xt * N
            logpdf(Poisson(expected_counts), observed_counts[i])
        end

        idxs = eachindex(observed_counts)
        ll_value = bin_log_likelihood(idxs[1])
        for i in idxs[2:end]
            ll_value += bin_log_likelihood(i)
        end

        return ll_value
    end)
end

# Run fit

posterior = PosteriorDensity(likelihood, prior);
samples = bat_sample(
    posterior,
    MCMCSampling(mcalg=MetropolisHastings(), nsteps=10^4, nchains=4)
).result;

# ### Visualise results

x_grid = range(0, stop=1, length=50)

plot()
for i in 1:200
    random_ind = rand(1:7000)
    s = samples[random_ind].v
    xt = [xtotx(x, s.λ_u, s.K_u, s.λ_d, s.K_d,
            s.λ_g1, s.λ_g2, s.K_g, s.λ_q, Vector(s.θ)) for x in bin_centers]
    plot!(bin_centers, xt .* bin_widths * N, alpha=0.1, lw=3,
        color="darkorange", label="")
end

xt = [xtotx(x, pdf_params) for x in bin_centers]
plot!(bin_centers, xt .* bin_widths * N, alpha=0.7,
    label="Expected", lw=3, color="red")

scatter!(bin_centers, observed_counts, lw=3, label="Observed", color="black")

plot!(xlabel="x")

# These first results are promising. We can also try changing the input
# parameters and priors to explore the performance of the fit.


text:
__precompile__()

module celerite

include("gp.jl")
include("sturms_theorem.jl")

end


text:
using MatrixCalculus
model = Dynamics.cartpole_uncertain
n = model.n; m = model.m; r = model.r

# Problem
x0 = [0.;0.;0.;0.]
xf = [0.;pi;0.;0.]

E0 = Diagonal(1.0e-6*ones(n))
H0 = zeros(n,r)
D = Diagonal([4])

N = 101 # knot points

Q = Diagonal(ones(n))
R = Diagonal(0.1*ones(m))
Rh = 0.
Qf = Diagonal(0.0*ones(n))

Q_lqr = Diagonal([10.;10.;1.;1.])
R_lqr = Diagonal(ones(m))
Qf_lqr = Diagonal([100.;100.;100.;100.])

Q_r = Q_lqr
R_r = R_lqr
Qf_r = Qf_lqr

tf0 = 5.
dt = tf0/(N-1)

h_max = dt
h_min = dt

u_max = 10.
u_min = -10.

eig_thr = 1.0e-3

# problem
obj = LQRObjective(Q,R,Qf,xf,N)

goal_con = goal_constraint(xf)
bnd_con = BoundConstraint(n,m,u_min=u_min,u_max=u_max,trim=true)
con = Constraints([bnd_con],N)

prob = TrajectoryOptimization.Problem(model, obj,constraints=con, N=N, tf=tf0, x0=x0, dt=dt)
prob.constraints[N] += goal_con # add goal constraint
copyto!(prob.X,line_trajectory(x0,xf,N)) # initialize state with linear interpolation
copyto!(prob.U,[0.01*rand(m) for k = 1:N-1]) # initialize control with rand

solver = DIRTRELSolver(E0,H0,D,Q,R,Rh,Qf,xf,Q_lqr,R_lqr,Qf_lqr,Q_r,R_r,Qf_r,eig_thr,h_max,h_min)

solve!(prob,solver)

plot(prob.X)
plot(prob.U)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
@testset "The SO(3)" begin
  import LinearAlgebra: det
  import Base: zeros
  import Random
  x = SOPoint([1. 0. 0.; 0. 1. 0.; 0. 0. 1.])
  xanti = SOPoint([1. 0. 0.; 0. -1. 0.; 0. 0. -1.])
  y = SOPoint([1 0 0; 0 cos(pi/5) -sin(pi/5); 0 sin(pi/5) cos(pi/5)])
  z = SOPoint([1 0 0; 0 cos(pi/3) -sin(pi/3); 0 sin(pi/3) cos(pi/3)])
  M = Rotations(3)
  ω = SOTVector(0.2*[0 1 2; -1 0 3; -2 -3 0])
  η = randomTVector(M,x)
  ξ = log(M, x ,y)
  μ = log(M, x, z)
  y2 = exp(M, x, ξ)
  z2 = exp(M, x, μ)
  x2 = addNoise(M,x,0.5)
  x3 = addNoise(M,x,500)
  w = randomMPoint(M)
  s = rand(Float64)
  Random.seed!(1)

  # Test unary operator
  # Test Dimension
  @test manifoldDimension(x)==3
  @test manifoldDimension(M)==3
  # Test Minus operator
  @test getValue( (-η) ) == getValue(-η)
  # Test Cauchy-Schwarz-Inequation
  @test dot(M,x,ω,η)^2 <= norm(M,x,ω)^2*norm(M,x,η)^2
  @test dot(M,x,ω,η) <= norm(M,x,ω)*norm(M,x,η)
  # Test dot
  @test dot(M,x,ω,η) ≈ dot(M,x,η,ω) atol = 10.0^(-16)
  @test dot(M,x,ω,η+ξ) ≈ dot(M,x,ω,η) + dot(M,x,ω,ξ) atol = 10.0^(-15)
  @test dot(M,x,ω,s*η) ≈ s*dot(M,x,ω,η) atol = 10.0^(-14)
  # Test norm
  @test norm(M,x,ω)^2 ≈ dot(M,x,ω,ω) atol = 10.0^(-14)
  @test norm(M,x,ω) ≈ sqrt(dot(M,x,ω,ω)) atol = 10.0^(-16)
  # Test distance
  @test distance(M,x,x) ≈ 0 atol = 10.0^(-16)
  @test distance(M,x,xanti) ≈ π atol=10.0^(-16) # antipodal ponts
  @test distance(M,x,z) ≈ norm(M,x,μ) atol=10.0^(-16)
  # Test Exponential and logarithm for both usual and antipodal points
  @test norm(getValue(log(M, x, y2)) - getValue(ξ)) ≈ 0 atol = 10.0^(-16)
  @test norm(getValue(log(M, x, z2)) - getValue(μ)) ≈ 0 atol = 10.0^(-16)
  @test norm(getValue(exp(M,x,log(M,x,y))) - getValue(y)) ≈ 0 atol = 10.0^(-15)
  @test norm(getValue(log(M,x,exp(M,x,ω))) - getValue(ω)) ≈ 0 atol = 10.0^(-15)
  @test_throws ErrorException log(M,x,xanti)
  # Test randomMPoint and randomTVector
  N = Rotations(1)
  @test norm(getValue(randomMPoint(N)) - ones(1,1)) ≈ 0 atol = 10.0^(-16)
  @test det(getValue(w)) ≈ 1 atol = 10.0^(-9)
  @test norm(transpose(getValue(w))*getValue(w) - one(getValue(w))) ≈ 0 atol = 10.0^(-14)
  @test norm(getValue(η) + transpose(getValue(η))) ≈ 0 atol = 10.0^(-16)
  @test getValue(randomTVector(Rotations(1), SOPoint(ones(1,1)))) == zeros(1,1)
  # Test addNoise
  @test det(getValue(x2)) ≈ 1 atol = 10.0^(-9)
  @test norm(transpose(getValue(x2))*getValue(x2) - one(getValue(x))) ≈ 0 atol=10.0^(-15)
  @test det(getValue(x3)) ≈ 1 atol = 10.0^(-9)
  @test norm(transpose(getValue(x3))*getValue(x3) - one(getValue(x))) ≈ 0 atol=10.0^(-12)
  #Test retraction
  @test norm(getValue(inverseRetractionQR(M,x,retractionQR(M,x,ω))) - getValue(ω)) ≈ 0 atol = 10.0^(-14)
  @test norm(getValue(inverseRetraction(M,x,retraction(M,x,ω))) - getValue(ω)) ≈ 0 atol = 10.0^(-14)
  @test norm(getValue(inverseRetractionPolar(M,x,retractionPolar(M,x,ω))) - getValue(ω)) ≈ 0 atol = 10.0^(-14)
  @test norm(transpose(getValue(retractionQR(M,x,ω))) * getValue(retractionQR(M,x,ω)) - one(getValue(x))) ≈ 0 atol = 10.0^(-15)
  @test norm(transpose(getValue(retraction(M,x,ω))) * getValue(retraction(M,x,ω)) - one(getValue(x))) ≈ 0 atol = 10.0^(-15)
  @test norm(transpose(getValue(retractionPolar(M,x,ω))) * getValue(retractionPolar(M,x,ω)) - one(getValue(x))) ≈ 0 atol = 10.0^(-15)
  #Test parallelTransport
  @test norm(getValue(parallelTransport(M,x,y,ξ)) - getValue(ξ)) ≈ 0 atol = 10.0^(-16)
  #Test zeroTVector
  @test norm(getValue(zeroTVector(M,x))) ≈ 0 atol = 10.0^(-16)
  #Test Lie Group capabilities
  @test distance(M, x⊗y, SOPoint( transpose(getValue(x))*getValue(y)) ) ≈ 0 atol = 10.0^(-16)
  #Test validateMPoint and validateTVector
  @test validateMPoint(M,w) == true
  @test validateTVector(M,x,η) == true
  ynot1 = SOPoint(0.5*[1 0 0; 0 cos(pi/5) -sin(pi/5); 0 sin(pi/5) cos(pi/5)])
  ynot2 = SOPoint([1 2 3; 0 2 1; 0 0 0.5])
  @test_throws ErrorException validateMPoint(M,ynot1)
  @test_throws ErrorException validateMPoint(M,ynot2)
  ωnot = SOTVector(0.2*[1 1 2; -1 1 3; -2 -3 1])
  @test_throws ErrorException validateTVector(M,x,ωnot)
end


text:
include("../src/nn.jl")
using Test
Random.seed!(1)

"""
charge distribution -> potential (Poisson's eqn), electric field (Gauss's law)
"""

# input Scalar field
inranks = [0]
# output scalar field, vector field
outranks = [0, 1]
sz=(8,8,8)
dx = 0.1
dV=dx^3
rmax = 2dx
lmax = 1

# charge distribution
x = [zeros(sz...)]
ix=[5, 5, 5]
x[1][ix...] = 1.0
X=[x]

# generate data
# Green's fn for Poisson, Gauss
f1 = LinearOperator(:potential,dx;rmax=rmax)
f2 = LinearOperator(:field,dx;rmax=rmax)

y1 = f1(X[1])
y2 = f2(X[1])

# check
v=[1,1,1]
ix = ix.+v
r=norm(v)*dx
@test y1[1][ix...]≈dV/r
@test [y2[i][ix...] for i = 1:3]≈dV/r^2*ones(3)/sqrt(3)

##
# train
# linear layer: tensor field convolution
L = EquivConv(inranks, outranks, dx; rmax = rmax)

function nn(X)
    L(X)
end

function loss()
    y1hat, y2hat = nn(X)
    l1 = Flux.mae(toArray(y1), toArray(y1hat))
    l2 = Flux.mae(toArray(y2), toArray(y2hat))
    l = l1 + l2
    println(l)
    l
end
loss()
##
ps = Flux.params(L)
data = [()]
opt = ADAM(0.1)

println("===\nTraining")
for i = 1:5
    # global doplot = i % 50 == 0
    Flux.train!(loss, ps, data, opt)
end

##
Random.seed!(1)
n=4
inranks=[0,0]
outranks=[0]
X=[[rand(n,n,n)],[rand(n,n,n)]]
y=[[X[1][1].*X[2][1]]]

# train
# linear layer: tensor field convolution

A = EquivAttn(inranks, outranks)

function nn(X)
    A(X)
end

function loss()
    yhat = nn(X)
    l = Flux.mae(toArray(y[1]), toArray(yhat[1]))
    println(l)
    l
end
loss()

ps = Flux.params(A)
data = [()]
opt = ADAM(0.1)

println("===\nTraining")
for i = 1:10
    # global doplot = i % 50 == 0
    Flux.train!(loss, ps, data, opt)
end


text:
using Keras
using PyCall
using Plots
@pyimport numpy

ip = numpy.load("testX.npy")
op = numpy.load("testY.npy")

op_obt = []
m = Keras.load("structure.json", "weights.h5")
println(m(ip[1,:,:]))
#for x=1:100
#    m = Keras.load("structure.json", "weights.h5")
#   push!(op_obt, m(ip[x,:,:])[1])
#end

#plot(op_obt)

max_stars_repo_path: examples/stock prediction/main.jl

max_issues_repo_path: examples/stock prediction/main.jl

max_forks_repo_path: examples/stock prediction/main.jl


text:
# Request handlers

function run_request(conn, state, msg_body, msg_id)
    @debug "run_request"

    state.debug_mode = :launch

    try
        include(msg_body)
    catch err
        Base.display_error(stderr, err, catch_backtrace())
    end

    send_notification(conn, "FINISHED")

    return :break
end

function debug_request(conn, state, msg_body, msg_id)
    @debug "debug_request"

    state.debug_mode = :launch
    index_of_sep = findfirst(';', msg_body)

    stop_on_entry_as_string = msg_body[1:index_of_sep-1]

    stop_on_entry = stop_on_entry_as_string=="stopOnEntry=true"

    filename_to_debug = msg_body[index_of_sep+1:end]

    @debug "We are debugging the file $filename_to_debug."

    ex = _parse_julia_file(filename_to_debug)

    state.top_level_expressions, _ = JuliaInterpreter.split_expressions(Main, ex)
    state.current_top_level_expression = 0

    state.frame = get_next_top_level_frame(state)

    if stop_on_entry
        send_notification(conn, "STOPPEDENTRY")
    elseif JuliaInterpreter.shouldbreak(state.frame, state.frame.pc)
        send_notification(conn, "STOPPEDBP")
    else
        ret = our_debug_command(:c, state)

        if ret===nothing
            send_notification(conn, "FINISHED")
            return :break
        else
            send_stopped_msg(conn, ret, state)
        end
    end

    return
end

function exec_request(conn, state, msg_body, msg_id)
    @debug "exec_request"

    state.debug_mode = :attach

    index_of_sep = findfirst(';', msg_body)

    stop_on_entry_as_string = msg_body[1:index_of_sep-1]

    stop_on_entry = stop_on_entry_as_string=="stopOnEntry=true"

    code_to_debug = msg_body[index_of_sep+1:end]

    state.sources[0] = code_to_debug

    ex = Meta.parse(code_to_debug)

    state.top_level_expressions, _ = JuliaInterpreter.split_expressions(Main, ex)
    state.current_top_level_expression = 0

    state.frame = get_next_top_level_frame(state)

    if stop_on_entry
        send_notification(conn, "STOPPEDENTRY")
    elseif JuliaInterpreter.shouldbreak(state.frame, state.frame.pc)
        send_notification(conn, "STOPPEDBP")
    else
        ret = our_debug_command(:c, state)

        if ret===nothing
            send_notification(conn, "FINISHED")
        else
            send_stopped_msg(conn, ret, state)
        end
    end
end

function setbreakpoints_request(conn, state, msg_body, msg_id)
    @debug "setbreakpoints_request"

    splitted_line = split(msg_body, ';')

    file = splitted_line[1]
    bps = map(split(i, ':') for i in splitted_line[2:end]) do i
        decoded_condition = String(Base64.base64decode(i[2]))
        # We handle conditions that don't parse properly as
        # no condition for now
        parsed_condition = try
            decoded_condition=="" ? nothing : Meta.parse(decoded_condition)
        catch err
            nothing
        end

        (line=parse(Int,i[1]), condition=parsed_condition)
    end

    for bp in JuliaInterpreter.breakpoints()
        if bp isa JuliaInterpreter.BreakpointFileLocation
            if bp.path==file
                JuliaInterpreter.remove(bp)
            end
        end
    end

    for bp in bps
        @debug "Setting one breakpoint at line $(bp.line) with condition $(bp.condition) in file $file."

        JuliaInterpreter.breakpoint(string(file), bp.line, bp.condition)
    end
end

function setexceptionbreakpoints_request(conn, state, msg_body, msg_id)
    @debug "setexceptionbreakpoints_request"

    opts = Set(split(msg_body, ';'))

    if "error" in opts
        JuliaInterpreter.break_on(:error)
    else
        JuliaInterpreter.break_off(:error)
    end

    if "throw" in opts
        JuliaInterpreter.break_on(:throw )
    else
        JuliaInterpreter.break_off(:throw )
    end

    if "compilemode" in opts
        state.compile_mode = JuliaInterpreter.Compiled()
    else
        state.compile_mode = JuliaInterpreter.finish_and_return!
    end
end

function setfunctionbreakpoints_request(conn, state, msg_body, msg_id)
    @debug "setfunctionbreakpoints_request"

    funcs = split(msg_body, ';', keepempty=false)

    bps = map(split(i, ':') for i in funcs) do i
        decoded_name = String(Base64.base64decode(i[1]))
        decoded_condition = String(Base64.base64decode(i[2]))

        parsed_condition = try
            decoded_condition=="" ? nothing : Meta.parse(decoded_condition)
        catch err
            nothing
        end

        try
            parsed_name = Meta.parse(decoded_name)

            if parsed_name isa Symbol
                return (mod=Main, name=parsed_name, signature=nothing, condition=parsed_condition)
            elseif parsed_name isa Expr
                if parsed_name.head==:.
                    # TODO Support this case
                    return nothing
                elseif parsed_name.head==:call
                    all_args_are_legit = true
                    if length(parsed_name.args)>1
                        for arg in parsed_name.args[2:end]
                            if !(arg isa Expr) || arg.head!=Symbol("::") || length(arg.args)!=1
                                all_args_are_legit =false
                            end
                        end
                        if all_args_are_legit

                            return (mod=Main, name=parsed_name.args[1], signature=map(j->j.args[1], parsed_name.args[2:end]), condition=parsed_condition)
                        else
                            return (mod=Main, name=parsed_name.args[1], signature=nothing, condition=parsed_condition)
                        end
                    else
                        return (mod=Main, name=parsed_name.args[1], signature=nothing, condition=parsed_condition)
                    end
                else
                    return nothing
                end
            else
                return nothing
            end
        catch err
            return nothing
        end

        return nothing
    end

    bps = filter(i->i!==nothing, bps)

    for bp in JuliaInterpreter.breakpoints()
        if bp isa JuliaInterpreter.BreakpointSignature
            JuliaInterpreter.remove(bp)
        end
    end

    state.not_yet_set_function_breakpoints = Set{Any}(bps)

    attempt_to_set_f_breakpoints!(state.not_yet_set_function_breakpoints)
end

function getstacktrace_request(conn, state, msg_body, msg_id)
    @debug "getstacktrace_request"

    fr = state.frame

    curr_fr = JuliaInterpreter.leaf(fr)

    frames_as_string = String[]

    id = 1
    while curr_fr!==nothing
        curr_scopeof = JuliaInterpreter.scopeof(curr_fr)
        curr_whereis = JuliaInterpreter.whereis(curr_fr)

        file_name = curr_whereis[1]
        lineno = curr_whereis[2]
        meth_or_mod_name = Base.nameof(curr_fr)

        if isfile(file_name)
            push!(frames_as_string, string(id, ";", meth_or_mod_name, ";path;", file_name, ";", lineno))
        elseif curr_scopeof isa Method
            ret = JuliaInterpreter.CodeTracking.definition(String, curr_fr.framecode.scope)
            if ret!==nothing
                state.sources[state.next_source_id], loc = ret
                s = string(id, ";", meth_or_mod_name, ";ref;", state.next_source_id, ";", file_name, ";", lineno)
                push!(frames_as_string, s)
                state.next_source_id += 1
            else
                src = curr_fr.framecode.src
                src = JuliaInterpreter.copy_codeinfo(src)
                JuliaInterpreter.replace_coretypes!(src; rev=true)
                code = Base.invokelatest(JuliaInterpreter.framecode_lines, src)

                state.sources[state.next_source_id] = join(code, '\n')

                s = string(id, ";", meth_or_mod_name, ";ref;", state.next_source_id, ";", file_name, ";", lineno)
                push!(frames_as_string, s)
                state.next_source_id += 1
            end
        else
            # For now we are assuming that this can only happen
            # for code that is passed via the @enter or @run macros,
            # and that code we have stored as source with id 0
            s = string(id, ";", meth_or_mod_name, ";ref;", 0, ";", "REPL", ";", lineno)
            push!(frames_as_string, s)
        end

        id += 1
        curr_fr = curr_fr.caller
    end

    send_response(conn, msg_id, join(frames_as_string, '\n'))
end

function getscope_request(conn, state::DebuggerState, msg_body, msg_id)
    @debug "getscope_request"

    frameId = parse(Int, msg_body)

    curr_fr = JuliaInterpreter.leaf(state.frame)

    i = 1

    while frameId > i
        curr_fr = curr_fr.caller
        i += 1
    end

    curr_scopeof = JuliaInterpreter.scopeof(curr_fr)
    curr_whereis = JuliaInterpreter.whereis(curr_fr)

    file_name = curr_whereis[1]
    code_range = curr_scopeof isa Method ? JuliaInterpreter.compute_corrected_linerange(curr_scopeof) : nothing

    push!(state.varrefs, VariableReference(:scope, curr_fr))

    var_ref_id = length(state.varrefs)

    if isfile(file_name) && code_range!==nothing
        send_response(conn, msg_id, "Local;$var_ref_id;$(code_range.start);$(code_range.stop);$file_name")
    else
        send_response(conn, msg_id, "Local;$var_ref_id")
    end
end

function getsource_request(conn, state, msg_body, msg_id)
    @debug "getsource_request"

    source_id = parse(Int, msg_body)

    send_response(conn, msg_id, state.sources[source_id])
end

function construct_return_msg_for_var(state::DebuggerState, name, value)
    v_type = typeof(value)
    v_type_as_string = string(v_type)
    v_value_as_string = Base.invokelatest(repr, value)
    v_value_encoded = Base64.base64encode(v_value_as_string)

    if (isstructtype(v_type) || value isa AbstractArray || value isa AbstractDict) && !(value isa String || value isa Symbol)
        push!(state.varrefs, VariableReference(:var, value))
        new_var_id = length(state.varrefs)

        named_count = if value isa Array || value isa Tuple
            0
        elseif value isa AbstractArray || value isa AbstractDict
            fieldcount(v_type) > 0 ? 1 : 0
        else
            fieldcount(v_type)
        end

        indexed_count = 0

        if value isa AbstractArray || value isa AbstractDict || value isa Tuple
            try
                indexed_count = Base.invokelatest(length, value)
            catch err
            end
        end

        return string(new_var_id, ";", name, ";", v_type_as_string, ";", named_count, ";", indexed_count, ";", v_value_encoded)
    else
        return string("0;", name, ";", v_type_as_string, ";0;0;", v_value_encoded)
    end
end

function get_keys_with_drop_take(value, skip_count, take_count)
    collect(Iterators.take(Iterators.drop(keys(value), skip_count), take_count))
end

function get_cartesian_with_drop_take(value, skip_count, take_count)
    collect(Iterators.take(Iterators.drop(CartesianIndices(value), skip_count), take_count))
end

function getvariables_request(conn, state::DebuggerState, msg_body, msg_id)
    @debug "getvariables_request"

    parts =split(msg_body, ';')

    var_ref_id = parse(Int, parts[1])

    filter_type = parts[2]
    skip_count = parts[3] == "" ? 0 : parse(Int, parts[3])
    take_count = parts[4] == "" ? typemax(Int) : parse(Int, parts[4])

    var_ref = state.varrefs[var_ref_id]

    vars_as_string = String[]

    if var_ref.kind==:scope
        curr_fr = var_ref.value

        vars = JuliaInterpreter.locals(curr_fr)

        for v in vars
            # TODO Figure out why #self# is here in the first place
            # For now we don't report it to the client
            if !startswith(string(v.name), "#") && string(v.name)!=""
                s = construct_return_msg_for_var(state, string(v.name), v.value)
                push!(vars_as_string, s)
            end
        end

        if JuliaInterpreter.isexpr(JuliaInterpreter.pc_expr(curr_fr), :return)
            ret_val = JuliaInterpreter.get_return(curr_fr)
            s = construct_return_msg_for_var(state, "Return Value", ret_val)

            push!(vars_as_string, s)
        end
    elseif var_ref.kind==:var
        container_type = typeof(var_ref.value)

        if filter_type=="" || filter_type=="named"
            if (var_ref.value isa AbstractArray || var_ref.value isa AbstractDict) && !(var_ref.value isa Array) &&
                fieldcount(container_type) > 0
                push!(state.varrefs, VariableReference(:fields, var_ref.value))
                new_var_id = length(state.varrefs)
                named_count = fieldcount(container_type)
                s = string(new_var_id, ";Fields;;", named_count, ";0;", Base64.base64encode(""))

                push!(vars_as_string, s)
            else
                for i=Iterators.take(Iterators.drop(1:fieldcount(container_type), skip_count), take_count)
                    s = construct_return_msg_for_var(state, string(fieldname(container_type, i)), getfield(var_ref.value, i) )
                    push!(vars_as_string, s)
                end
            end
        end

        if (filter_type=="" || filter_type=="indexed")
            if var_ref.value isa Tuple
                for i in Iterators.take(Iterators.drop(1:length(var_ref.value), skip_count), take_count)
                    s = construct_return_msg_for_var(state, join(string.(i), ','), var_ref.value[i])
                    push!(vars_as_string, s)
                end
            elseif var_ref.value isa AbstractArray
                for i in Base.invokelatest(get_cartesian_with_drop_take, var_ref.value, skip_count, take_count)
                    try
                        val = Base.invokelatest(getindex, var_ref.value, i)
                        s = construct_return_msg_for_var(state, join(string.(i.I), ','), val)
                    catch err
                        s = string("0;", join(string.(i.I), ','), ";;0;0;", Base64.base64encode("#error"))
                    end
                    push!(vars_as_string, s)
                end
            elseif var_ref.value isa AbstractDict
                for i in Base.invokelatest(get_keys_with_drop_take, var_ref.value, skip_count, take_count)
                    key_as_string = Base.invokelatest(repr, i)
                    try
                        val = Base.invokelatest(getindex, var_ref.value, i)
                        s = construct_return_msg_for_var(state, key_as_string, val)
                    catch err
                        s = string("0;", join(string.(i.I), ','), ";;0;0;", Base64.base64encode("#error"))
                    end
                    push!(vars_as_string, s)
                end
            end
        end
    elseif var_ref.kind==:fields
        container_type = typeof(var_ref.value)

        if filter_type=="" || filter_type=="named"
            for i=Iterators.take(Iterators.drop(1:fieldcount(container_type), skip_count), take_count)
                s = construct_return_msg_for_var(state, string(fieldname(container_type, i)), getfield(var_ref.value, i) )
                push!(vars_as_string, s)
            end
        end

    end

    send_response(conn, msg_id, join(vars_as_string, '\n'))
end

function setvariable_request(conn, state::DebuggerState, msg_body, msg_id)
    parts = split(msg_body, ';')

    varref_id = parse(Int, parts[1])
    var_name = String(Base64.base64decode(parts[2]))
    var_value = String(Base64.base64decode(parts[3]))

    val_parsed = try
        parsed = Meta.parse(var_value)

        if parsed isa Expr && !(parsed.head==:call || parsed.head==:vect || parsed.head==:tuple)
            send_response(conn, msg_id, "FAILED;$(Base64.base64encode("Only values or function calls are allowed."))")
            return
        end

        parsed
    catch err
        send_response(conn, msg_id, "FAILED;$(Base64.base64encode("Something went wrong in the eval."))")
        return
    end

    var_ref = state.varrefs[varref_id]

    if var_ref.kind==:scope
        try
            ret = JuliaInterpreter.eval_code(var_ref.value, "$var_name = $var_value");

            s = construct_return_msg_for_var(state::DebuggerState, "", ret)

            send_response(conn, msg_id, "SUCCESS;$s")
            return
        catch err
            send_response(conn, msg_id, "FAILED;$(Base64.base64encode("Something went wrong in the set: $err"))")
            return
        end
    elseif var_ref.kind==:var
        if isnumeric(var_name[1])
            try
                new_val = try
                    Core.eval(Main, val_parsed)
                catch err
                    send_response(conn, msg_id, "FAILED;$(Base64.base64encode("Expression cannot be evaluated."))")
                    return
                end

                idx = Core.eval(Main, Meta.parse("($var_name)"))

                setindex!(var_ref.value, new_val, idx...)

                s = construct_return_msg_for_var(state::DebuggerState, "", new_val)

                send_response(conn, msg_id, "SUCCESS;$s")
                return
            catch err
                send_response(conn, msg_id, "FAILED;$(Base64.base64encode("Something went wrong in the set: $err"))")
                return
            end
        else
            if Base.isimmutable(var_ref.value)
                send_response(conn, msg_id, "FAILED;$(Base64.base64encode("Cannot change the fields of an immutable struct."))")
                return
            else
                try
                    new_val = try
                        Core.eval(Main, val_parsed)
                    catch err
                        send_response(conn, msg_id, "FAILED;$(Base64.base64encode("Expression cannot be evaluated."))")
                        return
                    end

                    setfield!(var_ref.value, Symbol(var_name), new_val)

                    s = construct_return_msg_for_var(state::DebuggerState, "", new_val)

                    send_response(conn, msg_id, "SUCCESS;$s")
                    return
                catch err
                    send_response(conn, msg_id, "FAILED;$(Base64.base64encode("Something went wrong in the set: $err"))")
                    return
                end
            end
        end
    else
        error("Unknown var ref type.")
    end
end

function restartframe_request(conn, state::DebuggerState, msg_body, msg_id)
    frame_id = parse(Int, msg_body)

    curr_fr = JuliaInterpreter.leaf(state.frame)

    i = 1

    while frame_id > i
        curr_fr = curr_fr.caller
        i += 1
    end

    if curr_fr.caller===nothing
        # We are in the top level

        state.current_top_level_expression = 0

        state.frame = get_next_top_level_frame(state)
    else
        curr_fr.pc = 1
        curr_fr.assignment_counter = 1
        curr_fr.callee = nothing

        state.frame = curr_fr
    end

    ret = our_debug_command(:c, state)

    if ret===nothing
        send_notification(conn, "FINISHED")
        state.debug_mode==:launch && return :break
    else
        send_stopped_msg(conn, ret, state)
    end

    return
end

function getexceptioninfo_request(conn, state, msg_body, msg_id)
    exception_id = string(typeof(state.last_exception))
    exception_description = sprint(Base.showerror, state.last_exception)

    exception_stacktrace = sprint(Base.show_backtrace, state.frame)

    send_response(conn, msg_id, string(Base64.base64encode(exception_id), ';', Base64.base64encode(exception_description), ';', Base64.base64encode(exception_stacktrace)))
end

function evaluate_request(conn, state, msg_body, msg_id)
    @debug "evaluate_request"

    index_of_sep = findfirst(':', msg_body)

    stack_id = parse(Int, msg_body[1:index_of_sep-1])

    expression = msg_body[index_of_sep+1:end]

    curr_fr = state.frame
    curr_i = 1

    while stack_id > curr_i
        if curr_fr.caller!==nothing
            curr_fr = curr_fr.caller
            curr_i += 1
        else
            break
        end
    end

    try
        ret_val = JuliaInterpreter.eval_code(curr_fr, expression)

        send_response(conn, msg_id, Base.invokelatest(repr, ret_val))
    catch err
        send_response(conn, msg_id, "#error")
    end
end

function continue_request(conn, state, msg_body, msg_id)
    @debug "continue_request"

    ret = our_debug_command(:c, state)

    if ret===nothing
        send_notification(conn, "FINISHED")
        state.debug_mode==:launch && return :break
    else
        send_stopped_msg(conn, ret, state)
    end

    return
end

function next_request(conn, state, msg_body, msg_id)
    @debug "next_request"

    ret = our_debug_command(:n, state)

    if ret===nothing
        send_notification(conn, "FINISHED")
        state.debug_mode==:launch && return :break
    else
        send_stopped_msg(conn, ret, state)
    end

    return
end

function stepin_request(conn, state, msg_body, msg_id)
    @debug "stepin_request"

    ret = our_debug_command(:s, state)

    if ret===nothing
        send_notification(conn, "FINISHED")
        state.debug_mode==:launch && return :break
    else
        send_stopped_msg(conn, ret, state)
    end

    return
end

function stepout_request(conn, state, msg_body, msg_id)
    @debug "stepout_request"

    ret = our_debug_command(:finish, state)

    if ret===nothing
        send_notification(conn, "FINISHED")
        state.debug_mode==:launch && return :break
    else
        send_stopped_msg(conn, ret, state)
    end

    return
end

function disconnect_request(conn, state, msg_body, msg_id)
    @debug "disconnect_request"

    return :break
end

function terminate_request(conn, state, msg_body, msg_id)
    @debug "terminate_request"

    send_notification(conn, "FINISHED")
    return :break
end

max_stars_repo_name: caleb-allen/julia-vscode

max_issues_repo_name: caleb-allen/julia-vscode

max_forks_repo_name: caleb-allen/julia-vscode


text:
module Pebblegames

using Graphtypes
using Graphsearch

export
    solve_basic_pebble_game,
    basic_pebble_game,
    add_edge!,
    slide_move!,
    go_get_pebbles!

type PebbleGame
    graph::Graph
    pebbles::Array{Int64,1}
    k::Int64
    l::Int64
end

#=
    BASIC (k, l) PEBBLE GAME
    from Pebble Game Algorithms and Sparse Graphs - Lee & Streinu
    http://arxiv.org/pdf/math/0702129.pdf

    Rules:
        - vertex can hold up to k pebbles
        - new edge is accepteble if at least l+1 pebbles at endpoints
    Moves:
        - slide move: if there's at a connection i,j and at least one
            pebble on j, reverse i,j and move a pebble from j to i
        - add edge move: add an acceptable edge i,j and remove a
        pebble from i
=#
function solve_basic_pebble_game(graph::Graph, k::Int64, l::Int64)

    vt= total_vertices(graph)
    et = total_edges(graph)

    game = basic_pebble_game(vt, k, l)
    edges = list_edges(graph)
    rejected_count = 0
    for n=1:et
        i = edges[n,1]
        j = edges[n,2]

        if add_edge!(game, i, j)
            #successfully added an edge, move on
            continue
        else
            if game.pebbles[i] + game.pebbles[j] < game.l+1
                go_get_pebbles!(game, i, [i, j])
                go_get_pebbles!(game, j, [i, j])
            end
            if game.pebbles[i] + game.pebbles[j] < game.l+1
                rejected_count += 1
                continue
            end

            add_edge!(game, i, j)
    end end
    return sum(game.pebbles), rejected_count
end

function basic_pebble_game (s::Int64, k::Int64, l::Int64)
    g = simple_graph(s)
    p = PebbleGame(g, k*ones(s), k ,l)
    return p
end

function add_edge!(p::PebbleGame, i::Int64, j::Int64)
    if p.pebbles[i] + p.pebbles[j] < p.l+1
        return false
    end
    p.pebbles[i] -= 1
    connect!(p.graph,i,j)
    return true
end

function slide_move!(p::PebbleGame, i::Int64, j::Int64)
    if !connected(p.graph,i,j)
        return false
    end
    if p.pebbles[j]<1
        return false
    end
    if p.pebbles[i]>=p.k
        return false
    end
    reverse!(p.graph,i,j)
    p.pebbles[j]-=1
    p.pebbles[i]+=1
    return true
end

function go_get_pebbles!(game::PebbleGame, node::Int64,
                            blacklist::Array{Int64,1})
    graph = game.graph
    pebbles = game.pebbles
    #helper function that looks for available pebbles and moves them to node
    for i=1:(min(game.l+1,game.k) - pebbles[node]) #find missing pebbles for i
        path = depth_first_search(graph, node, total_vertices(graph),
                                (g,x)->(pebbles[x]>0), blacklist)

        path = reverse(path)
        #go fetch that pebble
        for pair in zip(path,path[2:end])
            x,y = pair
            slide_move!(game,y,x)
end end end

end

max_stars_repo_name: AaltoScienceInst/pebblegame

max_issues_repo_name: AaltoScienceInst/pebblegame

max_forks_repo_name: AaltoScienceInst/pebblegame


text:
# ---
# title: 373. Find K Pairs with Smallest Sums
# id: problem373
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Heap
# link: <https://leetcode.com/problems/find-k-pairs-with-smallest-sums/description/>
# hidden: true
# ---
#
# You are given two integer arrays **nums1** and **nums2** sorted in ascending
# order and an integer **k**.
#
# Define a pair **(u,v)** which consists of one element from the first array and
# one element from the second array.
#
# Find the k pairs **(u 1,v1),(u2,v2) ...(uk,vk)** with the smallest sums.
#
# **Example 1:**
#
#
#
#     Input: nums1 = [1,7,11], nums2 = [2,4,6], k = 3
#     Output: [[1,2],[1,4],[1,6]]
#     Explanation: The first 3 pairs are returned from the sequence:
#                  [1,2],[1,4],[1,6],[7,2],[7,4],[11,2],[7,6],[11,4],[11,6]
#
# **Example 2:**
#
#
#
#     Input: nums1 = [1,1,2], nums2 = [1,2,3], k = 2
#     Output: [1,1],[1,1]
#     Explanation: The first 2 pairs are returned from the sequence:
#                  [1,1],[1,1],[1,2],[2,1],[1,2],[2,2],[1,3],[1,3],[2,3]
#
# **Example 3:**
#
#
#
#     Input: nums1 = [1,2], nums2 = [3], k = 3
#     Output: [1,3],[2,3]
#     Explanation: All possible pairs are returned from the sequence: [1,3],[2,3]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
# module Crystal

# Examples

```jldoctest
julia>
```
"""
module Crystal

using StaticArrays: SVector, SMatrix, FieldVector

export CellParameters,
    CrystalFamily,
    Cubic, Hexagonal, Tetragonal, Orthorhombic, Monoclinic, Triclinic,
    CenteringType,
    Primitive, FaceCentered, BodyCentered, BaseCentered, RhombohedrallyCentered,
    CUBIC_LATTICE_SYSTEM,
    HEXAGONAL_LATTICE_SYSTEM,
    RHOMBOHEDRAL_LATTICE_SYSTEM,
    TETRAGONAL_LATTICE_SYSTEM,
    ORTHORHOMBIC_LATTICE_SYSTEM,
    MONOCLINIC_LATTICE_SYSTEM,
    TRICLINIC_LATTICE_SYSTEM

struct CellParameters
    coordinates::SMatrix{3, 3, Float64}
end

abstract type CrystalFamily end

struct Cubic <: CrystalFamily end
struct Hexagonal <: CrystalFamily end
struct Tetragonal <: CrystalFamily end
struct Orthorhombic <: CrystalFamily end
struct Monoclinic <: CrystalFamily end
struct Triclinic <: CrystalFamily end

abstract type CenteringType end

struct Primitive <: CenteringType end
struct FaceCentered <: CenteringType end
struct BodyCentered <: CenteringType end
struct BaseCentered <: CenteringType end
struct RhombohedrallyCentered <: CenteringType end

struct CellDimension <: FieldVector{6, Float64}
    a::Float64
    b::Float64
    c::Float64
    cosγ::Float64
    cosβ::Float64
    cosα::Float64
end

struct BravaisLattice{A <: CrystalFamily, B <: CenteringType}
    ibrav::Int
    celldm::CellDimension
end

BravaisLattice{A, B}(ibrav::Int) where {A, B} = celldm::AbstractArray{Float64, 1} -> BravaisLattice{A, B}(ibrav, CellDimension(celldm))

const CUBIC_LATTICE_SYSTEM = (BravaisLattice{Cubic, Primitive}(1), BravaisLattice{Cubic, FaceCentered}(2), BravaisLattice{Cubic, BodyCentered}(3))
const HEXAGONAL_LATTICE_SYSTEM = (BravaisLattice{Hexagonal, Primitive}(4),)
const RHOMBOHEDRAL_LATTICE_SYSTEM = (BravaisLattice{Hexagonal, RhombohedrallyCentered}(5), BravaisLattice{Hexagonal, RhombohedrallyCentered}(-5))
const TETRAGONAL_LATTICE_SYSTEM = (BravaisLattice{Tetragonal, Primitive}(6), BravaisLattice{Tetragonal, BodyCentered}(7))
const ORTHORHOMBIC_LATTICE_SYSTEM = (BravaisLattice{Orthorhombic, Primitive}(8), BravaisLattice{Orthorhombic, BaseCentered}(9),
    BravaisLattice{Orthorhombic, FaceCentered}(10), BravaisLattice{Orthorhombic, BodyCentered}(11))
const MONOCLINIC_LATTICE_SYSTEM = (BravaisLattice{Monoclinic, Primitive}(12), BravaisLattice{Monoclinic, Primitive}(-12),
    BravaisLattice{Monoclinic, BaseCentered}(13))
const TRICLINIC_LATTICE_SYSTEM = (BravaisLattice{Triclinic, Primitive}(14),)

function convert(b::BravaisLattice)::CellParameters
    ibrav = b.ibrav

end

[celldm.a 0 0; celldm.b * celldm.cosγ celldm.b * celldm.sinγ 0; 0 0 celldm.c]

end


text:
module Environment
export Bandit, step

struct Bandit
    p::Array{Real, 1}
end

function step(env::Bandit, action::Int64)
    k = length(env.p)
    rewards = env.p .>= rand(k)
    return rewards[action]
end

end#module


text:
using ReservoirComputing

#model parameters
const approx_res_size = 30
const radius = 1.2
const activation = tanh
const degree = 6
const sigma = 0.1
const beta = 0.0
const alpha = 1.0
nla_type = NLADefault()
const in_size = 3
const out_size = 3
const extended_states = false

const train_len = 50
const predict_len = 12
data = ones(Float64, in_size, 100)
train = data[:, 1:1+train_len-1]
test = data[:, train_len:train_len+predict_len-1]

#test non linear algos
nla = [NLAT1(), NLAT2(), NLAT3()]
for t in nla
    nla_type = t
    esn = ESN(approx_res_size,
        train,
        degree,
        radius,
        activation = activation,
        sigma = sigma,
        alpha = alpha,
        nla_type = nla_type,
        extended_states = extended_states)

    W_out = ESNtrain(esn, beta)
    @test size(W_out) == (out_size, esn.res_size)
    output = ESNpredict(esn, predict_len, W_out)
    @test size(output) == (out_size, predict_len)

end


text:
using MixedModels, Test

# include("src/RegressionTables.jl")

m1 = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false)
gm1 = fit(MixedModel, @formula(use ~ 1 + urban + livch + age + abs2(age) + (1|dist)),
              MixedModels.dataset(:contra), Bernoulli(); progress=false)
ni1 = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff), Normal(); progress=false)
gm2 = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff),
    Normal(), SqrtLink(); progress=false)

RegressionTables.regtable(m1,gm1,ni1,gm2; renderSettings = RegressionTables.asciiOutput(joinpath(dirname(@__FILE__), "tables", "mmtest1.txt")), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof],estimformat = "%0.1f", statisticformat = "%0.1f")
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "mmtest1.txt"), joinpath(dirname(@__FILE__), "tables", "mmtest1_reference.txt"))

rm(joinpath(dirname(@__FILE__), "tables", "mmtest1.txt"))


text:
module TwoLink
using FixedSizeArrays

export torque, forward_kin, inverse_kin, inverse_kin_up, inverse_kin_down, traj, V2, connect_points

typealias V2 Vec{2,Float64}

const v1 = v2 = k1 = k2 = m1 = m2 = l1 = l2 = 1
const g = 9.82

function torque(q,qd,qdd)

    q1,q2 = q[1],q[2]
    qd1,qd2 = qd[1],qd[2]
    qdd1,qdd2 = qdd[1],qdd[2]

    c1 = cos(q1)
    c2 = cos(q2)
    s1 = sin(q1)
    s2 = sin(q2)
    s12 = sin(q1+q2)

    τ1 = m2*l2^2*(qdd1+qdd2) + m2*l1*l2*c2*(2qdd1+qdd2) +
    (m1+m2)*l1^2+qdd1 - m2*l1*l2*s2*qd2^2 - 2m2*l1*l2*s2*qd1*qd2 +
    m2*l2*g*s12 + (m1+m2)*l1*g*s1 + v1*qd1 +
    k1*sign(qd1)

    τ2 = m2*l1*l2*c2*qdd1 + m2*l1*l2*s2*qd1^2 + m2*l2*g*s12 +
    m2*l2^2*(qdd1 + qdd2) + v2*qd2 + k2*sign(qd2)

    V2(τ1,τ2)

end

function forward_kin(q)
    q1,q2 = q[1],(q[1]+q[2])
    p1 = V2(l1*sin(q1), -l1*cos(q1))
    p = p1 + V2(l2*sin(q2), -l2*cos(q2))
    p1, p
end

function forward_kin(jtraj::Matrix)
    p = zeros(jtraj)
    for i = 1:size(jtraj,1)
        q = jtraj[i,:][:]
        p[i,:] = forward_kin(q)[2] |> collect
    end
    p
end

function inverse_kin(p)
    x,y = p[1], p[2]
    a2 = x^2 + y^2
    a = sqrt(a2)
    α = acos((l1^2 + l2^2 - a2)/(2l1*l2))
    α2 = 2π-α
    q21 = π - α
    q22 = π - α2

    β = asin(sin(-α)/a*l2)
    β2 = -β
    q11 = asin(x/a)-β2
    q12 = asin(x/a)-β
    V2(q11,q21), V2(q12,q22)
end

function inverse_kin_up(p)
    q1,q2 = inverse_kin(p)
    return q1[2] > q2[2] ? q1 : q2
end

function inverse_kin_down(p)
    q1,q2 = inverse_kin(p)
    return q1[2] < q2[2] ? q1 : q2
end

function inverse_kin(ctraj::Matrix, dir = :up)
    if dir == :up
        kin_fun = inverse_kin_up
    elseif dir == :down
        kin_fun = inverse_kin_down
    else
        error("Symbol direction unknown")
    end

    p = zeros(ctraj)
    for i = 1:size(ctraj,1)
        q = ctraj[i,:][:]
        p[i,:] = kin_fun(q) |> collect
    end
    p
end

function traj(q0,q1,t)
    tf = maximum(t)
    V = (q1-q0)/tf * 1.5
    traj(q0,q1,t, V)
end

function traj(q0,q1,t, V)

    tf = maximum(t)

    V = abs(V) * sign(q1-q0)
    if abs(V) < abs(q1-q0)/tf
        error("V too small")
    elseif abs(V) > 2*abs(q1-q0)/tf
        error("V too big")
    end

    if q0 == q1
        s = ones(Float64,size(t)) * q0
        sd = zeros(Float64,size(t))
        sdd = zeros(Float64,size(t))
        return s, sd, sdd
    end

    tb = (q0 - q1 + V*tf)/V
    a = V/tb

    p = Array(Float64,size(t))
    pd = Array(Float64,size(t))
    pdd = Array(Float64,size(t))

    for (i,t) = enumerate(t)
        if t <= tb
            # initial blend
            p[i] = q0 + a/2*t^2
            pd[i] = a*t
            pdd[i] = a
        elseif t <= tf-tb
            # linear motion
            p[i] = (q1+q0-V*tf)/2 + V*t
            pd[i] = V
            pdd[i] = 0
        else
            # final blend
            p[i] = q1 - a/2*tf^2 + a*tf*t - a/2*t^2
            pd[i] = a*tf - a*t;
            pdd[i] = -a;
        end
    end
    return p, pd, pdd

end

function connect_points(points,ni)
    fx(i) = traj(points[i,1],points[i+1,1],1:ni)
    fy(i) = traj(points[i,2],points[i+1,2],1:ni)
    n = size(points,1)
    x,xd,xdd = fx(1)
    y,yd,ydd = fy(1)
    p   = [x y]
    pd  = [xd yd]
    pdd = [xdd ydd]
    for i = 2:n-1
        x,xd,xdd = fx(i)
        y,yd,ydd = fy(i)
        p = cat(1,p,[x y])
        pd = cat(1,pd,[xd yd])
        pdd = cat(1,pdd,[xdd ydd])

    end
    p, pd, pdd

end

end


text:
# Inspired from example1.cpp from SDPA
p = SDPA.SDPAProblem()
mDIM   = 5
nBlock = 3
SDPA.inputConstraintNumber(p, mDIM)
SDPA.inputBlockNumber(p, nBlock)
SDPA.inputBlockSize(p, 1, 2)
SDPA.inputBlockSize(p, 2, 3)
SDPA.inputBlockSize(p, 3, -2)
SDPA.inputBlockType(p, 1, SDPA.SDP)
SDPA.inputBlockType(p, 2, SDPA.SDP)
SDPA.inputBlockType(p, 3, SDPA.LP)

SDPA.initializeUpperTriangleSpace(p)

SDPA.inputCVec(p, 1,1.1)
SDPA.inputCVec(p, 2,-10)
SDPA.inputCVec(p, 3,6.6)
SDPA.inputCVec(p, 4,19)
SDPA.inputCVec(p, 5,4.1)

# Objective
SDPA.inputElement(p, 0, 1, 1, 1, -1.4, false)
SDPA.inputElement(p, 0, 1, 1, 2, -3.2, false)
SDPA.inputElement(p, 0, 1, 2, 2,  -28, false)

SDPA.inputElement(p, 0, 2, 1, 1,  15, false)
SDPA.inputElement(p, 0, 2, 1, 2, -12, false)
SDPA.inputElement(p, 0, 2, 1, 3, 2.1, false)
SDPA.inputElement(p, 0, 2, 2, 2,  16, false)
SDPA.inputElement(p, 0, 2, 2, 3,-3.8, false)
SDPA.inputElement(p, 0, 2, 3, 3,  15, false)

SDPA.inputElement(p, 0, 3, 1, 1, 1.8, false)
SDPA.inputElement(p, 0, 3, 2, 2,-4.0, false)

# 1st constraint
SDPA.inputElement(p, 1, 1, 1, 1, 0.5, false)
SDPA.inputElement(p, 1, 1, 1, 2, 5.2, false)
SDPA.inputElement(p, 1, 1, 2, 2,-5.3, false)

SDPA.inputElement(p, 1, 2, 1, 1, 7.8, false)
SDPA.inputElement(p, 1, 2, 1, 2,-2.4, false)
SDPA.inputElement(p, 1, 2, 1, 3, 6.0, false)
SDPA.inputElement(p, 1, 2, 2, 2, 4.2, false)
SDPA.inputElement(p, 1, 2, 2, 3, 6.5, false)
SDPA.inputElement(p, 1, 2, 3, 3, 2.1, false)

SDPA.inputElement(p, 1, 3, 1, 1, -4.5, false)
SDPA.inputElement(p, 1, 3, 2, 2, -3.5, false)

# 2nd constraint
SDPA.inputElement(p, 2, 1, 1, 1,  1.7, false)
SDPA.inputElement(p, 2, 1, 1, 2,  7.0, false)
SDPA.inputElement(p, 2, 1, 2, 2, -9.3, false)

SDPA.inputElement(p, 2, 2, 1, 1, -1.9, false)
SDPA.inputElement(p, 2, 2, 1, 2, -0.9, false)
SDPA.inputElement(p, 2, 2, 1, 3, -1.3, false)
SDPA.inputElement(p, 2, 2, 2, 2, -0.8, false)
SDPA.inputElement(p, 2, 2, 2, 3, -2.1, false)
SDPA.inputElement(p, 2, 2, 3, 3,  4.0, false)

SDPA.inputElement(p, 2, 3, 1, 1, -0.2, false)
SDPA.inputElement(p, 2, 3, 2, 2, -3.7, false)

# 3rd constraint
SDPA.inputElement(p, 3, 1, 1, 1,  6.3, false)
SDPA.inputElement(p, 3, 1, 1, 2, -7.5, false)
SDPA.inputElement(p, 3, 1, 2, 2, -3.3, false)

SDPA.inputElement(p, 3, 2, 1, 1,  0.2, false)
SDPA.inputElement(p, 3, 2, 1, 2,  8.8, false)
SDPA.inputElement(p, 3, 2, 1, 3,  5.4, false)
SDPA.inputElement(p, 3, 2, 2, 2,  3.4, false)
SDPA.inputElement(p, 3, 2, 2, 3, -0.4, false)
SDPA.inputElement(p, 3, 2, 3, 3,  7.5, false)

SDPA.inputElement(p, 3, 3, 1, 1, -3.3, false)
SDPA.inputElement(p, 3, 3, 2, 2, -4.0, false)

# 4th constraint
SDPA.inputElement(p, 4, 1, 1, 1, -2.4, false)
SDPA.inputElement(p, 4, 1, 1, 2, -2.5, false)
SDPA.inputElement(p, 4, 1, 2, 2, -2.9, false)

SDPA.inputElement(p, 4, 2, 1, 1,  3.4, false)
SDPA.inputElement(p, 4, 2, 1, 2, -3.2, false)
SDPA.inputElement(p, 4, 2, 1, 3, -4.5, false)
SDPA.inputElement(p, 4, 2, 2, 2,  3.0, false)
SDPA.inputElement(p, 4, 2, 2, 3, -4.8, false)
SDPA.inputElement(p, 4, 2, 3, 3,  3.6, false)

SDPA.inputElement(p, 4, 3, 1, 1, 4.8, false)
SDPA.inputElement(p, 4, 3, 2, 2, 9.7, false)

# 5th constraint
SDPA.inputElement(p, 5, 1, 1, 1, -6.5, false)
SDPA.inputElement(p, 5, 1, 1, 2, -5.4, false)
SDPA.inputElement(p, 5, 1, 2, 2, -6.6, false)

SDPA.inputElement(p, 5, 2, 1, 1,  6.7, false)
SDPA.inputElement(p, 5, 2, 1, 2, -7.2, false)
SDPA.inputElement(p, 5, 2, 1, 3, -3.6, false)
SDPA.inputElement(p, 5, 2, 2, 2,  7.3, false)
SDPA.inputElement(p, 5, 2, 2, 3, -3.0, false)
SDPA.inputElement(p, 5, 2, 3, 3, -1.4, false)

SDPA.inputElement(p, 5, 3, 1, 1, 6.1, false)
SDPA.inputElement(p, 5, 3, 2, 2,-1.5, false)

SDPA.initializeUpperTriangle(p, false)
SDPA.initializeSolve(p)

SDPA.solve(p)

# See Section 6.3 of
# SDPA (SemiDefinite Programming Algorithm) User's Manual — Version 6.2.0
# https://pdfs.semanticscholar.org/0332/d0044b09e1212e181bc422f390de05df0c88.pdf Section 6.3
# for the values
@test SDPA.getIteration(p) == 13
@test isapprox(SDPA.getPrimalObj(p), 32.06269340482402)
@test isapprox(SDPA.getDualObj(p), 32.062692353573865)
@test SDPA.getPrimalError(p) < 1e-10
@test SDPA.getDualError(p) < 1e-10

X = SDPA.VarDualSolution(p)
@test isapprox(SDPA.block(X, 1), [+6.392e-08 -9.638e-09;
                      -9.638e-09 +4.539e-08], rtol=1e-4)
@test isapprox(SDPA.block(X, 2), [+7.119e+00 +5.025e+00 +1.916e+00;
                      +5.025e+00 +4.415e+00 +2.506e+00;
                      +1.916e+00 +2.506e+00 +2.048e+00], rtol=1e-4)
@test isapprox(SDPA.block(X, 3), Diagonal([+3.432e-01, +4.391e+00]), rtol=1e-4)

Y = SDPA.PrimalSolution(p)
@test size(Y) == (7, 7)
@test Y[1, 3] == 0
@test_throws BoundsError Y[1, 0]
@test_throws BoundsError SDPA.block(Y, 0)
@test_throws BoundsError SDPA.block(Y, 4)
@test_throws BoundsError SDPA.block(Y, 1)[0, 1]
@test_throws BoundsError SDPA.block(Y, 1)[1, 3]
@test isapprox(SDPA.block(Y, 1), [+2.640e+00 +5.606e-01;
                      +5.606e-01 +3.718e+00], rtol=1e-4)
@test isapprox(Y[6, 6], +4.087e-07, rtol=1e-4)
@test isapprox(SDPA.block(Y, 2), [+7.616e-01 -1.514e+00 +1.139e+00;
                      -1.514e+00 +3.008e+00 -2.264e+00;
                      +1.139e+00 -2.264e+00 +1.705e+00], rtol=1e-3)
@test isapprox(Y[4, 3], -1.514e+00, rtol=1e-3)
@test isapprox(SDPA.block(Y, 3), Diagonal([+4.087e-07, +3.195e-08]), rtol=1e-4)
@test isapprox(Y[6, 6], +4.087e-07, rtol=1e-4)
@test Y[6, 7] == 0

SDPA.terminate(p)