text:
function update_eta!(i::Int, j::Int, s::State, c::Constants, d::Data)
  counts = Dict(z => zeros(c.L[z]) for z in 0:1)

  for n in 1:d.N[i]
    k = s.lam[i][n]
    if k > 0
      z = s.Z[j, k]
      l = s.gam[i][n, j]
      counts[z][l] += 1
    end
  end

  updatedParam0 = c.eta_prior[0].alpha .+ counts[0]
  updatedParam1 = c.eta_prior[1].alpha .+ counts[1]

  s.eta[0][i, j, :] = rand(Dirichlet(updatedParam0))
  s.eta[1][i, j, :] = rand(Dirichlet(updatedParam1))
end

function update_eta!(s::State, c::Constants, d::Data)
  for i in 1:d.I
    for j in 1:d.J
      update_eta!(i, j, s, c, d)
    end
  end
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
using ModelingToolkit
using Test

const PType = ModelingToolkit.Parameter{Number}

@parameters t
@variables x(t) y(t) # test multi-arg
@variables z(t) # test single-arg
x1 = Variable(:x)(t)
y1 = Variable(:y)(t)
z1 = Variable(:z)(t)
@test isequal(x1, x)
@test isequal(y1, y)
@test isequal(z1, z)

@parameters begin
    t
    s
end
@parameters σ(..)

t1 = Variable{PType}(:t)()
s1 = Variable{PType}(:s)()
σ1 = Variable(:σ)
@test isequal(t1, t)
@test isequal(s1, s)
@test isequal(σ1, σ)

@derivatives D'~t
D1 = Differential(t)
@test D1 == D

@test isequal(x ≤ y + 1, (x < y + 1) | (x == y + 1))

@test @macroexpand(@parameters x, y, z(t)) == @macroexpand(@parameters x y z(t))
@test @macroexpand(@variables x, y, z(t)) == @macroexpand(@variables x y z(t))

# Test array expressions
@parameters begin
    t[1:2]
    s[1:2:4,1:2]
end
@parameters σ[1:2](..)

t1 = [Variable{PType}(:t, 1)(),
      Variable{PType}(:t, 2)()]
s1 = [Variable{PType}(:s, 1, 1)() Variable{PType}(:s, 1, 2)()
      Variable{PType}(:s, 3, 1)() Variable{PType}(:s, 3, 2)()]
σ1 = [Variable(:σ, 1),
      Variable(:σ, 2)]
@test isequal(t1, t)
@test isequal(s1, s)
@test isequal(σ1, σ)

@parameters t
@variables x[1:2](t)
x1 = [Variable(:x, 1)(t), Variable(:x, 2)(t)]
@test isequal(x1, x)

@variables a[1:11,1:2]
@variables a()

ext: jl


text:
#=
# Uses the elemental order to determine the nodal order.
# Order within elements is:
#   interior -> empty sides(z=0,y=0,x=0,z=1,y=1,x=1) or (y=0,x=0,y=1,x=1)
=#
export reorder_grid_element_first, get_element_first_order

function element_first_order_2d(grid, porder, elorder, invert = true)
    nnodes = size(grid.allnodes,2);
    eforder = zeros(Int, nnodes); # The ordering
    node_done = zeros(Bool, nnodes);

    # Loop over the elements in elorder, filling in the eforder as we go
    l2g = grid.loc2glb;
    efind = 0;
    for ei=elorder
        # Add interior first
        for j = 2:porder
            for i = 2:porder
                lind = i + (porder+1)*(j-1);
                gind = l2g[lind, ei];
                # These should not be done yet.
                efind = efind+1;
                if invert
                    eforder[gind] = efind;
                else
                    eforder[efind] = gind;
                end
                node_done[gind] = true;
            end
        end

        # y=0 side
        for i = 1:porder+1
            lind = i;
            gind = l2g[lind, ei];
            if !node_done[gind]
                efind = efind+1;
                if invert
                    eforder[gind] = efind;
                else
                    eforder[efind] = gind;
                end
                node_done[gind] = true;
            end
        end

        # x=0 side
        for i = 2:porder+1
            lind = (porder+1)*(i-1) + 1;
            gind = l2g[lind, ei];
            if !node_done[gind]
                efind = efind+1;
                if invert
                    eforder[gind] = efind;
                else
                    eforder[efind] = gind;
                end
                node_done[gind] = true;
            end
        end

        # y=1 side
        for i = 2:porder+1
            lind = i + (porder+1)*porder;
            gind = l2g[lind, ei];
            if !node_done[gind]
                efind = efind+1;
                if invert
                    eforder[gind] = efind;
                else
                    eforder[efind] = gind;
                end
                node_done[gind] = true;
            end
        end

        # x=1 side
        for i = 2:porder
            lind = (porder+1)*i;
            gind = l2g[lind, ei];
            if !node_done[gind]
                efind = efind+1;
                if invert
                    eforder[gind] = efind;
                else
                    eforder[efind] = gind;
                end
                node_done[gind] = true;
            end
        end

    end

    return eforder;
end

function element_first_order_3d(grid, porder, elorder, invert = true)
    nnodes = size(grid.allnodes,2);
    eforder = zeros(Int, nnodes); # The ordering
    node_done = zeros(Bool, nnodes);

    # Loop over the elements in elorder, filling in the eforder as we go
    l2g = grid.loc2glb;
    efind = 0;
    for ei=elorder
        # Add interior first
        for k=2:porder
            for j = 2:porder
                for i = 2:porder
                    lind = i + (porder+1)*((j-1) + (porder+1)*(k-1));
                    gind = l2g[lind, ei];
                    # These should not be done yet.
                    efind = efind+1;
                    if invert
                        eforder[gind] = efind;
                    else
                        eforder[efind] = gind;
                    end
                    node_done[gind] = true;
                end
            end
        end

        # z=0 side
        for j = 1:porder+1
            for i = 1:porder+1
                lind = i + (porder+1)*(j-1);
                gind = l2g[lind, ei];
                if !node_done[gind]
                    efind = efind+1;
                    if invert
                        eforder[gind] = efind;
                    else
                        eforder[efind] = gind;
                    end
                    node_done[gind] = true;
                end
            end
        end

        # y=0 side
        for j = 2:porder+1
            for i = 1:porder+1
                lind = i + (porder+1)*(porder+1)*(j-1);
                gind = l2g[lind, ei];
                if !node_done[gind]
                    efind = efind+1;
                    if invert
                        eforder[gind] = efind;
                    else
                        eforder[efind] = gind;
                    end
                    node_done[gind] = true;
                end
            end
        end

        # x=0 side
        for j = 2:porder+1
            for i = 2:porder+1
                lind = (porder+1)*(i-1) + 1 + (porder+1)*(porder+1)*(j-1);
                gind = l2g[lind,ei];
                if !node_done[gind]
                    efind = efind+1;
                    if invert
                        eforder[gind] = efind;
                    else
                        eforder[efind] = gind;
                    end
                    node_done[gind] = true;
                end
            end
        end

        # z=1 side
        for j = 2:porder+1
            for i = 2:porder+1
                lind = i + (porder+1)*(j-1) + (porder+1)*(porder+1)*porder;
                gind = l2g[lind, ei];
                if !node_done[gind]
                    efind = efind+1;
                    if invert
                        eforder[gind] = efind;
                    else
                        eforder[efind] = gind;
                    end
                    node_done[gind] = true;
                end
            end
        end

        # y=1 side
        for j = 2:porder
            for i = 2:porder+1
                lind = i + (porder+1)*(porder+1)*(j-1) + (porder+1)*porder;
                gind = l2g[lind, ei];
                if !node_done[gind]
                    efind = efind+1;
                    if invert
                        eforder[gind] = efind;
                    else
                        eforder[efind] = gind;
                    end
                    node_done[gind] = true;
                end
            end
        end

        # x=1 side
        for j = 2:porder
            for i = 2:porder
                lind = (porder+1)*(i-1) + 1 + (porder+1)*(porder+1)*(j-1) + porder;
                gind = l2g[lind, ei];
                if !node_done[gind]
                    efind = efind+1;
                    if invert
                        eforder[gind] = efind;
                    else
                        eforder[efind] = gind;
                    end
                    node_done[gind] = true;
                end
            end
        end

    end

    return eforder;
end

function reorder_grid_element_first(grid, porder, elorder)
    dim = size(grid.allnodes,1);
    nel = size(grid.loc2glb,2);
    nfc = size(grid.face2glb,2);

    norder= get_element_first_order(dim, grid, porder, elorder, true);

    newnodes = zeros(size(grid.allnodes));
    newbdry = copy(grid.bdry);
    newloc2glb = copy(grid.loc2glb);
    newglbvertex = copy(grid.glbvertex);
    newface2glb = copy(grid.face2glb);
    newfaceVertex2glb = copy(grid.faceVertex2glb);

    for mi=1:length(norder)
        for d=1:dim
            newnodes[d,norder[mi]] = grid.allnodes[d,mi];
        end
    end

    for bid=1:length(newbdry)
        for i=1:length(newbdry[bid])
            newbdry[bid][i] = norder[grid.bdry[bid][i]];
        end
    end

    for ei=1:nel
        for ni=1:size(newloc2glb,1)
            newloc2glb[ni,ei] = norder[grid.loc2glb[ni,ei]];
        end
        for ni=1:size(newglbvertex,1)
            newglbvertex[ni,ei] = norder[grid.glbvertex[ni,ei]];
        end
    end

    for fi=1:nfc
        for ni=1:size(newface2glb,1)
            newface2glb[ni,fi] = norder[grid.face2glb[ni,fi]];
        end
        for ni=1:size(newfaceVertex2glb,1)
            newfaceVertex2glb[ni,fi] = norder[grid.faceVertex2glb[ni,fi]];
        end
    end

    return Femshop.Grid(newnodes, newbdry, grid.bdryface, grid.bdrynorm, grid.bids, newloc2glb, newglbvertex, newface2glb, newfaceVertex2glb);
end

function get_element_first_order(dim, grid, porder, elorder, invert = true)
    if dim == 3
        return element_first_order_3d(grid, porder, elorder, invert);
    elseif dim == 2
        return eforder_order_2d(grid, porder, elorder, invert);
    else
        return 1:size(grid.allnodes,2); # 1D doesn't make sense
    end
end


text:
using Distributions, Random
using LinearAlgebra
using DataFrames
using Plots
using RCall
using CSV
include("funcs.jl")

data = DataFrame(
        m20 = [3, 11, 10, 7, 4, 3, 2, 1, 0, 0, 0, 1, 1, 101],
        f20 = [0, 2, 4, 8, 9, 3, 0, 0, 0, 0, 0, 0, 0, 126],
        m32 = [7, 10, 11, 16, 3, 2, 1, 0, 0, 0, 0, 0, 0, 19],
        f32 = [1, 5, 11, 10, 5, 1, 0, 1, 0, 0, 0, 0, 0, 47],
        m50 = [5, 8, 11, 15, 4, 2, 1, 1, 0, 0, 0, 0, 0, 7],
        f50 = [0, 4, 6, 6, 3, 1, 1, 4, 0, 0, 0, 1, 1, 17],
        m80 = [4, 10, 8, 14, 8, 2, 1, 0, 0, 1, 0, 0, 0, 2],
        f80 = [2, 7, 15, 9, 3, 4, 1, 1, 0, 1, 0, 0, 0, 4]
)

β_p = Exponential(0.25)
α_p = Uniform(0,1)
γ_p = Gamma(1.5,1)
p_p = Gamma(1.5,1)

α_i = rand(α_p)
β_i = 0.00000001
γ_i = rand(γ_p)
p_i = rand(p_p)

# test of simple model
test = beetle_sim(144, 0.6, 5, 2, 2, 0.2)
println(test)
println(sum(test))

function test_prop()
    tot = zeros(14)
    for i = 1:10000
        tot += beetle_sim(144, 0.5, 5, 0.08, 0.31, 0.2)
    end
    return tot ./ 10000
end
# proportions from test of simple model
test_proportions = test_prop()

α_i = rand(α_p)
p_i = rand(p_p)
# test of extended model
ext_test = beetle_sim_ext(144, α_i, β_i, γ_i, p_i, 0.2)
println(ext_test)
println(sum(ext_test))

diggle = [3, 11, 10, 7, 4, 3, 2, 1, 0,0, 0, 1, 1, 101]

# abc simulation
sim = abc_sim(diggle, 0.2, 10^5, 144, α_p, β_p, γ_p, p_p)

max_stars_repo_name: aidan-cg/undergrad-thesis-ABC

max_issues_repo_name: aidan-cg/undergrad-thesis-ABC

max_forks_repo_name: aidan-cg/undergrad-thesis-ABC


text:
"""
    generate_reduce(g, gen_func, comp, reps; parallel=:threads)

Compute `gen_func(g)` `reps` times and return the instance `best` for which
`comp(best, v)` is true where `v` is all the other instances of `gen_func(g)`.

For example, `comp(x, y) = length(x) < length(y) ? x : y` then instance with the smallest
length will be returned.
"""
generate_reduce(g::AbstractGraph{T}, gen_func::Function, comp::Comp, reps::Integer; parallel=:threads) where {T<:Integer,Comp} =
parallel == :threads ? threaded_generate_reduce(g, gen_func, comp, reps) : distr_generate_reduce(g, gen_func, comp, reps)

"""
    distr_generate_min_set(g, gen_func, comp, reps)

Distributed implementation of [`LightGraphs.generate_reduce`](@ref).
"""
function distr_generate_reduce(
    g::AbstractGraph{T},
    gen_func::Function,
    comp::Comp,
    reps::Integer
    ) where {T<:Integer,Comp}
    # Type assert required for type stability
    min_set::Vector{T} = @distributed ((x, y)->comp(x, y) ? x : y) for _ in 1:reps
        gen_func(g)
    end
    return min_set
end

"""
    threaded_generate_reduce(g, gen_func, comp reps)

Multi-threaded implementation of [`LightGraphs.generate_reduce`](@ref).
"""
function threaded_generate_reduce(
    g::AbstractGraph{T},
    gen_func::Function,
    comp::Comp,
    reps::Integer
    ) where {T<:Integer,Comp}
    n_t = Base.Threads.nthreads()
    is_undef = ones(Bool, n_t)
    min_set = [Vector{T}() for _ in 1:n_t]
    Base.Threads.@threads for _ in 1:reps
        t = Base.Threads.threadid()
        next_set = gen_func(g)
        if is_undef[t] || comp(next_set, min_set[t])
            min_set[t] = next_set
            is_undef[t] = false
        end
    end

    min_ind = 0
    for i in filter((j)->!is_undef[j], 1:n_t)
        if min_ind == 0 || comp(min_set[i], min_set[min_ind])
            min_ind = i
        end
    end

    return min_set[min_ind]
end


text:
# Probabilistic Inference Generalises Boolean SATisfiability

using Sigma
using Base.Test

A = Sigma.flip()
B = Sigma.flip(0.9)
C = Sigma.flip()

a, b, c = rand((A, B, C), !A & !B & !C)
@test !a & !b & !c

a, b, c = rand((A, B, C), (A & B) | C)
@test (a & b) | c

a, b, c = rand((A, B, C), (A & B) | !C & A)
@test (a & b) | !c & a


text:
function GenEuropeanBinomialCall(K, T, S, sig, r, N)

        dt = T/N ## Time to Maturity/ Number of Steps
        nu = r - 0.5*sig^2
        dxu = sqrt(sig^2*dt + (nu*dt)^2)
        dxd = -dxu
        pu = 1/2 + 1/2*(nu*dt/dxu)
        pd = 1 - pu

        disc = exp(-r*dt)

        St = zeros(Float64, (N+1, N+1))

        # Initialising asset prices at maturity

        St[1 , N+1] = S / exp(N * dxd)
        for j in 2:(N+1)
                St[j, N+1] = St[j-1, N+1] / exp(dxu - dxd)
        end

        # create Array for options value tree

        C = zeros(Float64, (N+1, N+1))#linspace(1, N+1, N+1)

        # Initialise Option Values at maturity

        for j in 1:(N+1)
                C[j, N+1] = max(0.0, St[j, N+1] - K)
        end

        # Stepping back through the tree

        for i = linrange(N, 1, N)
                for j in 1:i
                      C[j, i] = disc * (pu * C[j, i+1] + pd * C[j+1, i+1])
                      #St[j,i] = St[j,i+1] / d
                      #P[j, i] = max(P[j, i], K - St[j,i])
                end
        end

C
end


text:
# http://juliastats.github.io/StatsBase.jl/stable/weights.html
mutable struct Packed_AliasingScalarSampler
  samples::Vector{Float64}
  weights::Vector{Float64}
end

## avert your eyes, C++ below
# class MvNormal : public Distribution {
#   std::vector<double> mean_;
#   std::vector<std::vector<double>> cov_;
#
# public:
#   MvNormal(const std::vector<double> &mean,
#            const std::vector<std::vector<double>> &cov)
#       : mean_(mean), cov_(cov) {}
#   json ToJson(void) const {
#     json j;
#     j["mean"] = mean_;
#     j["cov"] = cov_;
#     return (j);
#   }
# };


text:
module MyPackage
   # Declare dependencies
   using ForwardDiff, DelimitedFiles

   # Fuctions to export for external use
   export my_f, my_f_der

   # define your function here
   my_f(x)=3x^2-5x+1;

   # automated differentiation using ForwardDiff Package.
   my_f_der(t)=ForwardDiff.derivative(x->my_f(x),t)

end

max_stars_repo_name: mcodevb/MyPackage

max_issues_repo_name: mcodevb/MyPackage

max_forks_repo_name: mcodevb/MyPackage


text:
@turing_testset "is.jl" begin
    function reference(n)
        as = Vector{Float64}(undef, n)
        bs = Vector{Float64}(undef, n)
        logps = Vector{Float64}(undef, n)

        for i in 1:n
            as[i], bs[i], logps[i] = reference()
        end
        logevidence = logsumexp(logps) - log(n)

        return (as = as, bs = bs, logps = logps, logevidence = logevidence)
    end

    function reference()
        x = rand(Normal(4,5))
        y = rand(Normal(x,1))
        loglik = logpdf(Normal(x,2), 3) + logpdf(Normal(y,2), 1.5)
        return x, y, loglik
    end

    @model function normal()
        a ~ Normal(4,5)
        3 ~ Normal(a,2)
        b ~ Normal(a,1)
        1.5 ~ Normal(b,2)
        a, b
    end

    alg = IS()
    seed = 0
    n = 10

    model = normal()
    for i in 1:100
        Random.seed!(seed)
        ref = reference(n)

        Random.seed!(seed)
        chain = sample(model, alg, n)
        sampled = get(chain, [:a, :b, :lp])

        @test vec(sampled.a) == ref.as
        @test vec(sampled.b) == ref.bs
        @test vec(sampled.lp) == ref.logps
        @test chain.logevidence == ref.logevidence
    end

    @turing_testset "logevidence" begin
        Random.seed!(100)

        @model function test()
            a ~ Normal(0, 1)
            x ~ Bernoulli(1)
            b ~ Gamma(2, 3)
            1 ~ Bernoulli(x / 2)
            c ~ Beta()
            0 ~ Bernoulli(x / 2)
            x
        end

        chains = sample(test(), IS(), 10000)

        @test all(isone, chains[:x])
        @test chains.logevidence ≈ - 2 * log(2)
    end
end


text:
include("../src/params.jl")
using Test, .Params

@testset "Parameters" begin
    @test Params.bign >= Params.bigk >= 1 # No fewer investors than funds
    @test Params.bigt >= 1
    @test -1< Params.drift <1
    @test Params.marketvol >= 0
    @test 0 < Params.perfwindow[1] <= Params.perfwindow[end] < Params.bigt
    @test Params.marketstartval > 0
    # Positive starting value required given the multiplicative formula
    @test Params.marketvol >= 0
    @test Params.betamean > 0 # Average stock positively correlated with index
    @test Params.betastd >= 0
    @test 0 <= Params.stockvolrange[1] <= Params.stockvolrange[end]
    @test 0 <=Params.invcaprange[1] <= Params.invcaprange[2]
    @test Params.stockstartval >= 0
    @test -1 < Params.thresholdmean < 1
    @test Params.thresholdstd >= 0
    @test 1 <=Params.portfsizerange[1]<=Params.portfsizerange[end] <=Params.bigm
    @test 0 <= Params.impactrange[1] <= Params.impactrange[2]

end


text:
function ctr(f::Function, box=[0;1], N = 10)
    h = (box[2]-box[1])/N
    val = 0
    val += f(box[1])/2
    for j = 1:N-1
        x = box[1] + j*h
        val += f(x)
    end
    val += f(box[2])/2
    val *= h
    return val
end


text:
"""
    change_basis(basis,changeofbasis::AbstractMatrix)

# Examples

Compute the Lagrangian basis associated with a set of nodes

```jldoctests
using LinearAlgebra
using Gridap.Fields
using Gridap.Polynomials

D = 2
order = 1
f = MonomialBasis{D}(Float64,order)

nodes = Point{2,Int}[(0,0),(1,0),(0,1),(1,1)]
change = inv(evaluate(f,nodes))

g = change_basis(f,change)
println(evaluate(g,nodes))

# output
[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]

```
"""
function change_basis(basis,changeofbasis::AbstractMatrix)
  BasisFromChangeOfBasis(basis,changeofbasis)
end

function return_type(::typeof(change_basis),prebasis,matrix_inv)
  typeof(change_basis(prebasis,matrix_inv))
end

struct BasisFromChangeOfBasis{B,M} <: AbstractVector{Field}
  basis::B
  change::M
  function BasisFromChangeOfBasis(basis,change::AbstractMatrix)
    B = typeof(basis)
    M = typeof(change)
    new{B,M}(basis,change)
  end
end

struct BasisTermFromChangeOfBasis end

Base.size(a::BasisFromChangeOfBasis) = (length(a.basis),)
Base.axes(a::BasisFromChangeOfBasis) = (axes(a.basis,1),)
# @santiagobadia : Not sure we want to create the real computation here
Base.getindex(a::BasisFromChangeOfBasis,i::Integer) = BasisTermFromChangeOfBasis()
Base.IndexStyle(::BasisFromChangeOfBasis) = IndexLinear()

function return_cache(b::BasisFromChangeOfBasis,x)
  cb = return_cache(b.basis,x)
  bx = evaluate!(cb,b.basis,x)
  c = CachedArray(bx*b.change)
  (c,cb)
end

function evaluate!(cache,b::BasisFromChangeOfBasis,x)
  c, cb = cache
  bx = evaluate!(cb,b.basis,x)
  setsize!(c,size(bx))
  mul!(c.array,bx,b.change)
  c.array
end

function return_gradient_cache(b::BasisFromChangeOfBasis,x)
  cb = return_gradient_cache(b.basis,x)
  bx = evaluate_gradient!(cb,b.basis,x)
  c = CachedArray(bx*b.change)
  (c,cb)
end

function evaluate_gradient!(cache,b::BasisFromChangeOfBasis,x)
  c, cb = cache
  bx = evaluate_gradient!(cb,b.basis,x)
  setsize!(c,size(bx))
  mul!(c.array,bx,b.change)
  c.array
end

function return_hessian_cache(b::BasisFromChangeOfBasis,x)
  cb = return_hessian_cache(b.basis,x)
  bx = evaluate_hessian!(cb,b.basis,x)
  c = CachedArray(bx*b.change)
  (c,cb)
end

function evaluate_hessian!(cache,b::BasisFromChangeOfBasis,x)
  c, cb = cache
  bx = evaluate_hessian!(cb,b.basis,x)
  setsize!(c,size(bx))
  mul!(c.array,bx,b.change)
  c.array
end


text:
module Solc

    """
        Solc
    A Julia wrapper tool for the solc Solidity compiler.
    """

    using HTTP, JSON

    export
    getOSname,
    getInstallableSolcVersions,
    installSolc,

    solcCompiler

    include("install.jl")
    include("compiler.jl")

    function __init__()

        nothing
    end

end


text:
DiffEqBase.@def iipnlcachefields begin
  nlcache = alg.nlsolve.cache
  @unpack κ,tol,max_iter,min_iter,new_W = nlcache
  z = similar(u)
  dz = similar(u); tmp = similar(u); b = similar(u)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  uToltype = real(uBottomEltypeNoUnits)
  ηold = one(uToltype)

  nf = f isa SplitFunction && alg isa SplitAlgorithms ? f.f1 : f
  islin = (f isa ODEFunction && islinear(f.f)) || (f isa SplitFunction && islinear(f.f1.f))
  # check if `nf` is linear
  if islin && alg.nlsolve isa NLNewton
    # get the operator
    J = nf.f
    W = WOperator(f.mass_matrix, dt, J, true)
    du1 = rate_prototype
    uf = nothing
    jac_config = nothing
    linsolve = alg.linsolve(Val{:init},nf,u)
    z₊ = z
  elseif alg.nlsolve isa NLNewton
    if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype !== nothing
      W = WOperator(f, dt, true)
      J = nothing # is J = W.J better?
    else
      J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
      W = similar(J)
    end
    du1 = zero(rate_prototype)
    # if the algorithm specializes on split problems the use `nf`
    uf = DiffEqDiffTools.UJacobianWrapper(nf,t,p)
    jac_config = build_jac_config(alg,nf,uf,du1,uprev,u,tmp,dz)
    linsolve = alg.linsolve(Val{:init},uf,u)
    z₊ = z
  elseif typeof(alg.nlsolve) <: NLFunctional
    J = nothing
    W = nothing
    du1 = rate_prototype
    uf = nothing
    jac_config = nothing
    linsolve = nothing
    z₊ = similar(z)
  end

  if κ !== nothing
    κ = uToltype(nlcache.κ)
  else
    κ = uToltype(1//100)
  end
  if tol !== nothing
    tol = uToltype(nlcache.tol)
  else
    tol = uToltype(min(0.03,first(reltol)^(0.5)))
  end
end
DiffEqBase.@def oopnlcachefields begin
  nlcache = alg.nlsolve.cache
  @unpack κ,tol,max_iter,min_iter,new_W = nlcache
  z = uprev
  nf = f isa SplitFunction && alg isa SplitAlgorithms ? f.f1 : f
  if alg.nlsolve isa NLNewton
    # only use `nf` if the algorithm specializes on split eqs
    uf = DiffEqDiffTools.UDerivativeWrapper(nf,t,p)
  else
    uf = nothing
  end
  islin = (f isa ODEFunction && islinear(f.f)) || (f isa SplitFunction && islinear(f.f1.f))
  if (islin || DiffEqBase.has_jac(f)) && typeof(alg.nlsolve) <: NLNewton
    # get the operator
    J = islin ? nf.f : f.jac(uprev, p, t)
    if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
      J = DiffEqArrayOperator(J)
    end
    W = WOperator(f.mass_matrix, dt, J, false)
  else
    W = typeof(u) <: Number ? u : Matrix{uEltypeNoUnits}(undef, 0, 0) # uEltype?
  end
  uToltype = real(uBottomEltypeNoUnits)
  ηold = one(uToltype)

  if κ !== nothing
    κ = uToltype(nlcache.κ)
  else
    κ = uToltype(1//100)
  end
  if tol !== nothing
    tol = uToltype(nlcache.tol)
  else
    tol = uToltype(min(0.03,first(reltol)^(0.5)))
  end
  z₊,dz,tmp,b,k = z,z,z,z,rate_prototype
end

DiffEqBase.@def iipnlsolve begin
  _nlsolve = typeof(alg.nlsolve).name.wrapper
  nlcache = NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,γ,c,ηold,z₊,dz,tmp,b,k)
  nlsolve = _nlsolve{true, typeof(nlcache)}(nlcache)
end
DiffEqBase.@def oopnlsolve begin
  _nlsolve = typeof(alg.nlsolve).name.wrapper
  nlcache = NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,γ,c,ηold,z₊,dz,tmp,b,k)
  nlsolve = _nlsolve{false, typeof(nlcache)}(nlcache)
end

abstract type SDIRKMutableCache <: OrdinaryDiffEqMutableCache end

@cache mutable struct ImplicitEulerCache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
end

function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  atmp = similar(u,uEltypeNoUnits)
  γ, c = 1, 1
  @iipnlsolve
  ImplicitEulerCache(u,uprev,uprev2,du1,fsalfirst,k,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,nlsolve)
end

mutable struct ImplicitEulerConstantCache{F,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  γ, c = 1, 1
  @oopnlsolve
  ImplicitEulerConstantCache(uf,nlsolve)
end

mutable struct ImplicitMidpointConstantCache{F,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::ImplicitMidpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  γ, c = 1//2, 1//2
  @oopnlsolve
  ImplicitMidpointConstantCache(uf,nlsolve)
end

@cache mutable struct ImplicitMidpointCache{uType,rateType,JType,WType,UF,JC,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
end

function alg_cache(alg::ImplicitMidpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  γ, c = 1//2, 1//2
  @iipnlsolve
  ImplicitMidpointCache(u,uprev,du1,fsalfirst,k,z,dz,b,tmp,J,W,uf,jac_config,linsolve,nlsolve)
end

mutable struct TrapezoidConstantCache{F,uType,tType,N} <: OrdinaryDiffEqConstantCache
  uf::F
  uprev3::uType
  tprev2::tType
  nlsolve::N
end

function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  uprev3 = u
  tprev2 = t
  γ, c = 1//2, 1
  @oopnlsolve
  TrapezoidConstantCache(uf,uprev3,tprev2,nlsolve)
end

@cache mutable struct TrapezoidCache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,tType,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  uprev3::uType
  tprev2::tType
  nlsolve::N
end

function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  atmp = similar(u,uEltypeNoUnits)
  uprev3 = similar(u)
  tprev2 = t
  γ, c = 1//2, 1
  @iipnlsolve
  TrapezoidCache(u,uprev,uprev2,du1,fsalfirst,k,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,uprev3,tprev2,nlsolve)
end

mutable struct TRBDF2ConstantCache{F,Tab,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
  tab::Tab
end

function alg_cache(alg::TRBDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  tab = TRBDF2Tableau(uToltype,real(tTypeNoUnits))
  γ, c = tab.d, tab.γ
  @oopnlsolve
  TRBDF2ConstantCache(uf,nlsolve,tab)
end

@cache mutable struct TRBDF2Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,Tab,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  zprev::uType
  zᵧ::uType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  tab::Tab
end

function alg_cache(alg::TRBDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  atmp = similar(u,uEltypeNoUnits); zprev = similar(u);
  zᵧ = similar(u)
  tab = TRBDF2Tableau(uToltype,real(tTypeNoUnits))
  γ, c = tab.d, tab.γ
  @iipnlsolve
  TRBDF2Cache(u,uprev,du1,fsalfirst,k,zprev,zᵧ,z,dz,b,tmp,atmp,J,
              W,uf,jac_config,linsolve,nlsolve,tab)
end

mutable struct SDIRK2ConstantCache{F,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::SDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  γ, c = 1, 1
  @oopnlsolve
  SDIRK2ConstantCache(uf,nlsolve)
end

@cache mutable struct SDIRK2Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
end

function alg_cache(alg::SDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  z₁ = similar(u); z₂ = z;
  atmp = similar(u,uEltypeNoUnits)
  γ, c = 1, 1
  @iipnlsolve
  SDIRK2Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,dz,b,tmp,atmp,J,
              W,uf,jac_config,linsolve,nlsolve)
end

mutable struct SSPSDIRK2ConstantCache{F,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::SSPSDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  γ, c = 1//4, 1//1
  @oopnlsolve
  SSPSDIRK2ConstantCache(uf,nlsolve)
end

@cache mutable struct SSPSDIRK2Cache{uType,rateType,JType,WType,UF,JC,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  dz::uType
  b::uType
  tmp::uType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
end

function alg_cache(alg::SSPSDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  z₁ = similar(u); z₂ = z;
  atmp = similar(u,uEltypeNoUnits)
  γ, c = 1//4, 1//1
  @iipnlsolve
  SSPSDIRK2Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,dz,b,tmp,J,
                W,uf,jac_config,linsolve,nlsolve)
end

mutable struct Kvaerno3ConstantCache{UF,Tab,N} <: OrdinaryDiffEqConstantCache
  uf::UF
  nlsolve::N
  tab::Tab
end

function alg_cache(alg::Kvaerno3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  tab = Kvaerno3Tableau(uToltype,real(tTypeNoUnits))
  γ, c = tab.γ, 2tab.γ
  @oopnlsolve
  Kvaerno3ConstantCache(uf,nlsolve,tab)
end

@cache mutable struct Kvaerno3Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,Tab,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  tab::Tab
end

function alg_cache(alg::Kvaerno3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  z₁ = similar(u); z₂ = similar(u);
  z₃ = similar(u); z₄ = z;
  atmp = similar(u,uEltypeNoUnits)
  tab = Kvaerno3Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, 2tab.γ
  @iipnlsolve

  Kvaerno3Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,dz,b,tmp,atmp,J,
                W,uf,jac_config,linsolve,nlsolve,tab)
end

mutable struct Cash4ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
  tab::Tab
end

function alg_cache(alg::Cash4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  tab = Cash4Tableau(uToltype,real(tTypeNoUnits))
  γ, c = tab.γ,tab.γ
  @oopnlsolve
  Cash4ConstantCache(uf,nlsolve,tab)
end

@cache mutable struct Cash4Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,N,Tab,F} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  tab::Tab
end

function alg_cache(alg::Cash4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  z₁ = similar(u); z₂ = similar(u);
  z₃ = similar(u); z₄ = similar(u)
  z₅ = z
  atmp = similar(u,uEltypeNoUnits)
  tab = Cash4Tableau(uToltype,real(tTypeNoUnits))
  γ, c = tab.γ,tab.γ
  @iipnlsolve
  Cash4Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,dz,b,tmp,atmp,J,
              W,uf,jac_config,linsolve,nlsolve,tab)
end

mutable struct Hairer4ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
  tab::Tab
end

function alg_cache(alg::Union{Hairer4,Hairer42},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  if typeof(alg) <: Hairer4
    tab = Hairer4Tableau(uToltype,real(tTypeNoUnits))
  else
    tab = Hairer42Tableau(uToltype,real(tTypeNoUnits))
  end
  γ, c = tab.γ, tab.γ
  @oopnlsolve

  Hairer4ConstantCache(uf,nlsolve,tab)
end

@cache mutable struct Hairer4Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,Tab,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  tab::Tab
end

function alg_cache(alg::Union{Hairer4,Hairer42},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  z₁ = similar(u); z₂ = similar(u);
  z₃ = similar(u); z₄ = similar(u)
  z₅ = z
  dz = similar(u)
  atmp = similar(u,uEltypeNoUnits)

  if typeof(alg) <: Hairer4
    tab = Hairer4Tableau(uToltype,real(tTypeNoUnits))
  else # Hairer42
    tab = Hairer42Tableau(uToltype,real(tTypeNoUnits))
  end
  γ, c = tab.γ, tab.γ
  @iipnlsolve

  Hairer4Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,dz,b,tmp,atmp,J,
               W,uf,jac_config,linsolve,nlsolve,tab)
end


text:
using ..Math:
    AffineTransform

import ..Math:
    copy

const STATE_STACK_DEPTH = 100

# Stack state
mutable struct State
    clear_color::Palette
    draw_color::Palette
    current::AffineTransform

    function State()
        new(Gray(),
            Gray(),
            AffineTransform{Float64}())
    end

    function State(clear::Palette, draw::Palette, curr::AffineTransform)
        new(clear,
            draw,
            copy(curr))
    end
end

# copy specializations
copy(s::State) = State(s.clear_color, s.draw_color, s.current)

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using PlanningVsCompensatoryCaching
using Test


text:
#=
NOTE: the outliers are drawn separately (see the boxplot function)

    =====   y: whigh
      |
     +=+    y: q75
     | |
     | |
     +=+    y: q25
      |
    =====   y: wlow

x-axis:
    box:      p-bwidth/2  <>  p+bwidth/2
    whisker:  p-wwidth/2  <>  p+wwidth/2
=#

const boxplot_box_lstyle  = "set lstyle blstyle lwidth blwidth color blcolor\$"
const boxplot_med_lstyle  = "set lstyle medlstyle lwidth medlwidth color medcolor\$"
const boxplot_mean_mstyle = """
    gsave
    \tset color mmcol\$
    \tmarker mmarker\$ mmsize
    \tgrestore
    """
const boxplot_core_vertical = """
    ! -----------------------------------------
    \t! SET LINE STYLE THEN DRAW BOX AND WHISKERS
    \t! -----------------------------------------
    \tgsave
    \tset cap round
    \t$boxplot_box_lstyle
    \t!
    \t! >> LOWER WHISKER
    \tamove xg(p-wwidth/2) yg(wlow)
    \taline xg(p+wwidth/2) yg(wlow)
    \t!
    \t! >> CONNECTION LOWER WHISKER-BOX
    \tamove xg(p) yg(wlow)
    \t!
    \t! >> BOX
    \taline xg(p) yg(q25)
    \tamove xg(p-bwidth/2) yg(q25)
    \tbox xg(p+bwidth/2)-xg(p-bwidth/2) yg(q75)-yg(q25)
    \t!
    \t! >> CONNECTION BOX-UPPER WHISKER
    \tamove xg(p) yg(q75)
    \taline xg(p) yg(whigh)
    \t!
    \t! >> UPPER WHISKER
    \tamove xg(p-wwidth/2) yg(whigh)
    \taline xg(p+wwidth/2) yg(whigh)
    \tgrestore
    \t!
    \t! ------------------------------------
    \t! SET LINE STYLE THEN DRAW MEDIAN LINE
    \t! ------------------------------------
    \tgsave
    \t$boxplot_med_lstyle
    \tamove xg(p-bwidth/2) yg(q50)
    \taline xg(p+bwidth/2) yg(q50)
    \tgrestore
    \t!
    \t! --------------------------------------
    \t! SET MARKER STYLE THEN DRAW MEAN MARKER
    \t! --------------------------------------
    \tif (mshow > 0) then
    \t    amove xg(p) yg(mean)
    \t$boxplot_mean_mstyle
    \tend if"""
const boxplot_core_horizontal = """
    ! -----------------------------------------
    \t! SET LINE STYLE THEN DRAW BOX AND WHISKERS
    \t! -----------------------------------------
    \tgsave
    \tset cap round
    \t$boxplot_box_lstyle
    \tamove xg(wlow) yg(p-wwidth/2)
    \taline xg(wlow) yg(p+wwidth/2)
    \tamove xg(wlow) yg(p)
    \taline xg(q25) yg(p)
    \tamove xg(q25) yg(p-bwidth/2)
    \tbox xg(q75)-xg(q25) yg(p+bwidth/2)-yg(p-bwidth/2)
    \tamove xg(q75) yg(p)
    \taline xg(whigh) yg(p)
    \tamove xg(whigh) yg(p-wwidth/2)
    \taline xg(whigh) yg(p+wwidth/2)
    \tgrestore
    \t!
    \t! ------------------------------------
    \t! SET LINE STYLE THEN DRAW MEDIAN LINE
    \t! ------------------------------------
    \tgsave
    \t$boxplot_med_lstyle
    \tamove xg(q50) yg(p-bwidth/2)
    \taline xg(q50) yg(p+bwidth/2)
    \tgrestore
    \t!
    \t! --------------------------------------
    \t! SET MARKER STYLE THEN DRAW MEAN MARKER
    \t! --------------------------------------
    \tif (mshow > 0) then
    \t    amove xg(mean) yg(p)
    \t    $boxplot_mean_mstyle
    \tend if"""
const boxplot_args = ( # NOTE: leave space at end of string
    "p wlow q25 q50 q75 whigh mean ",
    "bwidth wwidth blstyle blwidth blcolor\$ ",
    "medlstyle medlwidth medcolor\$ ",
    "mshow mmarker\$ mmsize mmcol\$ "
)

GLE_DRAW_SUB["bp_vert"] = """
    sub bp_vert $(prod(boxplot_args))
        $boxplot_core_vertical
    end sub
    """

GLE_DRAW_SUB["bp_horiz"] = """
    sub bp_horiz $(prod(boxplot_args))
        $boxplot_core_horizontal
    end sub
    """


text:
@muladd function perform_step!(integrator,cache::SRA1ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  gpdt = integrator.g(uprev,p,t+dt)
  sqrt3 = sqrt(3one(eltype(W.dW)))
  chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  k₁ = dt*integrator.f(uprev,p,t)
  k₂ = dt*integrator.f(uprev+3k₁/4 + 3chi2.*gpdt/2,p,t+3dt/4)
  E₂ = chi2.*(integrator.g(uprev,p,t)-gpdt) #Only for additive!
  u = @.. uprev + k₁ / 3 + 2 * k₂ / 3 + E₂ + W.dW * gpdt

  if integrator.opts.adaptive
    E₁ = k₁ .+ k₂

    resids = calculate_residuals(E₁, E₂, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
  integrator.u = u
end

#=
@muladd function perform_step!(integrator,cache::SRA1Cache,f=integrator.f)
  @unpack chi2,tmp1,E₁,E₂,gt,k₁,k₂,gpdt,tmp = cache
  @unpack t,dt,uprev,u,W,p = integrator
  integrator.g(t,uprev,gt)
  integrator.g(t+dt,uprev,gpdt)
  integrator.f(t,uprev,k₁); k₁*=dt
  @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  @.. tmp1 = uprev+3k₁/4 + 3chi2*gpdt/2

  integrator.f(t+3dt/4,tmp1,k₂); k₂*=dt

  @.. E₂ = chi2*(gt-gpdt) #Only for additive!

  @.. u = uprev + k₁/3 + 2k₂/3 + E₂ + W.dW*gpdt

  if integrator.opts.adaptive
    @.. E₁ = k₁ + k₂

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
  integrator.u = u
end
=#

@muladd function perform_step!(integrator,cache::SRA1Cache,f=integrator.f)
  @unpack chi2,tmp1,E₁,E₂,gt,k₁,k₂,gpdt,tmp = cache
  @unpack t,dt,uprev,u,W,p = integrator
  integrator.g(gt,uprev,p,t)
  integrator.g(gpdt,uprev,p,t+dt)
  integrator.f(k₁,uprev,p,t); k₁*=dt
  sqrt3 = sqrt(3one(eltype(W.dW)))
  if typeof(W.dW) <: Union{SArray,Number}
    chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  else
    @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  end

  if is_diagonal_noise(integrator.sol.prob)
    @.. E₁ = chi2*gpdt
    @.. E₂ = chi2*(gt-gpdt) #Only for additive!
  else
    mul!(E₁,gpdt,chi2)
    @.. gt -= gpdt
    mul!(E₂,gt,chi2)
  end

  @.. tmp1 = uprev+3k₁/4 + 3E₁/2

  integrator.f(k₂,tmp1,p,t+3dt/4); @.. k₂*=dt

  if is_diagonal_noise(integrator.sol.prob)
    @.. tmp1 = W.dW*gpdt
  else
    mul!(tmp1,gpdt,W.dW)
  end

  @.. u = uprev + k₁/3 + 2k₂/3 + E₂ + tmp1

  if integrator.opts.adaptive
    @.. E₁ = k₁ + k₂

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
end

@muladd function perform_step!(integrator,cache::SRA2ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack a21,b21,c02,c11,c12,α1,α2,beta12,beta21,beta22 = cache
  sqrt3 = sqrt(3one(eltype(W.dW)))
  chi2 = .5*(W.dW + W.dZ/sqrt3) #I_(1,0)/h

  g1 = integrator.g(uprev,p,t+c11*dt)
  k1 = integrator.f(uprev,p,t)
  H01 = uprev + dt*a21*k1 + chi2*b21*g1

  g2 = integrator.g(H01,p,t+c12*dt)
  k2 = integrator.f(H01,p,t+c02*dt)

  E₁ = dt*(α1*k1 + α2*k2)
  E₂ = chi2*(beta21*g1 + beta22*g2)

  u = uprev + E₁ + E₂ + W.dW*(beta12*g2)

  if integrator.opts.adaptive
    E₁ = @.. dt * (k1 + k2)

    resids = calculate_residuals(E₁, E₂, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SRA2Cache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack chi2,tab,g1,g2,k1,k2,E₁,E₂,tmp = cache
  @unpack a21,b21,c02,c11,c12,α1,α2,beta12,beta21,beta22 = cache.tab
  H01 = E₁
  sqrt3 = sqrt(3one(eltype(W.dW)))
  if typeof(W.dW) <: Union{SArray,Number}
    chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  else
    @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  end

  integrator.g(g1,uprev,p,t+c11*dt)
  integrator.f(k1,uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    @.. H01 = uprev + dt*a21*k1 + chi2*b21*g1
  else
    mul!(E₁,g1,chi2)
    @.. H01 = uprev + dt*a21*k1 + b21*E₁
  end

  integrator.g(g2,H01,p,t+c12*dt)
  integrator.f(k2,H01,p,t+c02*dt)

  if is_diagonal_noise(integrator.sol.prob)
    @.. E₂ = chi2*(beta21*g1 + beta22*g2)
    @.. u = uprev + dt*(α1*k1 + α2*k2) + E₂ + W.dW*(beta12*g2)
  else
    @.. g1 = beta21*g1 + beta22*g2
    mul!(E₂,g1,chi2)
    g2 .*= beta12
    mul!(E₁,g2,W.dW)
    @.. u = uprev + dt*(α1*k1 + α2*k2) + E₂ + E₁
  end

  if integrator.opts.adaptive
    @.. E₁ = dt * (k1 + k2)

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
end

@muladd function perform_step!(integrator,cache::ThreeStageSRAConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack a21,a31,a32,b21,b31,b32,c02,c03,c11,c12,c13,α1,α2,α3,beta11,beta12,beta13,beta21,beta22,beta23 = cache
  sqrt3 = sqrt(3one(eltype(W.dW)))
  chi2 = .5*(W.dW .+ W.dZ/sqrt3) #I_(1,0)/h

  g1 = integrator.g(uprev,p,t+c11*dt)
  k1 = integrator.f(uprev,p,t)

  H01 = uprev + dt*a21*k1 + b21*chi2.*g1

  g2 = integrator.g(H01,p,t+c12*dt)
  k2 = integrator.f(H01,p,t+c02*dt)

  H02 = uprev + dt*(a31*k1 + a32*k2) + chi2.*(b31*g1 + b32*g2)

  g3 = integrator.g(H02,p,t+c13*dt)
  k3 = integrator.f(H02,p,t+c03*dt)

  E₁ = dt*(α1*k1 + α2*k2 + α3*k3)
  E₂ = chi2.*(beta21*g1 + beta22*g2 + beta23*g3)

  u = uprev + E₁ + E₂ + W.dW.*(beta11*g1 + beta12*g2 + beta13*g3)

  if typeof(integrator.alg) <: StochasticCompositeAlgorithm && typeof(integrator.alg.algs[1]) <: SOSRA2
    ϱu = integrator.opts.internalnorm(k3 - k2, t)
    ϱd = integrator.opts.internalnorm(H02 - H01, t)
    integrator.eigen_est = ϱu/ϱd
  end

  if integrator.opts.adaptive
    E₁ = @.. dt * (k1 + k2 + k3)

    resids = calculate_residuals(E₁, E₂, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::ThreeStageSRACache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack chi2,tab,g1,g2,g3,k1,k2,k3,E₁,E₂,tmp,gtmp = cache
  @unpack a21,a31,a32,b21,b31,b32,c02,c03,c11,c12,c13,α1,α2,α3,beta11,beta12,beta13,beta21,beta22,beta23 = cache.tab

  H01 = E₁; H02 = E₁
  sqrt3 = sqrt(3one(eltype(W.dW)))
  if typeof(W.dW) <: Union{SArray,Number}
    chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  else
    @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  end

  integrator.g(g1,uprev,p,t+c11*dt)
  integrator.f(k1,uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    @.. H01 = uprev + dt*a21*k1 + chi2*b21*g1
  else
    mul!(E₁,g1,chi2)
    @.. H01 = uprev + dt*a21*k1 + b21*E₁
  end

  integrator.g(g2,H01,p,t+c12*dt)
  integrator.f(k2,H01,p,t+c02*dt)

  if is_diagonal_noise(integrator.sol.prob)
    @.. H02 = uprev + dt*(a31*k1 + a32*k2) + chi2*(b31*g1 + b32*g2)
  else
    @.. gtmp = b31*g1 + b32*g2
    mul!(E₁,gtmp,chi2)
    @.. H02 = uprev + dt*(a31*k1 + a32*k2) + E₁
  end

  integrator.g(g3,H02,p,t+c13*dt)
  integrator.f(k3,H02,p,t+c03*dt)

  if is_diagonal_noise(integrator.sol.prob)
    @.. E₂ = chi2*(beta21*g1 + beta22*g2 + beta23*g3)
    @.. u = uprev + dt*(α1*k1 + α2*k2 + α3*k3) + E₂ + W.dW*(beta11*g1 + beta12*g2 + beta13*g3)
  else
    @.. gtmp = beta21*g1 + beta22*g2 + beta23*g3
    mul!(E₂,gtmp,chi2)
    @.. gtmp = beta11*g1 + beta12*g2 + beta13*g3
    mul!(E₁,gtmp,W.dW)
    u = uprev + dt*(α1*k1 + α2*k2 + α3*k3) + E₂ + E₁
  end

  if typeof(integrator.alg) <: StochasticCompositeAlgorithm && typeof(integrator.alg.algs[1]) <: SOSRA2
    @.. tmp = k3 - k2
    ϱu = integrator.opts.internalnorm(tmp, t)
    @.. tmp = H02 - H01
    ϱd = integrator.opts.internalnorm(tmp, t)
    integrator.eigen_est = ϱu/ϱd
  end

  if integrator.opts.adaptive
    @.. E₁ = dt * (k1 + k2 + k3)

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
  nothing
end

#=
@muladd function perform_step!(integrator,cache::SRACache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack H0,A0temp,B0temp,ftmp,gtmp,chi2,atemp,btemp,E₁,E₁temp,E₂,tmp = cache
  @unpack c₀,c₁,A₀,B₀,α,β₁,β₂,stages = cache.tab
  @.. chi2 = .5*(W.dW + W.dZ/sqrt3) #I_(1,0)/h
  for i in 1:stages
    fill!(H0[i],zero(eltype(integrator.u)))
  end
  for i = 1:stages
    fill!(A0temp,zero(eltype(integrator.u)))
    fill!(B0temp,zero(eltype(integrator.u)))
    for j = 1:i-1
      integrator.f((t + c₀[j]*dt),H0[j],ftmp)
      integrator.g((t + c₁[j]*dt),H0[j],gtmp)
      @.. A0temp = A0temp + A₀[j,i] * ftmp
      @.. B0temp = B0temp + B₀[j,i] * gtmp
    end
    @.. H0[i] = uprev + dt * A0temp + chi2 * B0temp
  end
  fill!(atemp ,zero(eltype(integrator.u)))
  fill!(btemp ,zero(eltype(integrator.u)))
  fill!(E₂    ,zero(eltype(integrator.u)))
  fill!(E₁temp,zero(eltype(integrator.u)))

  for i = 1:stages
    integrator.f((t+c₀[i]*dt),H0[i],ftmp)
    integrator.g((t+c₁[i]*dt),H0[i],gtmp)
    @.. atemp  =  atemp  + α[i] * ftmp
    @.. btemp  =  btemp  + (β₁[i] * W.dW) * gtmp
    @.. E₂     =  E₂     + (β₂[i] * chi2) * gtmp
    @.. E₁temp =  E₁temp +  ftmp
  end
  @.. u = uprev + dt * atemp + btemp + E₂

  if integrator.opts.adaptive
    @.. E₁ = dt*E₁temp

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
  integrator.u = u
end
=#

@muladd function perform_step!(integrator,cache::SRACache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack H0,A0temp,B0temp,ftmp,gtmp,chi2,atemp,btemp,E₁,E₁temp,E₂,tmp = cache
  @unpack c₀,c₁,A₀,B₀,α,β₁,β₂,stages = cache.tab
  sqrt3 = sqrt(3one(eltype(W.dW)))
  if typeof(W.dW) <: Union{SArray,Number}
    chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  else
    @.. chi2 = (W.dW + W.dZ/sqrt3)/2 #I_(1,0)/h
  end

  for i in 1:stages
    fill!(H0[i],zero(eltype(integrator.u)))
  end
  for i = 1:stages
    fill!(A0temp,zero(eltype(integrator.u)))
    fill!(B0temp,zero(eltype(integrator.u)))
    for j = 1:i-1
      integrator.f(ftmp,H0[j],p,t + c₀[j]*dt)
      integrator.g(gtmp,H0[j],p,t + c₁[j]*dt)
      @.. A0temp = A0temp + A₀[j,i] * ftmp
      if is_diagonal_noise(integrator.sol.prob)
        @.. B0temp = B0temp + B₀[j,i] * chi2 * gtmp
      else
        mul!(E₁temp,gtmp,chi2)
        @.. B0temp = B0temp + B₀[j,i] * E₁temp
      end
    end

    @.. H0[i] = uprev + dt * A0temp + B0temp
  end
  fill!(atemp ,zero(eltype(integrator.u)))
  fill!(btemp ,zero(eltype(integrator.u)))
  fill!(E₂    ,zero(eltype(integrator.u)))
  fill!(E₁temp,zero(eltype(integrator.u)))

  for i = 1:stages
    integrator.f(ftmp,H0[i],p,t+c₀[i]*dt)
    integrator.g(gtmp,H0[i],p,t+c₁[i]*dt)
    if is_diagonal_noise(integrator.sol.prob)
      @.. btemp = btemp + β₁[i]*W.dW*gtmp
    else
      mul!(E₁temp,gtmp,W.dW)
      @.. btemp = btemp + β₁[i]*E₁temp
    end
    if is_diagonal_noise(integrator.sol.prob)
      @.. E₂ = E₂ + β₂[i]*chi2*gtmp
    else
      mul!(E₁temp,gtmp,chi2)
      @.. E₂ = E₂ + β₂[i]*E₁temp
    end

    @.. atemp  =  atemp  + α[i]*ftmp
    @.. E₁temp =  E₁temp +  ftmp
  end

  @.. u = uprev + dt*atemp + btemp + E₂

  if integrator.opts.adaptive
    @.. E₁ = dt * E₁temp

    calculate_residuals!(tmp, E₁, E₂, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
end

@muladd function perform_step!(integrator,cache::SRAConstantCache,f=integrator.f)
  @unpack c₀,c₁,A₀,B₀,α,β₁,β₂,stages,H0 = cache
  @unpack t,dt,uprev,u,W,p = integrator
  sqrt3 = sqrt(3one(eltype(W.dW)))
  chi2 = .5*(W.dW + W.dZ/sqrt3) #I_(1,0)/h
  H0[:]=fill(zero(u),stages)

  for i = 1:stages
    A0temp = zero(u)
    B0temp = zero(u)
    for j = 1:i-1
      A0temp = A0temp .+ A₀[j,i] .* integrator.f(H0[j],p,t + c₀[j]*dt)
      B0temp = B0temp .+ B₀[j,i] .* integrator.g(H0[j],p,t + c₁[j]*dt) #H0[..,i] argument ignored
    end

    H0[i] = uprev + A0temp*dt + B0temp.*chi2
  end

  atemp = zero(u)
  btemp = zero(u)
  E₂    = zero(u)
  E₁temp= zero(u)

  for i = 1:stages
    ftmp = integrator.f(H0[i],p,t+c₀[i]*dt)
    E₁temp = E₁temp .+ ftmp
    atemp  = @.. atemp + α[i] * ftmp

    gtmp = integrator.g(H0[i], p, t + c₁[i] * dt) #H0[i] argument ignored
    btemp  = @.. btemp + (β₁[i] * W.dW) * gtmp
    E₂     = @.. E₂ + (β₂[i] * chi2) * gtmp
  end

  u = @.. uprev + dt * atemp + btemp + E₂

  if integrator.opts.adaptive
    E₁ = @.. dt * E₁temp

    resids = calculate_residuals(E₁, E₂, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end
  integrator.u = u
end


text:
#=

    This is a port of the R2 model LinearRegression.cs

    with the datafiles:
    - linear_regression_r2_x.csv
    - linear_regression_r2_y.csv
    - linear_regression_r2_x1.csv
    - linear_regression_r2_y1.csv

    * Using linear_regression_r2_x1.csv and linear_regression_r2_y1.csv:

  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

           a   -0.4114    0.0001     0.0000    0.0000   329.9865    1.0070       31.0751
           b    3.0589    0.0806     0.0025    0.0045   273.3509    1.0100       25.7417
    invNoise    1.2838    0.0280     0.0009    0.0012   341.9381    1.0049       32.2006

    * Using linear_regression_r2_x1.csv and linear_regression_r2_y1.csv:

  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

           a   -0.5363    0.0370     0.0012    0.0008   1937.8035    0.9998      793.8564
           b    3.0138    0.0427     0.0014    0.0009   1209.1662    0.9990      495.3569
    invNoise    1.4142    0.0314     0.0010    0.0007   1540.4202    0.9995      631.0611

=#

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

@model function linear_regression_r2(x,y)
    a ~ Normal(0,1)
    b ~ Normal(5.0,1.82574185835055371152)

    invNoise ~ Gamma(1,1)
    mu = a * x .+ b
    y ~ MvNormal(mu, invNoise)
end

x = parse.(Float64,split(readline("linear_regression_r2_x.csv"),","))
y = parse.(Float64,split(readline("linear_regression_r2_y.csv"),","))

# x = parse.(Float64,split(readline("linear_regression_r2_x1.csv"),","))
# y = parse.(Float64,split(readline("linear_regression_r2_y1.csv"),","))

model = linear_regression_r2(x,y)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 1_000)

chns = sample(model, NUTS(), 1_000)
# chns = sample(model, HMC(0.1,6), 1_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
check_libpath(PyCall) = begin
    if realpath(PyCall.libpython) == realpath(CONFIG.libpath)
        # @info "libpython path agrees between Python and PyCall" Python.CONFIG.libpath PyCall.libpython
    else
        @warn "Python and PyCall are using different versions of libpython. This will probably go badly." Python.CONFIG.libpath PyCall.libpython
    end
end

@init begin
    # Check if libpython is already loaded (i.e. if the Julia interpreter was started from a Python process)
    CONFIG.isembedded = haskey(ENV, "PYTHONJL_LIBPTR")

    if CONFIG.isembedded
        # In this case, getting a handle to libpython is easy
        CONFIG.libptr = Ptr{Cvoid}(parse(UInt, ENV["PYTHONJL_LIBPTR"]))
        # Check Python is initialized
        C.Py_IsInitialized() == 0 && error("Python is not already initialized.")
        CONFIG.isinitialized = CONFIG.preinitialized = true
    else
        # Find Python executable
        exepath = something(
            CONFIG.exepath,
            get(ENV, "PYTHONJL_EXE", nothing),
            Sys.which("python3"),
            Sys.which("python"),
            Some(nothing),
        )
        if exepath === nothing
            error(
                """
              Could not find Python executable.

              Ensure 'python3' or 'python' is in your PATH or set environment variable 'PYTHONJL_EXE'
              to the path to the Python executable.
              """,
            )
        end
        if exepath == "CONDA" || startswith(exepath, "CONDA:")
            CONFIG.isconda = true
            CONFIG.condaenv = exepath == "CONDA" ? Conda.ROOTENV : exepath[7:end]
            Conda._install_conda(CONFIG.condaenv)
            exepath = joinpath(
                Conda.python_dir(CONFIG.condaenv),
                Sys.iswindows() ? "python.exe" : "python",
            )
        end
        if isfile(exepath)
            CONFIG.exepath = exepath
        else
            error("""
                Python executable $(repr(exepath)) does not exist.

                Ensure either:
                - python3 or python is in your PATH
                - PYTHONJL_EXE is "CONDA" or "CONDA:<env>"
                - PYTHONJL_EXE is the path to the Python executable
                """)
        end

        # For calling Python with UTF-8 IO
        function python_cmd(args)
            env = copy(ENV)
            env["PYTHONIOENCODING"] = "UTF-8"
            setenv(`$(CONFIG.exepath) $args`, env)
        end

        # Find Python library
        libpath =
            something(CONFIG.libpath, get(ENV, "PYTHONJL_LIB", nothing), Some(nothing))
        if libpath !== nothing
            libptr = dlopen_e(path, CONFIG.dlopenflags)
            if libptr == C_NULL
                error("Python library $(repr(libpath)) could not be opened.")
            else
                CONFIG.libpath = libpath
                CONFIG.libptr = libptr
            end
        else
            for libpath in readlines(
                python_cmd([joinpath(@__DIR__, "find_libpython.py"), "--list-all"]),
            )
                libptr = dlopen_e(libpath, CONFIG.dlopenflags)
                if libptr == C_NULL
                    @warn "Python library $(repr(libpath)) could not be opened."
                else
                    CONFIG.libpath = libpath
                    CONFIG.libptr = libptr
                    break
                end
            end
            CONFIG.libpath === nothing && error("""
                Could not find Python library for Python executable $(repr(CONFIG.exepath)).

                If you know where the library is, set environment variable 'PYTHONJL_LIB' to its path.
                """)
        end

        # Compare libpath with PyCall
        PyCall = get(Base.loaded_modules, Base.PkgId(Base.UUID("438e738f-606a-5dbb-bf0a-cddfbfd45ab0"), "PyCall"), nothing)
        if PyCall === nothing
            @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" check_libpath(PyCall)
        else
            check_libpath(PyCall)
        end

        # Initialize
        with_gil() do
            if C.Py_IsInitialized() != 0
                # Already initialized (maybe you're using PyCall as well)
            else
                # Find ProgramName and PythonHome
                script = if Sys.iswindows()
                    """
                    import sys
                    print(sys.executable)
                    if hasattr(sys, "base_exec_prefix"):
                        sys.stdout.write(sys.base_exec_prefix)
                    else:
                        sys.stdout.write(sys.exec_prefix)
                    """
                else
                    """
                    import sys
                    print(sys.executable)
                    if hasattr(sys, "base_exec_prefix"):
                        sys.stdout.write(sys.base_prefix)
                        sys.stdout.write(":")
                        sys.stdout.write(sys.base_exec_prefix)
                    else:
                        sys.stdout.write(sys.prefix)
                        sys.stdout.write(":")
                        sys.stdout.write(sys.exec_prefix)
                    """
                end
                CONFIG.pyprogname, CONFIG.pyhome = readlines(python_cmd(["-c", script]))

                # Set PythonHome
                CONFIG.pyhome_w = Base.cconvert(Cwstring, CONFIG.pyhome)
                C.Py_SetPythonHome(pointer(CONFIG.pyhome_w))

                # Set ProgramName
                CONFIG.pyprogname_w = Base.cconvert(Cwstring, CONFIG.pyprogname)
                C.Py_SetProgramName(pointer(CONFIG.pyprogname_w))

                # Start the interpreter and register exit hooks
                C.Py_InitializeEx(0)
                atexit() do
                    CONFIG.isinitialized = false
                    CONFIG.version < v"3.6" ? C.Py_Finalize() : checkm1(C.Py_FinalizeEx())
                end
            end
            CONFIG.isinitialized = true
            check(
                C.Py_AtExit(
                    @cfunction(() -> (CONFIG.isinitialized = false; nothing), Cvoid, ())
                ),
            )
        end
    end

    C.PyObject_TryConvert_AddRules(
        "builtins.object",
        [(PyObject, CTryConvertRule_wrapref, -100), (PyRef, CTryConvertRule_wrapref, -200)],
    )
    C.PyObject_TryConvert_AddRules(
        "collections.abc.Sequence",
        [(PyList, CTryConvertRule_wrapref, 100)],
    )
    C.PyObject_TryConvert_AddRules(
        "collections.abc.Set",
        [(PySet, CTryConvertRule_wrapref, 100)],
    )
    C.PyObject_TryConvert_AddRules(
        "collections.abc.Mapping",
        [(PyDict, CTryConvertRule_wrapref, 100)],
    )
    C.PyObject_TryConvert_AddRules("_io._IOBase", [(PyIO, CTryConvertRule_trywrapref, 100)])
    C.PyObject_TryConvert_AddRules("io.IOBase", [(PyIO, CTryConvertRule_trywrapref, 100)])
    C.PyObject_TryConvert_AddRules(
        "<buffer>",
        [
            (PyArray, CTryConvertRule_trywrapref, 200),
            (PyBuffer, CTryConvertRule_wrapref, -200),
        ],
    )
    C.PyObject_TryConvert_AddRules(
        "<arrayinterface>",
        [(PyArray, CTryConvertRule_trywrapref, 200)],
    )
    C.PyObject_TryConvert_AddRules(
        "<arraystruct>",
        [(PyArray, CTryConvertRule_trywrapref, 200)],
    )
    C.PyObject_TryConvert_AddRules("<array>", [(PyArray, CTryConvertRule_trywrapref, 0)])

    with_gil() do

        @pyg `import sys, os`

        pywordsize = (@pyv `sys.maxsize > 2**32`::Bool) ? 64 : 32
        pywordsize == Sys.WORD_SIZE || error("Julia is $(Sys.WORD_SIZE)-bit but Python is $(pywordsize)-bit (at $(CONFIG.exepath ? "unknown location" : CONFIG.exepath))")

        if !CONFIG.isembedded
            @py ```
            # Some modules expect sys.argv to be set
            sys.argv = [""]
            sys.argv.extend($ARGS)

            # Some modules test for interactivity by checking if sys.ps1 exists
            if $(isinteractive()) and not hasattr(sys, "ps1"):
                sys.ps1 = ">>> "
            ```
        end

        # Is this the same Python as in Conda?
        if !CONFIG.isconda &&
           haskey(ENV, "CONDA_PREFIX") &&
           isdir(ENV["CONDA_PREFIX"]) &&
           haskey(ENV, "CONDA_PYTHON_EXE") &&
           isfile(ENV["CONDA_PYTHON_EXE"]) &&
           realpath(ENV["CONDA_PYTHON_EXE"]) == realpath(
               CONFIG.exepath === nothing ? @pyv(`sys.executable`::String) : CONFIG.exepath,
           )

            CONFIG.isconda = true
            CONFIG.condaenv = ENV["CONDA_PREFIX"]
            CONFIG.exepath === nothing && (CONFIG.exepath = @pyv(`sys.executable`::String))
        end

        # Get the python version
        CONFIG.version =
            let (a, b, c, d, e) = @pyv(`sys.version_info`::Tuple{Int,Int,Int,String,Int})
                VersionNumber(a, b, c, (d,), (e,))
            end
        v"3" ≤ CONFIG.version < v"4" || error(
            "Only Python 3 is supported, this is Python $(CONFIG.version) at $(CONFIG.exepath===nothing ? "unknown location" : CONFIG.exepath).",
        )

        # set up the 'juliaaa' module
        @py ```
        import sys
        if $(CONFIG.isembedded):
            jl = sys.modules["juliaaa"]
        elif "juliaaa" in sys.modules:
            raise ImportError("'juliaaa' module already exists")
        else:
            jl = sys.modules["juliaaa"] = type(sys)("juliaaa")
            jl.CONFIG = dict()
        jl.Main = $(pyjl(Main))
        jl.Base = $(pyjl(Base))
        jl.Core = $(pyjl(Core))
        code = """
        def newmodule(name):
            "A new module with the given name."
            return Base.Module(Base.Symbol(name))
        class As:
            "Interpret 'value' as type 'type' when converting to Julia."
            __slots__ = ("value", "type")
            def __init__(self, value, type):
                self.value = value
                self.type = type
            def __repr__(self):
                return "juliaaa.As({!r}, {!r})".format(self.value, self.type)
        """
        exec(code, jl.__dict__)
        ```

        # EXPERIMENTAL: hooks to perform actions when certain modules are loaded
        if !CONFIG.isembedded
            @py ```
            import sys
            class JuliaCompatHooks:
                def __init__(self):
                    self.hooks = {}
                def find_module(self, name, path=None):
                    hs = self.hooks.get(name)
                    if hs is not None:
                        for h in hs:
                            h()
                def add_hook(self, name, h):
                    if name not in self.hooks:
                        self.hooks[name] = [h]
                    else:
                        self.hooks[name].append(h)
                    if name in sys.modules:
                        h()
            JULIA_COMPAT_HOOKS = JuliaCompatHooks()
            sys.meta_path.insert(0, JULIA_COMPAT_HOOKS)

            # Before Qt is loaded, fix the path used to look up its plugins
            qtfix_hook = $(() -> (CONFIG.qtfix && fix_qt_plugin_path(); nothing))
            JULIA_COMPAT_HOOKS.add_hook("PyQt4", qtfix_hook)
            JULIA_COMPAT_HOOKS.add_hook("PyQt5", qtfix_hook)
            JULIA_COMPAT_HOOKS.add_hook("PySide", qtfix_hook)
            JULIA_COMPAT_HOOKS.add_hook("PySide2", qtfix_hook)
            ```

            @require IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" begin
                IJulia.push_postexecute_hook() do
                    CONFIG.pyplotautoshow && pyplotshow()
                end
            end
        end

        # EXPERIMENTAL: IPython integration
        if CONFIG.isembedded && CONFIG.ipythonintegration
            if !CONFIG.isipython
                @py ```
                try:
                    ok = "IPython" in sys.modules and sys.modules["IPython"].get_ipython() is not None
                except:
                    ok = False
                $(CONFIG.isipython::Bool) = ok
                ```
            end
            if CONFIG.isipython
                # Set `Base.stdout` to `sys.stdout` and ensure it is flushed after each execution
                @eval Base stdout = $(@pyv `sys.stdout`::PyIO)
                pushdisplay(TextDisplay(Base.stdout))
                pushdisplay(IPythonDisplay())
                @py ```
                mkcb = lambda cb: lambda: cb()
                sys.modules["IPython"].get_ipython().events.register("post_execute", mkcb($(() -> flush(Base.stdout))))
                ```
            end
        end
    end

    @debug "Initialized Python.jl" CONFIG.isembedded CONFIG.isinitialized CONFIG.exepath CONFIG.libpath CONFIG.libptr CONFIG.pyprogname CONFIG.pyhome CONFIG.version CONFIG.isconda CONFIG.condaenv
end


text:
module DSPopt

using Pkg
using Libdl
using SparseArrays
using StructJuMP
using JuMP # To reexport, should be using (not import)
using MPI
import MathOptInterface

const SJ = StructJuMP
const MOI = MathOptInterface

@SJ.exportall JuMP
export WassersteinSet

include("DSPCInterface.jl")

function __init__()
    try
        if Libdl.find_library("libDsp") == ""
            @warn("Could not load DSP shared library. Make sure it is in your library path.")
            global dspenv = DSPProblem(true)
        else
            Libdl.dlopen("libDsp", Libdl.RTLD_GLOBAL)
            global dspenv = DSPProblem(false)
        end
    catch
        @warn("Failed to initialize the package.")
        rethrow()
    end
end

function check_dsp()
    if dspenv.p == C_NULL
        @error("DSP is not initialied.")
    end
end

include("Core.jl")
include("Model.jl")
include("Dro.jl")
include("Parallel.jl")
include("JuMP_Wrapper.jl")
include("FileIO.jl")

end # module


text:
using HTTP
using JSON

function call(url)

    try
        response = HTTP.get(url)
        return JSON.parse(String(response.body))
    catch e
        return "ERROR: $e"
    end

end

function main()

    url = "https://cat-fact.herokuapp.com/facts/"

    r = call(url)

    println(r)
    println(keys(r[1]))
    println(r[1]["text"])

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
"""
    load_dcm_array(dcm_data::Vector{DICOM.DICOMData})
Given some DICOM.DICOMData, `load_dcm_array` loads the pixels
of each slice into a 3D array and returns the array
"""
function load_dcm_array(dcm_data::Vector{DICOM.DICOMData})
    return array = cat(
        [dcm_data[i][(0x7fe0, 0x0010)] for i in 1:length(dcm_data)]...; dims=3
    )
end

"""
    get_pixel_size(header)

Get the pixel information of the DICOM image given the `header` info.
Returns the x, y, and z values, where `z` corresponds to slice thickness

"""
function get_pixel_size(header)
	head = copy(header)
	pixel_size =
		try
			pixel_size = (head[(0x0028, 0x0030)])
            push!(pixel_size, head[(0x0018, 0x0050)])
		catch
			FOV = (head[(0x0018, 0x1100)])
			matrix_size = head[(0x0028, 0x0010)]

			pixel_size = FOV / matrix_size
            push!(pixel_size, head[(0x0018, 0x0050)])
		end
	return pixel_size
end


text:
using ApproxBayes
using Distributions
using DelimitedFiles
using Distances

# Run ABC assuming neutral model on copy number data only, with two summary statistics

# export JULIA_NUM_THREADS=32

# read from file
dataset=ARGS[1]
ncell=ARGS[2]
epsilon=ARGS[3]
birth_rate = ARGS[4]

# dataset="dataset_04022020"
# ncell=37
# epsilon=0.01
# birth_rate = 0.5

println(dataset)
println(ncell)

# change directory accordingly
dir = "../../"
progSim=dir * "/code/simpdo"
simdir=dir * "/data/abcsmc"
freal = dir * "/data/4trees/" * dataset * "_rvec.txt"

targetdata = readdlm(freal)[:, 1]
println(targetdata)

# Get the distance between simulated and real data (read from files) from copy numbers
function tumourABCneutralcn(params, constants, targetdata)
  # call simulation program to generate data
  # print(params)
  mu1=params[1]
  mu2=params[2]
  death_rate = 0

  ncell=constants[1]

  seed=rand(UInt32, 1)[1]
  suffix=string(seed) * "_" * dataset * "_n" * string(ncell)

  skip = 0

  cmdSim=`$progSim -o $dir -e $ncell --chr_prob $mu1 --arm_prob $mu2 --birth_rate $birth_rate --death_rate $death_rate --skip $skip --seed $seed --suffix $suffix`
  res = chomp(read(cmdSim, String))
  arr = split(res,"\t")
  simdata = map(x->parse(Float64,x), arr)
  simdata = simdata[1:2]

  # r = euclidean(targetdata, simdata)
  r = evaluate(Euclidean(1e-16), targetdata, simdata)

  return r, 1
end

setup = ABCSMC(tumourABCneutralcn, #simulation function
    2, # number of parameters
    parse(Float64, epsilon), # target ϵ
    Prior([Uniform(0, 1), Uniform(0, 1)]), #Prior for each of the parameters
    maxiterations = 9*10^6, #Maximum number of iterations before the algorithm terminates
    nparticles = 500,
    α = 0.3,
    ϵ1 = 10000.0,
    convergence = 0.05,
    constants = [ncell],
  )

smc = runabc(setup, targetdata, verbose = true, progress = true, parallel = false)

fout="smc-cna-" * dataset * "_n" * string(ncell) * "_epsilon" * string(epsilon)* "_brate" * string(birth_rate)
writeoutput(smc, dir=simdir, file=fout)