text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using LinearAlgebra
using Random
using Distributions
using StatsPlots

# %%
function simlogmlr(dist_true, n, logmaxlik0, logmaxlik1, L=10^5)
    sample = zeros(n)
    logmlr = zeros(L)
    for i in 1:L
        rand!(dist_true, sample)
        logmlr[i] = 2(logmaxlik1(sample) - logmaxlik0(sample))
    end
    logmlr
end

# %%
@time logmlr = simlogmlr(Normal(), 100,
    sample -> loglikelihood(Normal(), sample),
    sample -> loglikelihood(fit_mle(Normal, sample), sample),
)
histogram(logmlr; norm=true, alpha=0.3, bin=0:0.2:15, label="2log(maximum likelihood ratio)")
plot!(Chisq(2), 0, 15; label="Chisq(2)", lw=1.5)
title!("""
    dist_true = Normal(0, 1),  n = 100
    model0 = Normal(0, 1)
    model1 = Normal(μ, σ)
    """; titlefontsize=10)

# %%
@time logmlr = simlogmlr(Normal(), 100,
    sample -> loglikelihood(Normal(mean(sample), 1), sample),
    sample -> loglikelihood(fit_mle(Normal, sample), sample),
)
histogram(logmlr; norm=true, alpha=0.3, bin=0:0.1:6, label="2log(maximum likelihood ratio)")
plot!(Chisq(1), 0.05, 6; label="Chisq(1)", lw=1.5)
title!("""
    dist_true = Normal(0, 1),  n = 100
    model0 = Normal(μ, 1)
    model1 = Normal(μ, σ)
    """; titlefontsize=10)

# %%
@time logmlr = simlogmlr(Normal(), 100,
    sample -> loglikelihood(Normal(0.0, std(sample; corrected=false, mean=0.0)), sample),
    sample -> loglikelihood(fit_mle(Normal, sample), sample),
)
histogram(logmlr; norm=true, alpha=0.3, bin=0:0.1:6, label="2log(maximum likelihood ratio)")
plot!(Chisq(1), 0.05, 6; label="Chisq(1)", lw=1.5)
title!("""
    dist_true = Normal(0, 1),  n = 100
    model0 = Normal(0, σ)
    model1 = Normal(μ, σ)
    """; titlefontsize=10)

# %%
@time logmlr = simlogmlr(Normal(), 100,
    sample -> loglikelihood(Normal(), sample),
    sample -> loglikelihood(Normal(mean(sample), 1), sample),
)
histogram(logmlr; norm=true, alpha=0.3, bin=0:0.1:6, label="2log(maximum likelihood ratio)")
plot!(Chisq(1), 0.05, 6; label="Chisq(1)", lw=1.5)
title!("""
    dist_true = Normal(0, 1),  n = 100
    model0 = Normal(0, 1)
    model1 = Normal(μ, 1)
    """; titlefontsize=10)

# %%
@time logmlr = simlogmlr(Gamma(10, 1), 100,
    sample -> loglikelihood(Gamma(10, 1), sample),
    sample -> loglikelihood(fit_mle(Gamma, sample), sample),
)
histogram(logmlr; norm=true, alpha=0.3, bin=0:0.2:15, label="2log(maximum likelihood ratio)")
plot!(Chisq(2), 0, 15; label="Chisq(2)", lw=1.5)
title!("""
    dist_true = Gamma(10, 1),  n = 100
    model0 = Gamma(10, 1)
    model1 = Gamma(α, θ)
    """; titlefontsize=10)

# %%
n = 5
@time logmlr = simlogmlr(Normal(), n,
    sample -> loglikelihood(Normal(mean(sample), 1.0), sample),
    sample -> sum(x -> logpdf(Normal(x, 1.0), x), sample),
)
histogram(logmlr; norm=true, alpha=0.3, bin=range(0, 4n; length=61), label="2log(maximum likelihood ratio)")
plot!(Chisq(n-1), 0, 4n; label="Chisq(n-1)", lw=1.5)
title!("""
    dist_true = Normal(0, 1),  n = $n
    model0 = MvNormal(fill(μ, n), I)
    model1 = MvNormal(μ⃗, I)
    """; titlefontsize=10)

# %%
n = 5
@time logmlr = simlogmlr(Normal(), n,
    sample -> loglikelihood(Normal(0.0, std(sample; corrected=false, mean=0.0)), sample),
    sample -> loglikelihood(fit_mle(Normal, sample), sample),
)
histogram(logmlr; norm=true, alpha=0.3, bin=0:0.1:6, label="2log(maximum likelihood ratio)")
plot!(Chisq(1), 0.05, 6; label="Chisq(1)", lw=1.5)
title!("""
    dist_true = Normal(0, 1),  n = $n
    model0 = Normal(0, σ)
    model1 = Normal(μ, σ)
    """; titlefontsize=10)

# %%
n = 5
@time logmlr = simlogmlr(Normal(), n,
    sample -> loglikelihood(Normal(0.0, std(sample; corrected=false, mean=0.0)), sample),
    sample -> loglikelihood(fit_mle(Normal, sample), sample),
)
fstat = @. (n - 1)*(exp(logmlr/n) - 1)

P = histogram(fstat; norm=true, alpha=0.3, bin=range(0, 6; length=61), label="(t-statistic)²")
plot!(FDist(1, n-1), 0.05, 6; label="FDist(1, n-1)", lw=1.5)
plot!(Chisq(1), 0.05, 6; label="Chisq(1)", lw=1.5, ls=:dot)
title!("dist_true = Normal(0, 1),  n = $n"; titlefontsize=10)

Q = histogram(fstat; norm=true, alpha=0.3, bin=range(0, 30; length=61), label="(t-statistic)²")
plot!(FDist(1, n-1), 4, 30; label="FDist(1, n-1)", lw=1.5)
plot!(Chisq(1), 4, 30; label="Chisq(1)", lw=1.5, ls=:dash)
title!("dist_true = Normal(0, 1),  n = $n"; titlefontsize=10)
plot!(; xlim=(4, 30), ylim=(0, 0.034))

plot(P, Q; size=(800, 300))

# %%
function simfstat(dist_true, n, L=10^5)
    sample = zeros(n)
    fstat = zeros(L)
    μ₀ = mean(dist_true)
    σ₀² = var(dist_true)
    for i in 1:L
        rand!(dist_true, sample)
        fstat[i] = (mean(sample) - μ₀)^2/(var(sample)/n)
    end
    fstat
end

n = 5
@time fstat = simfstat(Normal(2, 1), n)

P = histogram(fstat; norm=true, alpha=0.3, bin=range(0, 6; length=61), label="(t-statistic)²")
plot!(FDist(1, n-1), 0.05, 6; label="FDist(1, n-1)", lw=1.5)
plot!(Chisq(1), 0.05, 6; label="Chisq(1)", lw=1.5, ls=:dot)
title!("dist_true = Normal(0, 1),  n = $n"; titlefontsize=10)

Q = histogram(fstat; norm=true, alpha=0.3, bin=range(0, 30; length=61), label="(t-statistic)²")
plot!(FDist(1, n-1), 4, 30; label="FDist(1, n-1)", lw=1.5)
plot!(Chisq(1), 4, 30; label="Chisq(1)", lw=1.5, ls=:dash)
title!("dist_true = Normal(0, 1),  n = $n"; titlefontsize=10)
plot!(; xlim=(4, 30), ylim=(0, 0.034))

plot(P, Q; size=(800, 300))

# %%

max_stars_repo_path: 0022/degree of freedom.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0022/degree of freedom.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0022/degree of freedom.jl

max_forks_repo_name: genkuroki/public


text:
### Load in with
### include(download("https://raw.githubusercontent.com/mth229/MTH229.jl/master/src/229.jl"))

using Pkg
pkg"add SymPy SpecialFunctions Roots ForwardDiff Plots QuadGK"

using Plots
using SpecialFunctions

import ForwardDiff
import QuadGK: quadgk

" f'(x) will find the derivative of `f` using Automatic Differentation from the `ForwardDiff` package "
Base.adjoint(f::Function) = x -> ForwardDiff.derivative(f, float(x))
D(f, n=1) = n > 1 ? D(D(f), n-1) : x -> ForwardDiff.derivative(f, float(x))

"""
Returns a function describing the tangent line to the graph of f at x=c.

Example. Where does the tangent line intersect the y axis?
```
f(x) = sin(x)
tl(x) = tangent(f, pi/4)(x)  # or tl = tangent(f, pi/3) to use a non-generic function
tl(0)
```

Uses the automatic derivative of `f` to find the slope of the tangent line at `x=c`.

"""
tangent(f,c) = x -> f(c) + f'(c) * (x-c)

"""
Returns a function describing the secant line to the graph of f at x=a and x=b.

Example. Where does the secant line intersect the y axis?
```
f(x) = sin(x)
a, b = pi/4, pi/3
sl(x) = secant(f, a, b)(x)  # or sl = sl(f, a, b) to use a non-generic function
sl(0)
```

"""
secant(f, a, b) = x -> f(a) + (f(b) - f(a)) / (b-a) * (x - a)

"""

`lim(f, c, n, dir="+")`: means to generate numeric table of values of `f` as `h` gets close to `c`.

Example:
```
f(x) = sin(x) / x
lim(f, 0)
```
"""
function lim(f::Function, c::Real; n::Int=6, dir="+")
	 hs = [(1/10)^i for i in 1:n] # close to 0
	 if dir == "+"
	   xs = c .+ hs
	 else
	   xs = c .- hs
	 end
	 ys = map(f, xs)
	 [xs ys]
end

"""

Simple implementation of the bisection method.

Example:

```julia
bisection(sin, 3, 4)
f(x) = x^5 - x^4 - x^3 - x^2 - x - 1
a = bisection(f, 1, 2)
f(a)
```

The display shows a simple graphic illustrating the method's division for the first few steps.

An easier-to-understand alternative to `Roots.find_zero(f, (a,b), Bisection())`.

"""
function bisection(f::Function, a, b)
    a,b = sort([a,b])

    if f(a) * f(b) > 0
        error("[a,b] is not a bracket. A bracket means f(a) and f(b) have different signs!")
    end

    M = a + (b-a) / 2

    i, j = 0, 64
    ss = fill("#", 65)
    ss[i+1]="a"; ss[j+1]="b"
    println("")
    println(join(ss))
    flag = true

    while a < M < b
        if flag && j-i == 1
            ss = fill(" ", 65)
            ss[j:(j+1)] = "⋮"
            println(join(ss))
            println("")
            flag = false
        end

        if f(M) == 0.0
            println("... exact answer found ...")
	    break
        end
        ## update step
	if f(a) * f(M) < 0
	    a, b = a, M

            if flag
                j = div(i + j, 2)
            end

	else
	    a, b = M, b

            if flag
                i = div(i + j, 2)
            end

	end

        if flag
            ss = fill(".", 65)
            ss[i+1]="a"; ss[j+1]="b"; ss[(i+2):j]="#"
            println(join(ss))
        end

        M = a + (b-a) / 2
    end
    M
end

import Roots
import Roots: newton, find_zero, find_zeros
newton(f, fp, x0; kwargs...) = Roots.find_zero((f,fp), x0, Roots.Newton(); kwargs...)
newton(f, x0; kwargs...) = newton(f, D(f), x0; kwargs...)
fzero(f, x0; kwargs...) = Roots.find_zero(f, x0; kwargs...)
fzero(f, a, b; kwargs...) = Roots.find_zero(f, (a, b); kwargs...)
fzeros(f, a, b; kwargs...) = Roots.find_zeros(f, a, b; kwargs...)

# some plotting utilities

"""
   trimplot(f, a, b, c=20; kwargs...)

Plot f over [a,b] but break graph if it exceeds c in absolute value.
"""
function trimplot(f, a, b, c=20; kwargs...)
  xs = linspace(a, b, 251) #range(a, stop=b, length=251)
  ys = f.(xs)

  us, vs = Real[], Real[]
  p = plot(us, vs, xlim=(a, b), legend=false, kwargs...)
  for (x,y) in zip(xs, ys)
    if abs(y) <= c
       push!(us, x); push!(vs, y)
    else
      length(us) > 0 && plot!(p, us, vs, color=:blue)
      empty!(us); empty!(vs)
    end
 end
 length(us) > 0 && plot!(p, us, vs, color=:blue)
 p
end

"""
    plotif(f, g, a, b)

Plot f colored depending on g < 0 or not.
"""
function plotif(f, g, a, b)
  plot([f, x -> g(x) > 0.0 ? f(x) : NaN], a, b, linewidth=5)
end

"""
   signchart(f, a, b)

Plot f over a,b with different color when negative.
"""
function signchart(f, a, b)
    p = plotif(f, f, a, b)
    plot!(p, zero)
    p
end

# visualize newtons method
function newton_vis(f, x0, a=Inf,b=-Inf; steps=5, kwargs...)
    xs = Float64[x0]
    for i in 1:steps
        push!(xs, xs[end] - f(xs[end]) / f'(xs[end]))
    end

    m,M = extrema(xs)
    m = min(m, a)
    M = max(M, b)

    p = plot(f, m, M; linewidth=3, legend=false, kwargs...)
    plot!(p, zero, m, M)
    for i in 1:steps
        plot!(p, [xs[i],xs[i],xs[i+1]], [0,f(xs[i]), 0])
        scatter!(p, xs[i:i],[0])
    end
    scatter!(p, [xs[steps+1]], [0])
    p
end

"""
riemann: compute Riemann sum approximations to a definite integral. As well, implement trapezoid and Simpson's rule.

Example:
```
f(x) = exp(x^2)
riemann(f, 0, 1, 1000)   # default right-Riemann sums
riemann(f, 0, 1, 1000, method="left")       # left sums
riemann(f, 0, 1, x1000, method="trapezoid")  # use trapezoid rule
riemann(f, 0, 1, 1000, method="simpsons")   # use Simpson's rule
```

"""
function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
  if method == "right"
     meth = (f,l,r) -> f(r) * (r-l)
  elseif method == "left"
     meth= (f,l,r) -> f(l) * (r-l)
  elseif method == "trapezoid"
     meth = (f,l,r) -> (1/2) * (f(l) + f(r)) * (r-l)
  elseif method == "simpsons"
     meth = (f,l,r) -> (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l)
  end

  xs = a + (0:n) * (b-a)/n
  as = [meth(f, l, r) for (l,r) in zip(xs[1:end-1], xs[2:end])]
  sum(as)
end

ext: jl


text:
for i in one(Int64):typemax(Int64)
    print(oct(i), " ")
    sleep(0.1)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestMLJUserReport

using MLJBase
using Test

import LightGBM

# Setup stuff
nrows = 500
nfeatures = 10

x = randn(nrows, nfeatures)
y = randn(nrows)

@testset "MLJ user report structure" begin

    # by settoing metric freq to 2 we can check the frequency is occurring correctly,
    # first time we expect only 2 metrics and then (metric @iter0 and metric @iter2)
    # if we go for 3 more iterations we expect 1 more metric
    estimator = LightGBM.LGBMRegression(;num_iterations=3, is_training_metric=true, metric=["l2"], metric_freq=2)
    metrics = LightGBM.fit!(estimator, x, y; verbosity=-1)

    report = LightGBM.MLJInterface.user_fitreport(estimator, metrics)

    # check types (loosely)
    @test report isa NamedTuple{(:training_metrics, :importance)}
    @test report.training_metrics isa Dict
    @test report.importance isa NamedTuple{(:gain, :split)}
    # check value properties
    @test length(report.importance.gain) == nfeatures
    @test length(report.importance.split) == nfeatures
    @test length(report.training_metrics["training"]["l2"]) == 2

    # continue for another 3 iterations, and check that the metrics lengths update properly
    new_metrics = LightGBM.train!(estimator, 3, String[], -1, LightGBM.Dates.now())

    new_report = LightGBM.MLJInterface.user_fitreport(estimator, metrics, new_metrics)

    # repeat all same tests except training metrics is now 3, not 1
    @test new_report isa NamedTuple{(:training_metrics, :importance)}
    @test new_report.training_metrics isa Dict
    @test new_report.importance isa NamedTuple{(:gain, :split)}
    # check value properties
    @test length(new_report.importance.gain) == nfeatures
    @test length(new_report.importance.split) == nfeatures
    @test length(new_report.training_metrics["training"]["l2"]) == 3

    # check metrics are right when freq is just 1
    new_estimator = LightGBM.LGBMRegression(;num_iterations=3, is_training_metric=true, metric=["l2"], metric_freq=1)
    freq_metrics = LightGBM.fit!(new_estimator, x, y; verbosity=-1)

    report = LightGBM.MLJInterface.user_fitreport(new_estimator, freq_metrics)
    @test length(report.training_metrics["training"]["l2"]) == 3

    new_freq_metrics = LightGBM.train!(new_estimator, 7, String[], -1, LightGBM.Dates.now())
    new_report = LightGBM.MLJInterface.user_fitreport(new_estimator, freq_metrics, new_freq_metrics)
    @test length(new_report.training_metrics["training"]["l2"]) == 10

end

end # module


text:
"""

    coordinatwisemedians(datamat)

    Return vector of medians of each column in a matrix.

# Arguments
- `datamat::Array{Float64, 2}`: A matrix.

# Example
```julia-repl
julia> mat = [1.0 2.0; 3.0 4.0; 5.0 6.0]
3×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
 5.0  6.0

julia> coordinatwisemedians(mat)
2-element Array{Float64,1}:
 3.0
 4.0
```
"""
function coordinatwisemedians(datamat::Array{Float64,2})::Array{Float64,1}
    n, p = size(datamat)
    meds = map(i -> median(datamat[:, i]), 1:p)
    return meds
end

"""

    bch(setting; alpha = 0.05, maxiter = 1000, epsilon = 0.000001)

Perform the Billor & Chatterjee & Hadi (2006) algorithm for the given regression setting.

# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
- `alpha::Float64`: Optional argument of the probability of rejecting the null hypothesis.
- `maxiter::Int`: Maximum number of iterations for calculating iterative weighted least squares estimates.
- `epsilon::Float64`: Accuracy for determining convergency.

# Examples
```julia-repl
julia> reg  = createRegressionSetting(@formula(calls ~ year), phones);
julia> Dict{Any,Any} with 7 entries:
"betas"                               => [-55.9205, 1.15572]
"squared.normalized.robust.distances" => [0.104671, 0.0865052, 0.0700692, 0.0553633, 0.0423875, 0.03…
"weights"                             => [0.00186158, 0.00952088, 0.0787321, 0.0787321, 0.0787321, 0…
"outliers"                            => [1, 14, 15, 16, 17, 18, 19, 20, 21]
"squared.normalized.residuals"        => [5.53742e-5, 2.42977e-5, 2.36066e-6, 2.77706e-6, 1.07985e-7…
"residuals"                           => [2.5348, 1.67908, 0.523367, 0.567651, 0.111936, -0.343779, …
"basic.subset"                        => [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  15, 16, 17, 18, 19, 20, …
```

# References
Billor, Nedret, Samprit Chatterjee, and Ali S. Hadi. "A re-weighted least squares method
for robust regression estimation." American journal of mathematical and management sciences 26.3-4 (2006): 229-252.
"""
function bch(setting::RegressionSetting; alpha=0.05, maxiter=1000, epsilon=0.000001)
    X = designMatrix(setting)
    y = responseVector(setting)
    return bch(X, y, alpha=alpha, maxiter=maxiter, epsilon=epsilon)
end

function bch(Xdesign::Array{Float64,2}, y::Array{Float64,1}; alpha=0.05, maxiter=1000, epsilon=0.000001)
    n, p = size(Xdesign)
    h = Int(floor((n + p + 1.0) / 2))
    X = Xdesign
    if Xdesign[:, 1] == ones(n)
        X = Xdesign[:, 2:end]
        n, p = size(X)
    end

    cnp = 1 + ((p + 1) / (n - p)) + (2 / (n - 1 - (3 * p)))
    chidist = Chisq(p)
    chicrit = quantile(chidist, 1 - alpha / n)
    crit = cnp * chicrit

    c = (2 * (n + 2 * p)) / (n - 2 * p)

    estimatedvariance = 0.0
    estimatedvariances = []

    # Algorithm 2 - Step 0.a
    coordmeds = coordinatwisemedians(X)
    A = ((X .- coordmeds')' * (X .- coordmeds')) / (n - 1)
    dsquared = diag(mahalanobisSquaredMatrix(DataFrame(X), meanvector=coordmeds, covmatrix=A))
    d = sqrt.(dsquared)

    # Algorithm 2 - Step 0.b
    bestindicesofd = sortperm(d)[1:h]
    colmeansofh = map(i -> mean(X[bestindicesofd, i]), 1:p)
    covmatofh = cov(X[bestindicesofd,:])
    newdsquared = diag(mahalanobisSquaredMatrix(DataFrame(X), meanvector=colmeansofh, covmatrix=covmatofh))
    newd = sqrt.(newdsquared)

    # Algorithm 2 - Steps 1, 2, 3
    basicsubsetindices = sortperm(newd)[1:(p + 1)]
    while length(basicsubsetindices) < h
        colmeanofbasicsubset = map(i -> mean(X[basicsubsetindices, i]), 1:p)
        covmatofbasicsubset = cov(X[basicsubsetindices,:])
        newdsquared = diag(mahalanobisSquaredMatrix(DataFrame(X), meanvector=colmeanofbasicsubset, covmatrix=covmatofbasicsubset))
        newd = sqrt.(newdsquared)
        basicsubsetindices = sortperm(newd)[1:(length(basicsubsetindices) + 1)]
    end

    # Algorithm 2 - Steps 4
    while length(basicsubsetindices) < n
        r = length(basicsubsetindices)
        colmeanofbasicsubset = map(i -> mean(X[basicsubsetindices, i]), 1:p)
        covmatofbasicsubset = cov(X[basicsubsetindices,:])
        newdsquared = diag(mahalanobisSquaredMatrix(DataFrame(X), meanvector=colmeanofbasicsubset, covmatrix=covmatofbasicsubset))
        newd = sqrt.(newdsquared)
        sortednewdsquared = sort(newdsquared)
        if sortednewdsquared[r + 1] >= crit
            break
        end
        basicsubsetindices = sortperm(newd)[1:(r + 1)]
    end

    # Algorithm 3 - Fitting
    squared_normalized_robust_distances = newd.^2.0 / sum(newd.^2.0)
    md = median(newd)
    newdmd = [newd[i] / maximum([newd[i], md]) for i in 1:n]
    newdmd2 = newdmd.^2.0
    sumnewdmd2 = sum(newdmd2)
    weights = newdmd2 / sumnewdmd2

    # Algorithm 3 - Step j
    betas = []
    squared_normalized_resids = []
    resids = []
    for i in 1:maxiter
        wols = wls(Xdesign, y, weights)
        betas = coef(wols)
        resids = residuals(wols)
        squared_normalized_resids = (resids.^2.0) / (sum(resids.^2.0))
        abssnresids = abs.(squared_normalized_resids)
        medsnresids = median(squared_normalized_resids)
        a = [(1 - weights[i]) / maximum([abssnresids[i], medsnresids]) for i in 1:n]
        weights = (a.^2.0) / (sum(a.^2.0))
        estimatedvariance = n * c * sum(resids.^2) / (n - p + 1)
        push!(estimatedvariances, estimatedvariance)
        if length(estimatedvariances) > 2
            if abs(estimatedvariance - estimatedvariances[end - 1]) < epsilon || abs(estimatedvariance - estimatedvariances[end - 2]) < epsilon
                break
            end
        end
    end

    # Report
    result = Dict()
    result["weights"] = weights
    result["betas"] = betas
    result["squared.normalized.residuals"] = squared_normalized_resids
    result["squared.normalized.robust.distances"] = squared_normalized_robust_distances
    result["residuals"] = resids
    result["outliers"] = filter(i -> abs(resids[i]) > 2.5, 1:n)
    result["basic.subset"] = sort(basicsubsetindices)
    return result
end

"""

    bchplot(setting::RegressionSetting; alpha=0.05, maxiter=1000, epsilon=0.00001)

Perform the Billor & Chatterjee & Hadi (2006) algorithm and generates outlier plot
for the given regression setting.

# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
- `alpha::Float64`: Optional argument of the probability of rejecting the null hypothesis.
- `maxiter::Int`: Maximum number of iterations for calculating iterative weighted least squares estimates.
- `epsilon::Float64`: Accuracy for determining convergency.

# References
Billor, Nedret, Samprit Chatterjee, and Ali S. Hadi. "A re-weighted least squares method
for robust regression estimation." American journal of mathematical and management sciences 26.3-4 (2006): 229-252.
"""
function bchplot(setting::RegressionSetting; alpha=0.05, maxiter=1000, epsilon=0.00001)
    X = designMatrix(setting)
    y = responseVector(setting)
    return bchplot(X, y, alpha=alpha, maxiter=maxiter, epsilon=epsilon)
end

function bchplot(Xdesign::Array{Float64,2}, y::Array{Float64,1}; alpha=0.05, maxiter=1000, epsilon=0.00001)
    result = bch(Xdesign, y, alpha=alpha, maxiter=maxiter, epsilon=epsilon)
    squared_normalized_residuals = result["squared.normalized.residuals"]
    squared_normalized_robust_distances = result["squared.normalized.robust.distances"]
    n = length(squared_normalized_robust_distances)
    scplot = scatter(squared_normalized_robust_distances,
            squared_normalized_residuals,
            legend=false,
            series_annotations=text.(1:n, :bottom),
            tickfont=font(10), guidefont=font(10), labelfont=font(10)
            )
    title!("Billor & Chatterjee & Hadi Plot")
    xlabel!("Squared Normalized Robust Distances")
    ylabel!("Squared Normalized Residuals")
end

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
module Models

using MLJModelInterface
import MLJBase # needed for UnivariateFinite in ConstantClassifier

const MMI = MLJModelInterface

include("DecisionTree.jl")

end


text:
using HomotopyContinuation, LinearAlgebra, DynamicPolynomials

# See https://www.juliahomotopycontinuation.org/examples/cyclooctane/
c² = 2
@polyvar z[1:3, 1:6]
z_vec = vec(z)[1:17] # the 17 variables in a vector
Z = [zeros(3) z[:,1:5] [z[1,6]; z[2,6]; 0] [√c²; 0; 0]] # the eight points in a matrix

# define the functions for cyclooctane:
F1 = [(Z[:, i] - Z[:, i+1]) ⋅ (Z[:, i] - Z[:, i+1]) - c² for i in 1:7]
F2 = [(Z[:, i] - Z[:, i+2]) ⋅ (Z[:, i] - Z[:, i+2]) - 8c²/3 for i in 1:6]
F3 = (Z[:, 7] - Z[:, 1]) ⋅ (Z[:, 7] - Z[:, 1]) - 8c²/3
F4 = (Z[:, 8] - Z[:, 2]) ⋅ (Z[:, 8] - Z[:, 2]) - 8c²/3
f = [F1; F2; F3; F4]

n = 2 # dimension of the cyclooctane variety
N = 17 # ambient dimension
@polyvar Aᵥ[1:n, 1:N] bᵥ[1:n] # variables for the linear equations
p = [vec(Aᵥ); bᵥ] # parameters
F = [f; Aᵥ * z_vec - bᵥ] # the polynomial system we have to solve

# now we solve one particular instance for A,b complex. we use this as start system
A₀ = randn(ComplexF64, n, N)
b₀ = randn(ComplexF64, n)
p₀ = [vec(A₀); b₀]

F₀ = [subs(Fᵢ, p => [vec(A₀); b₀]) for Fᵢ in F]
complex_result = solve(F₀)
println("Degree of the variety: ", nfinite(complex_result)) # should be 1408


text:
struct TriEdge
    ie::Int # Edge index in range [1,3]
    it::Int # Triangle index in range [1,nt]
end

struct TriMesh{X,Y,T}
    x::X # x-coordinates of the vertices, vector of length nv
    y::Y # y-coordinates of the verties, vector of length nv
    t::T # Vertex indices of the triangles, array of size (3,nt)
    neighbors::Array{Int,2} # Adjacency list of the mesh, array of size (3,nt)
    # Note: neighbors[ie,it] is index of triangle neighboring edge `ie` of triangle `it`
    bdr_loops::Vector{Vector{TriEdge}}
    bdr_map::Dict{TriEdge,Tuple{Int,Int}} # Map from edge index to boundary indices
    visited::Vector{Bool} # Keep track of visited elements for internal calculations
    bdr_visited::Vector{Vector{Bool}} # Keep track of visited boundary edges
    bdr_used::Vector{Bool} # Keep track of which boundary loops are used for filled contours
end

function TriMesh(x, y, t)
    nv = length(x)
    nt = size(t,2)

    # Compute adjacency of the mesh
    neighbors = zeros(Int,3,nt)
    edgemap = Dict{Tuple{Int,Int},TriEdge}()
    for it=1:nt
        for ie=1:3
            edge = (t[ie,it],t[mod1(ie+1,3),it])
            # Assume mesh is properly oriented, i.e. edges of neighboring triangles are
            # numbered in reversed order
            r_edge = (edge[2], edge[1])
            if haskey(edgemap,r_edge)
                tri_edge = edgemap[r_edge]
                neighbors[ie,it] = tri_edge.it
                neighbors[tri_edge.ie,tri_edge.it] = it
                delete!(edgemap,r_edge)
            else
                edgemap[edge] = TriEdge(ie,it)
            end
        end
    end
    empty!(edgemap)

    # Find all boundary edges, i.e. those with no neighboring triangles
    bdr_i = findall(neighbors .== 0)
    bdr_edges = Set{TriEdge}()
    for i=bdr_i
        push!(bdr_edges, TriEdge(i[1], i[2]))
    end

    # Follow each boundary edge to create boundary loops
    bdr_map = Dict{TriEdge,Tuple{Int,Int}}()
    bdr_loops = Vector{Vector{TriEdge}}()
    while !isempty(bdr_edges)
        bdr = Vector{TriEdge}()
        push!(bdr_loops, bdr)
        first_edge = first(bdr_edges.dict)[1] # same as first(bdr_edges), slightly faster
        edge = first_edge
        while true
            push!(bdr, edge)
            bdr_map[edge] = (length(bdr_loops), length(bdr))
            delete!(bdr_edges, edge)
            it = edge.it
            ie = mod1(edge.ie+1,3)
            point = t[ie,it]
            # Traverse neighbors until we find one with a boundary edge
            while neighbors[ie,it] != 0
                it = neighbors[ie,it]
                ie = findfirst(t[:,it] .== point)
                @assert !isnothing(ie)
            end
            edge = TriEdge(ie,it)
            edge == first_edge && break # Returned to first edge, completed loop
        end
    end

    visited = fill(false,2*nt)
    bdr_visited = Vector{Vector{Bool}}()
    for bdr=bdr_loops
        push!(bdr_visited, fill(false, length(bdr)))
    end
    bdr_used = fill(false, length(bdr_loops))

    TriMesh(x,y,t,neighbors,bdr_loops,bdr_map,visited,bdr_visited,bdr_used)
end

function get_neighbor_edge(m::TriMesh, edge::TriEdge)
    it = m.neighbors[edge.ie,edge.it]
    if it == 0
        return TriEdge(0,0)
    else
        point = m.t[mod1(edge.ie+1,3),edge.it]
        ie = findfirst(m.t[:,it] .== point)
        return TriEdge(ie,it)
    end
end

function reset_visited!(m::TriMesh)
    m.visited .= false
    for bv=m.bdr_visited
        bv .= false
    end
    m.bdr_used .= false
end


text:
using DiffEqBase: EnsembleProblem, ODEProblem

"""
$(SIGNATURES)

Compute the amplitude of the periodic orbit associated to `x`. The keyword argument `ratio = 1` is used as follows. If `length(x) = ratio * n`, the call returns the amplitude over `x[1:n]`.
"""
function getAmplitude(prob::AbstractShootingProblem, x::AbstractVector, p; ratio = 1)
	_max = _getExtremum(prob, x, p; ratio = ratio)
	_min = _getExtremum(prob, x, p; ratio = ratio, op = (min, minimum))
	return maximum(_max .- _min)
end

"""
$(SIGNATURES)

Compute the maximum of the periodic orbit associated to `x`. The keyword argument `ratio = 1` is used as follows. If `length(x) = ratio * n`, the call returns the amplitude over `x[1:n]`.
"""
function getMaximum(prob::AbstractShootingProblem, x::AbstractVector, p; ratio = 1)
	mx = _getExtremum(prob, x, p; ratio = ratio)
	return maximum(mx)
end
####################################################################################################
# Standard Shooting functional
"""
	pb = ShootingProblem(flow::Flow, ds, section; parallel = false)

Create a problem to implement the Standard Simple / Parallel Multiple Standard Shooting method to locate periodic orbits. More details (maths, notations, linear systems) can be found [here](https://rveltz.github.io/BifurcationKit.jl/dev/periodicOrbitShooting/). The arguments are as follows
- `flow::Flow`: implements the flow of the Cauchy problem though the structure [`Flow`](@ref).
- `ds`: vector of time differences for each shooting. Its length is written `M`. If `M==1`, then the simple shooting is implemented and the multiple one otherwise.
- `section`: implements a phase condition. The evaluation `section(x, T)` must return a scalar number where `x` is a guess for **one point** the periodic orbit and `T` is the period of the guess. The type of `x` depends on what is passed to the newton solver. See [`SectionSS`](@ref) for a type of section defined as a hyperplane.
- `parallel` whether the shooting are computed in parallel (threading). Available through the use of Flows defined by `EnsembleProblem`.

A functional, hereby called `G`, encodes the shooting problem. For example, the following methods are available:

- `pb(orbitguess, par)` evaluates the functional G on `orbitguess`
- `pb(orbitguess, par, du; δ = 1e-9)` evaluates the jacobian `dG(orbitguess).du` functional at `orbitguess` on `du`. The optional argument `δ` is used to compute a finite difference approximation of the derivative of the section.
- `pb`(Val(:JacobianMatrixInplace), J, x, par)` compute the jacobian of the functional analytically. This is based on ForwardDiff.jl. Useful mainly for ODEs.
- `pb(Val(:JacobianMatrix), x, par)` same as above but out-of-place.

You can then call `pb(orbitguess, par)` to apply the functional to a guess. Note that `orbitguess::AbstractVector` must be of size M * N + 1 where N is the number of unknowns of the state space and `orbitguess[M * N + 1]` is an estimate of the period `T` of the limit cycle. This form of guess is convenient for the use of the linear solvers in `IterativeSolvers.jl` (for example) which accepts only `AbstractVector`s. Another accepted guess is of the form `BorderedArray(guess, T)` where `guess[i]` is the state of the orbit at the `i`th time slice. This last form allows for non-vector state space which can be convenient for 2d problems for example, use `GMRESKrylovKit` for the linear solver in this case.

## Simplified constructors
- A simpler way to build the functional is to use
	pb = ShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, centers::AbstractVector; kwargs...)
where `F(x,p)` is the vector field, `p` is a parameter (to be passed to the vector field and the flow), `prob` is an `ODEProblem` (resp. `EnsembleProblem`) which is used to create a flow using the ODE solver `alg` (for example `Tsit5()`). `centers` is list of `M` points close to the periodic orbit, they will be used to build a constraint for the phase. `parallel = false` is an option to use Parallel simulations (Threading) to simulate the multiple trajectories in the case of multiple shooting. This is efficient when the trajectories are relatively long to compute. Finally, the arguments `kwargs` are passed to the ODE solver defining the flow. Look at `DifferentialEquations.jl` for more information. Note that, in this case, the derivative of the flow is computed internally using Finite Differences.

- Another way to create a Shooting problem with more options is the following where in particular, one can provide its own scalar constraint `section(x)::Number` for the phase
	pb = ShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, M::Int, section; parallel = false, kwargs...)
or

	pb = ShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, ds, section; parallel = false, kwargs...)
- The next way is an elaboration of the previous one
	pb = ShootingProblem(F, p, prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2, M::Int, section; parallel = false, kwargs...)
or

	pb = ShootingProblem(F, p, prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2, ds, section; parallel = false, kwargs...)
where we supply now two `ODEProblem`s. The first one `prob1`, is used to define the flow associated to `F` while the second one is a problem associated to the derivative of the flow. Hence, `prob2` must implement the following vector field ``\\tilde F(x,y,p) = (F(x,p),dF(x,p)\\cdot y)``.
"""
@with_kw struct ShootingProblem{Tf <: Flow, Ts, Tsection} <: AbstractShootingProblem
	M::Int64 = 0							# number of sections
	flow::Tf = Flow()						# should be a Flow
	ds::Ts = diff(LinRange(0, 1, M + 1))	# difference of times for multiple shooting
	section::Tsection = nothing				# sections for phase condition
	parallel::Bool = false					# whether we use DE in Ensemble mode for multiple shooting
end

# this constructor takes into accound a parameter passed to the vector field
# if M = 1, we disable parallel processing
function ShootingProblem(F, p, prob::ODEProblem, alg, ds, section; parallel = false, kwargs...)
	_M = length(ds)
	parallel = _M == 1 ? false : parallel
	_pb = parallel ? EnsembleProblem(prob) : prob
	return ShootingProblem(M = _M, flow = Flow(F, p, _pb, alg; kwargs...),
			ds = ds, section = section, parallel = parallel)
end

ShootingProblem(F, p, prob::ODEProblem, alg, M::Int, section; parallel = false, kwargs...) = ShootingProblem(F, p, prob, alg, diff(LinRange(0, 1, M + 1)), section; parallel = parallel, kwargs...)

ShootingProblem(F, p, prob::ODEProblem, alg, centers::AbstractVector; parallel = false, kwargs...) = ShootingProblem(F, p, prob, alg, diff(LinRange(0, 1, length(centers) + 1)), SectionSS(F(centers[1], p)./ norm(F(centers[1], p)), centers[1]); parallel = parallel, kwargs...)

# this is the "simplest" constructor to use in automatic branching from Hopf
ShootingProblem(M::Int, par, prob::ODEProblem, alg; parallel = false, kwargs...) = ShootingProblem(nothing, par, prob, alg, M, nothing; parallel = parallel, kwargs...)

ShootingProblem(M::Int, par, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; parallel = false, kwargs...) = ShootingProblem(nothing, par, prob1, alg1, prob2, alg2, M, nothing; parallel = parallel, kwargs...)

# idem but with an ODEproblem to define the derivative of the flow
function ShootingProblem(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2, ds, section; parallel = false, kwargs...)
	_M = length(ds)
	parallel = _M == 1 ? false : parallel
	_pb1 = parallel ? EnsembleProblem(prob1) : prob1
	_pb2 = parallel ? EnsembleProblem(prob2) : prob2
	return ShootingProblem(M = _M, flow = Flow(F, p, _pb1, alg1, _pb2, alg2; kwargs...), ds = ds, section = section, parallel = parallel)
end

ShootingProblem(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2, M::Int, section; parallel = false, kwargs...) = ShootingProblem(F, p, prob1, alg1, prob2, alg2, diff(LinRange(0, 1, M + 1)), section; parallel = parallel, kwargs...)

ShootingProblem(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2, centers::AbstractVector; parallel = false, kwargs...) = ShootingProblem(F, p, prob1, alg1, prob2, alg2, diff(LinRange(0, 1, length(centers) + 1)), SectionSS(F(centers[1], p)./ norm(F(centers[1], p)), centers[1]); parallel = parallel, kwargs...)

@inline isSimple(sh::ShootingProblem) = sh.M == 1
@inline isParallel(sh::ShootingProblem) = sh.parallel

# this function extracts the last component of the periodic orbit
@inline extractPeriodShooting(x::AbstractVector) = x[end]
@inline extractPeriodShooting(x::BorderedArray)  = x.p

# this function updates the section during the continuation run
function updateSection!(prob::ShootingProblem, x, par)
	xt = extractTimeSlices(x, getM(prob))
	@views update!(prob.section, prob.flow.F(xt[:, 1], par), xt[:, 1])
	prob.section.normal ./= norm(prob.section.normal)
	return true
end

"""
$(SIGNATURES)

Compute the period of the periodic orbit associated to `x`.
"""
@inline getPeriod(sh::ShootingProblem, x, par = nothing) = extractPeriodShooting(x)

function extractTimeSlices(x::AbstractVector, M::Int)
	N = div(length(x) - 1, M)
	return @views reshape(x[1:end-1], N, M)
end
extractTimeSlices(x::BorderedArray, M::Int) = x.u

@inline extractTimeSlice(x::AbstractMatrix, ii::Int) = @view x[:, ii]
@inline extractTimeSlice(x::AbstractVector, ii::Int) = xc[ii]

####################################################################################################
# Standard shooting functional using AbstractVector, convenient for IterativeSolvers.
function (sh::ShootingProblem)(x::AbstractVector, par)
	# Sundials does not like @views :(
	T = extractPeriodShooting(x)
	M = getM(sh)
	N = div(length(x) - 1, M)

	# extract the orbit guess and reshape it into a matrix as it's more convenient to handle it
	xc = extractTimeSlices(x, M)

	# variable to hold the computed result
	out = similar(x)
	outc = extractTimeSlices(out, M)

	if ~isParallel(sh)
		for ii in 1:M
			ip1 = (ii == M) ? 1 : ii+1
			# we can use views but Sundials will complain
			outc[:, ii] .= sh.flow(xc[:, ii], par, sh.ds[ii] * T) .- xc[:, ip1]
		end
	else
		solOde = sh.flow(xc, par, sh.ds .* T)
		for ii in 1:M
			ip1 = (ii == M) ? 1 : ii+1
			# we can use views but Sundials will complain
			outc[:, ii] .= @views solOde[ii][2] .- xc[:, ip1]
		end
	end

	# add constraint
	out[end] = @views sh.section(xc[:, 1], T)

	return out
end

# shooting functional, this allows for AbstractArray state space
function (sh::ShootingProblem)(x::BorderedArray, par)
	# period of the cycle
	T = extractPeriodShooting(x)
	M = getM(sh)

	# extract the orbit guess and reshape it into a matrix as it's more convenient to handle it
	xc = extractTimeSlices(x, M)

	# variable to hold the computed result
	out = similar(x)

	if ~isParallel(sh)
		for ii in 1:M
			# we can use views but Sundials will complain
			ip1 = (ii == M) ? 1 : ii+1
			out.u[ii] .= sh.flow(xc[ii], par, sh.ds[ii] * T) .- xc[ip1]
		end
	else
		@assert 1==0 "Not implemented yet. Try to use AbstractVectors instead"
	end

	# add constraint
	out.p = sh.section(xc[1], T)

	return out
end

# jacobian of the shooting functional
function (sh::ShootingProblem)(x::AbstractVector, par, dx::AbstractVector; δ = 1e-9)
	# period of the cycle
	# Sundials does not like @views :(
	dT = extractPeriodShooting(dx)
	T  = extractPeriodShooting(x)
	M = getM(sh)
	N = div(length(x) - 1, M)

	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)

	dxv = @view dx[1:end-1]
	dxc = reshape(dxv, N, M)

	# variable to hold the computed result
	out = similar(x)
	outv = @view out[1:end-1]
	outc = reshape(outv, N, M)

	if ~isParallel(sh)
		for ii in 1:M
			ip1 = (ii == M) ? 1 : ii+1
			# call jacobian of the flow
			tmp = sh.flow(xc[:, ii], par, dxc[:, ii], sh.ds[ii] * T)
			outc[:, ii] .= @views tmp.du .+ sh.flow.F(tmp.u, par) .* sh.ds[ii] .* dT .- dxc[:, ip1]
		end
	else
		# call jacobian of the flow
		solOde = sh.flow(xc, par, dxc, sh.ds .* T)
		for ii in 1:M
			ip1 = (ii == M) ? 1 : ii+1
			outc[:, ii] .= solOde[ii].du .+ sh.flow.F(solOde[ii].u, par) .* sh.ds[ii] .* dT .- dxc[:, ip1]
		end
	end

	# add constraint
	out[end] = @views (sh.section(x[1:N] .+ δ .* dx[1:N], T + δ * dT ) - sh.section(x[1:N], T)) / δ

	return out
end

# jacobian of the shooting functional, this allows for Array state space
function (sh::ShootingProblem)(x::BorderedArray, par, dx::BorderedArray; δ = 1e-9)
	dT = extractPeriodShooting(dx)
	T  = extractPeriodShooting(x)
	M = getM(sh)

	# variable to hold the computed result
	out = BorderedArray{typeof(x.u), typeof(x.p)}(similar(x.u), typeof(x.p)(0))

	if ~isParallel(sh)
		for ii in 1:M
			ip1 = (ii == M) ? 1 : ii+1
			# call jacobian of the flow
			tmp = sh.flow(x.u[ii], par, dx.u[ii], sh.ds[ii] * T)
			out.u[ii] .= tmp.du .+ sh.flow.F(tmp.u, par) .* sh.ds[ii] .* dT .- dx.u[ip1]
		end
	else
		@assert 1==0 "Not implemented yet. Try using AbstractVectors instead"
	end

	# add constraint
	x_tmp = similar(x.u); copyto!(x_tmp, x.u)
	axpy!(δ , dx.u, x_tmp)
	out.p = (sh.section(BorderedArray(x_tmp, T + δ * dT), T + δ * dT ) - sh.section(x, T)) / δ

	return out
end

# inplace computation of the matrix of the jacobian of the shooting problem, only serial for now
function (sh::ShootingProblem)(::Val{:JacobianMatrixInplace}, J::AbstractMatrix, x::AbstractVector, par)
	T = extractPeriodShooting(x)
	M = getM(sh)
	N = div(length(x) - 1, M)

	# extract the orbit guess and reshape it into a matrix as it's more convenient to handle it
	xc = extractTimeSlices(x, M)

	# jacobian of the flow
	dflow = (_J, _x, _T) -> ForwardDiff.jacobian!(_J, z -> sh.flow(Val(:SerialTimeSol), z, par, _T).u, _x)

	# put the matrices by blocks
	In = I(N)
	for ii=1:M
		@views dflow(J[(ii-1)*N+1:(ii-1)*N+N, (ii-1)*N+1:(ii-1)*N+N], xc[:, ii], sh.ds[ii] * T)
		# we put the identity matrices
		ip1 = (ii == M) ? 1 : ii+1
		if M == 1
			J[(ii-1)*N+1:(ii-1)*N+N, (ip1-1)*N+1:(ip1-1)*N+N] .+= -In
		else
			J[(ii-1)*N+1:(ii-1)*N+N, (ip1-1)*N+1:(ip1-1)*N+N] .= -In
		end
		# we fill the last column
		tmp = @views sh.flow(Val(:SerialTimeSol), xc[:, ii], par, sh.ds[ii] * T).u
		J[(ii-1)*N+1:(ii-1)*N+N, end] .= sh.flow.F(tmp, par) .* sh.ds[ii]
	end

	# we fill the last row
	@views ForwardDiff.gradient!(J[end, 1:N], z -> sh.section(z, T), x[1:N])
	J[end, end] = @views ForwardDiff.derivative(z -> sh.section(x[1:N], z), T)

	return J
end

# out of place version
(sh::ShootingProblem)(::Val{:JacobianMatrix}, x::AbstractVector, par) = sh(Val(:JacobianMatrixInplace), zeros(length(x), length(x)), x, par)
####################################################################################################

function _getExtremum(prob::ShootingProblem, x::AbstractVector, p; ratio = 1, op = (max, maximum))
	# this function extracts the amplitude of the cycle
	T = extractPeriodShooting(x)
	M = getM(prob)
	N = div(length(x) - 1, M)
	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)
	Th = eltype(x)
	n = div(N, ratio)

	# !!!! we could use @views but then Sundials will complain !!!
	if ~isParallel(prob)
		sol = prob.flow(Val(:Full), xc[:, 1], p, T)
		mx = @views op[2](sol[1:n, :], dims = 1)
	else # threaded version
		sol = prob.flow(Val(:Full), xc, p, prob.ds .* T)
		mx = op[2](sol[1][1:n, :] , dims = 2)
		for ii = 2:M
			mx = op[1].(mx, op[2](sol[ii][1:n, :], dims = 2))
		end
	end
	return mx
end

"""
$(SIGNATURES)

Compute the full trajectory associated to `x`. Mainly for plotting purposes.
"""
function getTrajectory(prob::ShootingProblem, x::AbstractVector, p)
	T = extractPeriodShooting(x)
	M = getM(prob)
	N = div(length(x) - 1, M)
	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)
	Th = eltype(x)

	# !!!! we could use @views but then Sundials will complain !!!
	if ~isParallel(prob)
		return prob.flow(Val(:Full), xc[:, 1], p, T)
	else # threaded version
		sol = prob.flow(Val(:Full), xc[:, 1:1], p, [T])
		return sol[1]
	end
end

####################################################################################################
# functions needed for Branch switching from Hopf bifurcation point
function problemForBS(prob::ShootingProblem, F, dF, par, hopfpt, ζr, orbitguess_a, period)
	# append period at the end of the initial guess
	orbitguess_v = reduce(vcat, orbitguess_a)
	orbitguess = vcat(vec(orbitguess_v), period) |> vec

	# update the problem but not the section if the user passed one
	probSh = setproperties(prob, section = isnothing(prob.section) ? SectionSS(F(orbitguess_a[1], hopfpt.params), copy(orbitguess_a[1])) : prob.section)
	probSh.section.normal ./= norm(probSh.section.normal)

	# be sure that the vector field is correctly inplace in the Flow structure
	probSh = @set probSh.flow.F = F

	return probSh, 	orbitguess
end


text:
using ComplexMixtures
using PDBTools
using Plots
using LaTeXStrings
using EasyFit

# Plot defaults
plot_font = "Computer Modern"
default(
    fontfamily=plot_font,
    linewidth=2.5,
    framestyle=:box,
    label=nothing,
    grid=false,
    palette=:tab10
)
scalefontsizes(); scalefontsizes(1.3)

# Load system PDB file
system = readPDB("./simulation/equilibrated.pdb")

# Select the atoms corresponding to glycerol and water
dmf = select(system,"resname DMF")
acr = select(system,"resname FACR or resname ACR or resname LACR")

# Set ComplexMixtures selections
solute = Selection(acr,nmols=1)
solvent = Selection(dmf,natomspermol=12)

# Load previously computed data
mddf_dmf_acr = load("./results/mddf_dmf_acr.json")

# Plot the MDDF of DMF relative to PolyACR and its corresponding KB integral
plot(layout=(2,1))
x = mddf_dmf_acr.d
y = movavg(mddf_dmf_acr.mddf,n=10).x
plot!(x,y,ylabel="MDDF",subplot=1)
plot!(xlim=(0,20),subplot=1)

y = movavg(mddf_dmf_acr.kb,n=10).x
plot!(x,y,xlabel=L"\textrm{Distance / \AA}",ylabel=L"\textrm{KB~/~cm^3~mol^{-1}}",subplot=2)
plot!(xlim=(0,20),subplot=2)
savefig("./results/mddf_kb.png")

# Plot DMF group contributions to the MDDF. We use a vector of tuples here,
# the first element of the tuple is a list of atom names, and the second element
# is the label of the group, to be used in the plot.
groups = [
    (["C","O"],"CO"), # carbonyl
    (["HA"],"HA"),
    (["N"],"N"),
    (["CC","CT","HC1","HC2","HC3","HT1","HT2","HT3"],"Methyl groups")
]
plot(layout=(2,1))
x = mddf_dmf_acr.d
y = movavg(mddf_dmf_acr.mddf,n=10).x
plot!(x,y,label="Total",subplot=1)
for group in groups
    # Retrieve the contributions of the atoms of this group
    group_contrib = contrib(solvent,mddf_dmf_acr.solvent_atom,group[1])
    y = movavg(group_contrib,n=10).x
    plot!(x,y,label=group[2],subplot=1)
end
plot!(
    xlim=(1,5),
    ylabel="MDDF",
    subplot=1
)

# Plot ACR group contributions to the MDDF
groups = [
    (["CF","HF1","HF2","HF3",
      "CL","HL1","HL2","HL3"],L"\textrm{CH_3}"), # terminal methyles
    (["OE1","CD"],"CO"), # carbonyl
    (["NE2","HE22","HE21"],L"\textrm{NH_2}"), # amine
    (["C","H2","H1","CA","HA"],L"\textrm{CHCH_2}"), # backbone
]
x = mddf_dmf_acr.d
y = movavg(mddf_dmf_acr.mddf,n=10).x
plot!(x,y,label="Total",subplot=2)
for group in groups
    group_contrib = contrib(solute,mddf_dmf_acr.solute_atom,group[1])
    y = movavg(group_contrib,n=10).x
    plot!(x,y,label=group[2],subplot=2)
end
plot!(
    xlim=(1,5),
    xlabel=L"\textrm{Distance / \AA}",
    ylabel="MDDF", subplot=2
)
savefig("./results/mddf_groups.png")

# 2D plot of group contributions
groups = [
    (["CF","HF1","HF2","HF3"],L"\textrm{CH_3}"),
    (["OE1","CD"],"CO"),
    (["NE2","HE22","HE21"],L"\textrm{NH_2}"),
    (["C","H2","H1","CA","HA"],L"\textrm{CHCH_2}"),
    (["CL","HL1","HL2","HL3"],L"\textrm{CH_3}")
]

# Here we split the polymer in residues, to extract the contribution of
# each chemical group of each polymer mer independently
mers = collect(eachresidue(acr))
group_contribs = zeros(length(mddf_dmf_acr.d),0)
labels = String[]
for (imer, mer) in pairs(mers)
    # The mer.atoms array is just the original array of atoms,
    # and the atoms of this residue are retrieved by choosing the
    # corresponding range (the eachresidue function does not copy)
    mer_atoms = mer.atoms[mer.range]
    for igroup in 1:length(groups)
        group = groups[igroup]
        (imer != 1 && igroup == 1) && continue # only first residue has a first CH3
        (imer != 5 && igroup == 5) && continue # only last residue has a terminal CH3
        # And from the mer_atoms atoms, filter the ones corresponding to this group
        atoms = filter(at -> at.name in group[1], mer_atoms)
        # Retrive the contribution of these atoms to the MDDF
        contribs = movavg(contrib(solute,mddf_dmf_acr.solute_atom,atoms),n=10).x
        # Concatenate the results to build the 2D matrix
        group_contribs = hcat(group_contribs,contribs)
        # Push label to label list
        push!(labels,group[2])
    end
end

# Find the indices of the limits of the map we want
idmin = findfirst( d -> d > 1.5, mddf_dmf_acr.d)
idmax = findfirst( d -> d > 3.2, mddf_dmf_acr.d)

# Plot contour map
contourf(
    1:length(labels),
    mddf_dmf_acr.d[idmin:idmax],
    group_contribs[idmin:idmax,:],
    color=cgrad(:tempo),linewidth=1,linecolor=:black,
    colorbar=:none,levels=10,
    xlabel="Group",ylabel=L"r/\AA",xrotation=60,
    xticks=(1:length(labels),labels),
    margin=5Plots.Measures.mm # adjust margin
)
savefig("./results/map2D_acr.png")


text:
function canFgBg(effect::Maple.Effect)
    return true, true
end

function editingOptions(effect::Maple.Effect)
    return Dict{String, Any}()
end

function registerPlacements!(placements::Array{Type{Maple.Effect{T}} where T, 1}, loaded::Array{String, 1})
    empty!(placements)

    for modul in loaded
        if hasModuleField(modul, "placements")
            res = getModuleField(modul, "placements")

            if isa(res, Array)
                append!(placements, res)

            else
                push!(placements, res)
            end
        end
    end
end

const loadedEffects = joinpath.(abs"effects", readdir(abs"effects"))
loadModule.(loadedEffects)
const effectPlacements = Array{Type{Maple.Effect{T}} where T, 1}()


text:
#using .Bio3DView

"""
viewstruct

based on Bio3Dview.jl
"""
function viewstruc(ta::TrjArray; kwargs...)
    io = IOBuffer()
    writepdb(io, ta)
    return Bio3DView.view("data-type='pdb'", String(take!(io)); style=Bio3DView.defaultstyle("pdb"), kwargs...)
end

#cpu(m) = fmap(x -> adapt(Array, x), m)
#gpu(x) = use_cuda[] ? fmap(CuArrays.cu, x) : x

function gpu(ta::TrjArray{T, U}) where {T, U}
    TrjArray{T, U}(CuArray(ta.x), CuArray(ta.y), CuArray(ta.z), CuArray(ta.boxsize),
      ta.chainname, CuArray(ta.chainid),
      ta.resname, CuArray(ta.resid),
      ta.atomname, CuArray(ta.atomid),
      CuArray(ta.mass), CuArray(ta.charge),
      CuArray(ta.list_bond), CuArray(ta.list_angle), CuArray(ta.list_dihedral),
      CuArray(ta.list_improper), CuArray(ta.list_cmap))
end

function gpu32(ta::TrjArray)
    TrjArray{Float32, Int64}(CuArrays.cu(ta.x), CuArrays.cu(ta.y), CuArrays.cu(ta.z), CuArrays.cu(ta.boxsize),
      ta.chainname, CuArrays.cu(ta.chainid),
      ta.resname, CuArrays.cu(ta.resid),
      ta.atomname, CuArrays.cu(ta.atomid),
      CuArrays.cu(ta.mass), CuArrays.cu(ta.charge),
      CuArrays.cu(ta.list_bond), CuArrays.cu(ta.list_angle), CuArrays.cu(ta.list_dihedral),
      CuArrays.cu(ta.list_improper), CuArrays.cu(ta.list_cmap))
end

function cpu(ta::TrjArray{T, U}) where {T, U}
    TrjArray{T, U}(Array(ta.x), Array(ta.y), Array(ta.z), Array(ta.boxsize),
      ta.chainname, Array(ta.chainid),
      ta.resname, Array(ta.resid),
      ta.atomname, Array(ta.atomid),
      Array(ta.mass), Array(ta.charge),
      Array(ta.list_bond), Array(ta.list_angle), Array(ta.list_dihedral),
      Array(ta.list_improper), Array(ta.list_cmap))
end

function cpu64(ta::TrjArray)
    TrjArray{Float64, Int64}(Array{Float64}(ta.x), Array{Float64}(ta.y), Array{Float64}(ta.z), Array{Float64}(ta.boxsize),
      ta.chainname, Array{Int64}(ta.chainid),
      ta.resname, Array{Int64}(ta.resid),
      ta.atomname, Array{Int64}(ta.atomid),
      Array{Float64}(ta.mass), Array{Float64}(ta.charge),
      Array{Int64}(ta.list_bond), Array{Int64}(ta.list_angle), Array{Int64}(ta.list_dihedral),
      Array{Int64}(ta.list_improper), Array{Int64}(ta.list_cmap))
end


text:
"""
    FEMMDeforWinklerModule

Module for operations on boundaries of domains to construct system matrices and
system vectors for linear deformation models with distributed-spring supports
(Winkler foundation model).
"""
module FEMMDeforWinklerModule

__precompile__(true)

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, nodesperelem, manifdim
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobiansurface
import FinEtools.FieldModule: ndofs, gatherdofnums!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.FEMMBaseModule: AbstractFEMM
import FinEtools.AssemblyModule: AbstractSysvecAssembler, AbstractSysmatAssembler, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!
import FinEtools.MatrixUtilityModule: add_nnt_ut_only!, complete_lt!, locjac!
import FinEtools.SurfaceNormalModule: SurfaceNormal, updatenormal!
import LinearAlgebra: norm, cross

"""
    FEMMDeforWinkler{S<:AbstractFESet, F<:Function} <: AbstractFEMM

Type for normal spring support  (Winkler).
"""
mutable struct FEMMDeforWinkler{S<:AbstractFESet, F<:Function} <: AbstractFEMM
    integdomain::IntegDomain{S, F} # geometry data
end

"""
    surfacenormalspringstiffness(self::FEMMDeforWinkler, assembler::A,
        geom::NodalField{FFlt}, u::NodalField{T},
        springconstant::FFlt, surfacenormal::SurfaceNormal) where {A<:AbstractSysmatAssembler, T<:Number}

Compute the stiffness matrix of surface normal spring.

Rationale: consider continuously distributed springs between the surface of the
solid body and the 'ground', in the direction normal to the surface. If the
spring coefficient becomes large, we have an approximate method of enforcing the
normal displacement to the surface.gas
"""
function surfacenormalspringstiffness(self::FEMMDeforWinkler, assembler::A,
    geom::NodalField{FFlt}, u::NodalField{T},
    springconstant::FFlt, surfacenormal::SurfaceNormal) where {A<:AbstractSysmatAssembler, T<:Number}
    integdomain = self.integdomain
    # Constants
    nfes = count(integdomain.fes); # number of finite elements in the set
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(integdomain.fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(integdomain.fes); # manifold dimension of the element
    Kedim = ndn*nne;             # dimension of the element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc = integrationdata(integdomain);
    # Prepare assembler and temporaries
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    Ke = zeros(FFlt,Kedim,Kedim);                # element matrix -- used as a buffer
    dofnums = zeros(FInt, Kedim); # degree of freedom array -- used as a buffer
    loc = zeros(FFlt, 1,sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim,mdim); # Jacobian matrix -- used as a buffer
    Nn = zeros(FFlt, Kedim); # column vector
    startassembly!(assembler, Kedim, Kedim, nfes, u.nfreedofs, u.nfreedofs);
    for i = 1:nfes # Loop over elements
        gathervalues_asmat!(geom, ecoords, integdomain.fes.conn[i]);
        fill!(Ke, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(integdomain, J, loc, integdomain.fes.conn[i], Ns[j]);
            n = updatenormal!(surfacenormal, loc, J, integdomain.fes.label[i])
            for k= 1:nne
                for r = 1:sdim
                    Nn[(k-1)*sdim+r] = n[r] * Ns[j][k];
                end
            end
            add_nnt_ut_only!(Ke, Nn, springconstant*Jac*w[j])
        end # Loop over quadrature points
        complete_lt!(Ke)
        gatherdofnums!(u, dofnums, integdomain.fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, Ke, dofnums, dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function surfacenormalspringstiffness(self::FEMMDeforWinkler,
              geom::NodalField{FFlt}, u::NodalField{T},
              springconstant::FFlt, surfacenormal::SurfaceNormal) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return surfacenormalspringstiffness(self, assembler, geom, u, springconstant, surfacenormal);
end

end


text:
struct Foo

end

function bar(self::Foo)::Int64
    return baz(self)
end

function baz(self::Foo)::Int64
    return 10
end

function main()
    f = Foo()
    b = bar(f)
    println(join([b], " "))
    @assert(b == 10)
end

main()


text:
export get_agent_view, AbstractGridWorld
export get_actions
abstract type AbstractGridWorld end

function get_agent_view end
function get_agent end

Base.convert(::Type{GridWorldBase}, w::AbstractGridWorld) = w.world
get_object(w::AbstractGridWorld) = get_object(convert(GridWorldBase, w))
get_object(w::AbstractGridWorld, x::Type{<:AbstractObject}) = filter(o -> o isa x, get_object(w))
get_object(w::AbstractGridWorld, x::Type{Agent}) = w.agent
get_pos(w::AbstractGridWorld, ::Type{Agent}) = w.agent_pos

get_agent(w::AbstractGridWorld) = get_object(w, Agent)
get_agent_pos(w::AbstractGridWorld) = get_pos(w, Agent)
get_agent_dir(w::AbstractGridWorld) = w |> get_agent |> get_dir

function get_agent_view(w::AbstractGridWorld, agent_view_size=(7,7))
    wb = convert(GridWorldBase, w)
    v = BitArray{3}(undef, size(wb, 1), agent_view_size...)
    fill!(v, false)
    get_agent_view!(v, w)
end

function (w::AbstractGridWorld)(action::Union{TurnRight, TurnLeft})
    a = get_agent(w)
    set_dir!(a, action(get_dir(a)))
    w
end

get_actions(w::AbstractGridWorld) = (MOVE_FORWARD, TURN_LEFT, TURN_RIGHT)

get_agent_view_inds(w::AbstractGridWorld, s=(7,7)) = get_agent_view_inds(get_agent_pos(w).I, s, get_agent_dir(w))

get_agent_view!(v::BitArray{3}, w::AbstractGridWorld) = get_agent_view!(v, convert(GridWorldBase, w), get_agent_pos(w), get_agent_dir(w))


text:
# If want to explore ternary plots and line fitting:

using JLD2, DataFrames, CSV, Dates, Query
using LinearAlgebra, GLM
using RCall
using Gadfly, Colors

include("/home/kristen/Documents/coda_utilities/coda_ops_transforms.jl")

#load in the data:
@load "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_df_zero_plus.jld2"
#should return dataframe oligo_df_plus

## Ternary Plots, courtesy of the lovely Ternary package in R:

R"require(Ternary)"
R"require(png)"

x136=closure(convert(Array,oligo_df_plus[!,[:Oligo1, :Oligo3,:Oligo6]])) #convert Dataframe columns to array
R"""
TernaryPlot()
TernaryPlot(atip = "O1-IC", btip = "O3-II/XV", ctip = "O6-I?", lab.col=c("#0000FF","#FF7F00","#000099"),tip.col=c("#0000FF","#FF7F00","#000099"))
TernaryPoints($x136, pch=16, cex=0.8, col='black')
"""

x245=closure(convert(Array,oligo_df_plus[!,[:csO2, :csO4,:csO5]]))
R"""
TernaryPlot()
TernaryPlot(atip = "O2-CB5",btip = "O4-III/IV", ctip = "O5-IE", lab.col=c("#CC334C","#009900","#00CCFF"),tip.col=c("#CC334C","#009900","#00CCFF"),grid.lines = 10, grid.minor.lines = 0, grid.col = "lightgrey")
TernaryPoints($x245, pch=16, cex=0.8, col='black')
"""

## Fancier year day colors:

x123=closure(convert(Array,oligo_df_plus[!,[:csO1,:csO2,:csO3]]))
yearday=dayofyear.(oligo_df_plus[!,:date])

R"""
TernaryPlot()
jet <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                     "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
jmap <- jet(365)[$yearday] #first argument is max....I think...
TernaryPlot(atip = "O1-I", btip = "O2-CB5", ctip = "O3-II/XV", lab.col=c("#0000FF","#CC334C","#FF7F00"),tip.col=c("#0000FF","#CC334C","#FF7F00"),grid.lines = 10, grid.minor.lines = 0, grid.col = "lightgrey")
TernaryPoints($x123, pch=16, cex=0.8, col=jmap) #whoo hoo!!!!!
"""

## Fit a line through this data in the simplex

#Note, this is easiest to do with an ilr tranform, fitting the line with standard linear regression and then back transfrom to simplex:

#Line throuhg oligotype 1,2,3 data:
x123=closure(convert(Array,oligo_df_plus[!,[:csO1,:csO2,:csO3]]))
sbp=[-1 -1 1; 1 -1 -0] #to form basis
ilr_x, bal=ilr(x123,sbp) #ilr transform

#see what this transformed data looks like:
Gadfly.set_default_plot_size(35cm,25cm)
layer1=layer(x=ilr_x[:,2], y=ilr_x[:,1],Geom.point, style(default_color=RGB(0,0,1))) #okay, that seems to be working...
plt=plot(layer1)

#fit line:
T1=DataFrame(x = ilr_x[:,2], y = ilr_x[:,1])
test1 = lm(@formula(y ~ x), T1)
bf=zeros(2,1)
bf[2] = GLM.coef(test1)[2]
bf[1] = GLM.coef(test1)[1]

xf=collect(-2:0.1:8) #and see what this fit looks like:
yf=bf[1] .+ bf[2] .* xf
layer2=layer(x=xf, y=yf, Geom.line)
push!(plt,layer2) #check that line goes through points!

## Now plot this line back in the simplex and save:

ll=invilr([yf xf],bal) #back transform fitted line

R"""
TernaryPlot()
jet <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                     "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
jmap <- jet(365)[$yearday] #first argument is max....I think...
TernaryPlot(atip = "O1-I", btip = "O2-CB5", ctip = "O3-II/XV", lab.col=c("#0000FF","#CC334C","#FF7F00"),tip.col=c("#0000FF","#CC334C","#FF7F00"),grid.lines = 10, grid.minor.lines = 0, grid.col = "lightgrey")
TernaryPoints($x123, pch=16, cex=0.8, col=jmap) #whoo hoo!!!!!
TernaryLines($ll, col='black')
"""

#and if wanted to save:
R"""
png(file="/home/kristen/Documents/V6V8_analysis/paper_V6V8/V6V8_Syn_CoDa_paper/figures/tern1.png")
TernaryPlot()
jet <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                     "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
jmap <- jet(365)[$yearday] #first argument is max....I think...
TernaryPlot(atip = "O1-I", btip = "O2-CB5", ctip = "O3-II/XV", lab.col=c("#0000FF","#CC334C","#FF7F00"),tip.col=c("#0000FF","#CC334C","#FF7F00"),grid.lines = 10, grid.minor.lines = 0, grid.col = "lightgrey")
TernaryPoints($x123, pch=16, cex=0.8, col=jmap) #whoo hoo!!!!!
TernaryLines($ll, col='black')
dev.off()
"""

## Again for the remaining set:

x456=closure(convert(Array,oligo_df_plus[!,[:csO4,:csO5,:csO6]]))
sbp=[-1 1 -1; 1 0 -1]
ilr_x, bal=ilr(x456,sbp)
layer1=layer(x=ilr_x[:,2], y=ilr_x[:,1],Geom.point, style(default_color=RGB(0,0,1))) #okay, that seems to be working...
plt=plot(layer1)

#fit line just within date:
T1=DataFrame(x = ilr_x[:,2], y = ilr_x[:,1])
test1 = lm(@formula(y ~ x), T1)
bf=zeros(2,1)
bf[2] = GLM.coef(test1)[2]
bf[1] = GLM.coef(test1)[1]

xf=collect(-3:0.1:5)
yf=bf[1] .+ bf[2] .* xf
layer2=layer(x=xf, y=yf, Geom.line)
push!(plt,layer2) #check that line goes through points!

ll=invilr([yf xf],bal) #back transform fitted line...

R"""
png(file="/home/kristen/Documents/V6V8_analysis/paper_V6V8/V6V8_Syn_CoDa_paper/figures/tern2.png")
TernaryPlot()
jet <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                     "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
jmap <- jet(365)[$yearday] #first argument is max....I think...
TernaryPlot(atip = "O4-III/IV", btip = "O5-IE", ctip = "O6-I",lab.col=c("#009900","#00CCFF","#000099"),tip.col=c("#009900","#00CCFF","#000099"),grid.lines = 10, grid.minor.lines = 0, grid.col = "lightgrey")
TernaryPoints($x456, pch=16, cex=0.8, col=jmap) #whoo hoo!!!!!
TernaryLines($ll, col='black')
dev.off()
"""

## to get the simplified relationship:

a=ll[:,1]./ll[:,2]
b=ll[:,2]./ll[:,3]

# a=ll[:,1]./ll[:,3]
# b=ll[:,3]./ll[:,2]

T1=DataFrame(x = log.(a), y = log.(b))

test2 = lm(@formula(y ~ x), T1)
bf=zeros(2,1)
bf[2] = GLM.coef(test2)[2]
bf[1] = GLM.coef(test2)[1]

plot(x=log.(a),y=log.(b),Geom.point)


text:
@doc raw""" ```julia
struct ALSolver <: ConstrainedSolver{T}
```
Augmented Lagrangian (AL) is a standard tool for constrained optimization. For a trajectory optimization problem of the form:
```math
\begin{aligned}
  \min_{x_{0:N},u_{0:N-1}} \quad & \ell_f(x_N) + \sum_{k=0}^{N-1} \ell_k(x_k, u_k, dt) \\
  \textrm{s.t.}            \quad & x_{k+1} = f(x_k, u_k), \\
                                 & g_k(x_k,u_k) \leq 0, \\
                                 & h_k(x_k,u_k) = 0.
\end{aligned}
```
AL methods form the following augmented Lagrangian function:
```math
\begin{aligned}
    \ell_f(x_N) + &λ_N^T c_N(x_N) + c_N(x_N)^T I_{\mu_N} c_N(x_N) \\
           & + \sum_{k=0}^{N-1} \ell_k(x_k,u_k,dt) + λ_k^T c_k(x_k,u_k) + c_k(x_k,u_k)^T I_{\mu_k} c_k(x_k,u_k)
\end{aligned}
```
This function is then minimized with respect to the primal variables using any unconstrained minimization solver (e.g. iLQR).
    After a local minima is found, the AL method updates the Lagrange multipliers λ and the penalty terms μ and repeats the unconstrained minimization.
    AL methods have superlinear convergence as long as the penalty term μ is updated each iteration.
"""
struct ALSolver{T,S} <: ConstrainedSolver{T}
    opts::SolverOptions{T}
    stats::SolverStats{T}
    ilqr::S
end

function ALSolver(
        prob::Problem{T},
        opts::SolverOptions=SolverOptions(),
        stats::SolverStats=SolverStats(parent=solvername(ALSolver));
        use_static=usestaticdefault(get_model(prob)[1]),
        kwarg_opts...
    ) where {T}
    set_options!(opts; kwarg_opts...)

    # Build Augmented Lagrangian Objective
    alobj = ALObjective{T}(prob.obj)
    nx,nu = dims(prob)
    prob_al = Problem(prob.model, alobj, ConstraintList(nx,nu),
        prob.x0, prob.xf, prob.Z, prob.N, TO.get_initial_time(prob), TO.get_final_time(prob))

    # Instantiate the iLQR solver
    ilqr = iLQRSolver(prob_al, opts, stats, use_static=use_static)
    initialize!(alobj.conset, prob.constraints, ilqr.Z, ilqr.opts, alobj.alcost, ilqr.Efull)
    # settraj!(alobj.conset, get_trajectory(ilqr))

    # Build the solver
    solver = ALSolver(opts, stats, ilqr)
    reset!(solver)
    return solver
end

# Getters
for method in (:(RD.dims), :(RD.state_dim), :(RD.errstate_dim), :(RD.control_dim),
              :(TO.get_trajectory), :(TO.get_objective), :(TO.get_model),
              :(TO.get_initial_state), :getlogger)
    @eval $method(solver::ALSolver) = $method(solver.ilqr)
end
solvername(::Type{<:ALSolver}) = :AugmentedLagrangian
get_ilqr(solver::ALSolver) = solver.ilqr
TO.get_constraints(solver::ALSolver) = solver.ilqr.obj.conset
stats(solver::ALSolver) = solver.stats
options(solver::ALSolver) = solver.opts

# Methods
function reset!(solver::ALSolver)
    # reset_solver!(solver)
    opts = options(solver)::SolverOptions
    reset!(stats(solver), opts.iterations, solvername(solver))
    reset!(solver.ilqr)

    # Reset constraints
    conset = get_constraints(solver)
    reset!(conset)
end

function max_violation(solver::ALSolver)
    evaluate_constraints!(get_constraints(solver), get_trajectory(solver))
    max_violation(get_constraints(solver))
end


text:
# C. Cooperative Groups

export this_grid, sync_grid

"""
    this_grid()

Returns a `grid_handle` of the grid group this thread belongs to. Only available if a
cooperative kernel is launched.
"""
this_grid() = cudaCGGetIntrinsicHandle(cudaCGScopeGrid)

"""
    sync_grid(grid_handle::Culonglong)

Waits until all threads in all blocks in the grid `grid_handle` have reached this point and
all global memory accesses made by these threads prior to `sync_grid()` are visible to all
threads in the grid. A 32-bit integer `cudaError_t` is returned.
"""
sync_grid(handle) = cudaCGSynchronize(handle, 0)


text:
using CompScienceMeshes, BEAST
using LinearAlgebra
using Profile
using StaticArrays

function tau(x::SVector{U,T}) where {U,T}
    4.0-1.0
end

ttrc = X->ttrace(X,y)

T= tetmeshsphere(1.0,0.25)
X = nedelecc3d(T)
y = boundary(T)
@show numfunctions(X)

ϵ, μ, ω = 1.0, 1.0, 1.0; κ, η = ω * √(ϵ*μ), √(μ/ϵ)
ϵ_r = 4.0
χ = tau
K, I, B = VIE.singlelayer2(wavenumber=κ, tau=χ), Identity(), VIE.boundary2(wavenumber=κ, tau=χ)
E = VIE.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
H = -1/(im*κ*η)*curl(E)

@hilbertspace j
@hilbertspace k
α = ϵ_r
eq = @varform α*I[j,k]-K[j,k]+B[ttrc(j),k] == E[j]

evie = @discretise eq j∈X k∈X
u = solve(evie)

#postprocessing
Φ, Θ = [0.0], range(0,stop=2π,length=100)
pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ ]
ffe = potential(VIE.farfield(wavenumber=κ, tau=χ), pts, u, X)
ff = ffe

using Plots

#Farfield
plot(xlabel="theta")
plot!(Θ, norm.(ff), label="far field", title="E-VIE")

#Nearfield
Z = range(-1,1,length=100)
Y = range(-1,1,length=100)
nfpoints = [point(0,y,z) for  y in Y, z in Z]

Enear = BEAST.grideval(nfpoints,u,X)
Enear = reshape(Enear,100,100)

contour(real.(getindex.(Enear,1)))
heatmap(Z, Y, real.(getindex.(Enear,1)), clim=(0.0,1.0))


text:
module ShortStrings

using BitIntegers
using SortingAlgorithms

export fsort, fsort!, ShortString, ShortString3, ShortString7, ShortString15
export ShortString31, ShortString63, ShortString127, ShortString255
export @ss_str, @ss3_str, @ss7_str, @ss15_str, @ss31_str, @ss63_str, @ss127_str, @ss255_str

include("base.jl")

using MurmurHash3: mmhash128_a

function Base.hash(x::ShortString, h::UInt)
    h += Base.memhash_seed
    last(mmhash128_a(sizeof(x), bswap(x.size_content), h%UInt32)) + h
end

end # module


text:
let
    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)

    scene = Scene()
    get!(scene, trajdata, 1)
    @test length(scene) == 2
    for (i,veh) in enumerate(scene)
        @test scene[i].state == get_state(trajdata, i, 1)
        @test scene[i].def == get_def(trajdata, i)
    end

    scene2 = Scene(deepcopy(scene.entities), 2)
    @test length(scene2) == 2
    for (i,veh) in enumerate(scene2)
        @test scene2[i].state == get_state(trajdata, i, 1)
        @test scene2[i].def == get_def(trajdata, i)
    end

    @test get_by_id(scene, 1) == scene[1]

    empty!(scene2)
    @test length(scene2) == 0

    copy!(scene2, scene)
    @test length(scene2) == 2
    for (i,veh) in enumerate(scene2)
        @test scene2[i].state == get_state(trajdata, i, 1)
        @test scene2[i].def == get_def(trajdata, i)
    end

    delete!(scene2, scene2[1])
    @test length(scene2) == 1
    @test scene2[1].state == get_state(trajdata, 2, 1)
    @test scene2[1].def == get_def(trajdata, 2)
    scene2[1] = deepcopy(scene[1])
    @test scene2[1].state == get_state(trajdata, 1, 1)
    @test scene2[1].def == get_def(trajdata, 1)

    @test findfirst(scene, 1) == 1
    @test findfirst(scene, 2) == 2
    @test get_first_available_id(scene) == 3

    @test in(scene, 1)
    @test in(scene, 2)
    @test !in(scene, 3)

    veh = scene[2]
    @test veh.state == get_state(trajdata, 2, 1)
    @test veh.def == get_def(trajdata, 2)

    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 1), 1, roadway) == NeighborLongitudinalResult(2, 3.0078125)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 1), 2, roadway) == NeighborLongitudinalResult(0, 250.0)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 2), 1, roadway) == NeighborLongitudinalResult(2, 4.0)
    @test get_neighbor_fore_along_lane(get!(Scene(), trajdata, 2), 2, roadway) == NeighborLongitudinalResult(0, 250.0)

    push!(scene, get_state(trajdata, 1, 1))
end

let
    roadway = gen_stadium_roadway(1)
    scene = Scene(2)
    scene.n = 2
    def = VehicleDef(AgentClass.CAR, 2.0, 1.0)

    function place_at!(i, s)
        roadproj = proj(VecSE2(0.0,0.0,0.0), roadway)
        roadind = RoadIndex(roadproj.curveproj.ind, roadproj.tag)
        roadind = move_along(roadind, roadway, s)
        frenet = Frenet(roadind, roadway[roadind].s, 0.0, 0.0)
        state = VehicleState(frenet, roadway, 0.0)
        scene[i] = Vehicle(state, def, i)
    end

    place_at!(1, 0.0)
    place_at!(2, 0.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test isapprox(foreinfo.Δs, 0.0)

    place_at!(2, 100.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 100.0)

    foreinfo = get_neighbor_fore_along_lane(scene, 2, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 1
    @test isapprox(foreinfo.Δs, 277.07, atol=1e-2)

    place_at!(2, 145.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 145.0, atol=1e-5)

    place_at!(2, 240.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 240.0, atol=1e-5)

    place_at!(1, 240.0)
    foreinfo = get_neighbor_fore_along_lane(scene, 1, roadway, max_distance_fore=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 0.0, atol=1e-5)

    ################################################
    # get_frenet_relative_position

    place_at!(1, 0.0)
    place_at!(2, 10.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.origin == scene[1].state.posF.roadind
    @test frp.target.ind.i == 0
    @test isapprox(frp.target.ind.t, 0.1, atol=0.01)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 10.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 10.0)
    place_at!(2, 20.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.origin == scene[1].state.posF.roadind
    @test frp.target.ind.i == 0
    @test isapprox(frp.target.ind.t, 0.2, atol=0.01)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 10.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 0.0)
    place_at!(2, 120.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 120.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 0.0)
    place_at!(2, 250.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.target.tag != scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 250.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)
end


text:
# compute exp(pi*i*q),
# which is algebraic if q is rational.
function exp_alg(q::Rational)
       qh = q/2

       # first, obtain minimal polynomial
       s, x = PolynomialRing(AlgebraicNumbers.Nemo.ZZ, "x")
       poly = cyclotomic(denominator(qh), x)
       coeffs = [convert(BigInt, coeff(poly, i)) for i=0: denominator(qh)]

       # now, select root.
       apprx = exp(im*BigFloat(pi)*q)

       # Finally, return polynomial w.r.t. that root
       return AlgebraicNumber(coeffs, apprx)
end

cos_alg(q::Rational) = real(exp_alg(q))
sin_alg(q::Rational) = imag(exp_alg(q))

# checks if a polynomial is equal to any of the
# cyclotomic polynomials n in a list of candidates
function is_cyclotomic(poly, candidates, x)
	for n in candidates
		if poly == cyclotomic(n, x)
			return (true, n)
		end
	end
	return (false, 0)
end

# compute log(a)/(pi*i),
# which is rational if a is a root of unity.
# If a is not a root of unity, returns Nothing
function log_alg(a::AlgebraicNumber)
	s, x = PolynomialRing(Nemo.ZZ, "x")

	deg = length(a.coeff)-1
	poly = s(map(Nemo.ZZ, a.coeff))

	(is_cycl, denom) = is_cyclotomic(poly, inv_totient(deg), x)

	if is_cycl
		# TODO: should this be BigInt?
		num = round(Int, denom*imag(log(a.apprx)/pi))
		return num//denom
	else
		return Nothing
	end
end

function acos_alg(x::AlgebraicNumber)
	# TODO: check if the number can be made a root of unity.

	# First, make the number a root of unity.
	y = sqrt(1 - x^2)
	z = x + sqrt(AlgebraicNumber(-1))*y

	# Now take log
	log_alg(z)
end

function asin_alg(x::AlgebraicNumber)
	# TODO: check if the number can be made a root of unity.

	# First, make the number a root of unity.
	y = sqrt(1 - x^2)
	z = y + sqrt(AlgebraicNumber(-1))*x

	# Now take log
	log_alg(z)
end


text:
function testdistl1()
    println("\n in testdistl1")
    dim = 10
    hnsw = createHnswApi(Float32, 8, 16, "DistL1")
    # block // insertion
    datas = rand(Float32, (dim, 500));
    data_insert = map(i -> ( rand(Float32, dim) , UInt(i) ), 1:size(datas)[2])
    insert(hnsw, data_insert)
    # one insertion
    v = rand(Float32, dim)
    insert(hnsw, v, UInt(501))
    # testing search
    v1 = [Float32(0.5) for i in 1:dim]
    neighbours = search(hnsw, v1, 10, 16)
    # testing block // search
    datas = map(i -> rand(Float32, dim), 1:100)
    neighbours = search(hnsw, datas, 10, 16)
    true
end

#  rust calling so length will be passes as a usize.

function mydist(pa::Ptr{Float32}, pb::Ptr{Float32}, l::UInt64)
#    println("in mydist")
    #
   dist = 0
    #
    va = unsafe_wrap(Array{Float32,1}, pa, NTuple{1,Int64}(Int64(l)); own = false)
    vb = unsafe_wrap(Array{Float32,1}, pb, NTuple{1,Int64}(Int64(l)); own = false)
    for i in  1:l
        dist += abs(va[i] - vb[i])
    end
    return Float32(dist)
end

# test only serial insertion/search until julia 1.3

function testdistptr()
    println("\n\n in testdistptr")
    dim = 10
    mydist_ptr = Base.@cfunction(mydist, Cfloat, (Ptr{Cfloat}, Ptr{Cfloat}, Culonglong))
    hnsw = createHnswApi(Float32, 8, 16, mydist_ptr)
    # block // insertion
    nbinsert = 100
    for i in 1:nbinsert
        data = rand(Float32, dim);
        insert(hnsw, data, UInt64(i))
    end
    # testing search
    println("testing one search")
    v1 = [Float32(0.5) for i in 1:dim]
    neighbours = search(hnsw, v1, 8, 16)
    @debug "neighbours : " neighbours
    true
end


text:
include(joinpath(dirname(@__FILE__), "config.jl"))
include(joinpath(dirname(@__FILE__), "lstm.jl"))
include(joinpath(dirname(@__FILE__), "seq-data.jl"))

# build vocabulary
vocab   = build_vocabulary(INPUT_FILE, VOCAB_FILE)
n_class = length(vocab)

#--LSTM
# define LSTM
lstm = LSTM(LSTM_N_LAYER, SEQ_LENGTH, DIM_HIDDEN, DIM_EMBED,
            n_class, dropout=DROPOUT, name=NAME)
#--/LSTM

#--data
# load data
text_all  = readall(INPUT_FILE)
len_train = round(Int, length(text_all)*DATA_TR_RATIO)
text_tr   = text_all[1:len_train]
text_val  = text_all[len_train+1:end]

data_tr   = CharSeqProvider(text_tr, BATCH_SIZE, SEQ_LENGTH, vocab, NAME,
                            LSTM_N_LAYER, DIM_HIDDEN)
data_val  = CharSeqProvider(text_val, BATCH_SIZE, SEQ_LENGTH, vocab, NAME,
                            LSTM_N_LAYER, DIM_HIDDEN)
#--/data

# set up training
if USE_GPU
  context = [mx.gpu(i) for i = 0:N_GPU-1]
else
  context = [mx.cpu()]
end

#--train
model = mx.FeedForward(lstm, context=context)
optimizer = mx.ADAM(lr=BASE_LR, weight_decay=WEIGHT_DECAY, grad_clip=CLIP_GRADIENT)

mx.fit(model, optimizer, data_tr, eval_data=data_val, n_epoch=N_EPOCH,
       initializer=mx.UniformInitializer(0.1),
       callbacks=[mx.speedometer(), mx.do_checkpoint(CKPOINT_PREFIX)], eval_metric=NLL())
#--/train


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

# -----------------
# BASIC TRANSFORMS
# -----------------

include("transforms/closure.jl")
include("transforms/remainder.jl")

# ---------------------
# LOG-RATIO TRANSFORMS
# ---------------------

"""
    LogRatio

A log-ratio transform following the
[TableTransforms.jl](https://github.com/JuliaML/TableTransforms.jl)
interface.

See also [`ALR`](@ref), [`CLR`](@ref), [`ILR`](@ref).
"""
abstract type LogRatio <: Stateless end

assertions(::Type{<:LogRatio}) = [TableTransforms.assert_continuous]

isrevertible(::Type{<:LogRatio}) = true

function apply(transform::LogRatio, table)
  # basic checks
  for assertion in assertions(transform)
    assertion(table)
  end

  # original variable names
  vars = Tables.columnnames(table)

  # reference variable
  rvar = refvar(transform, vars)
  @assert rvar ∈ vars "invalid reference variable"
  rind = first(indexin([rvar], collect(vars)))

  # permute columns if necessary
  ovars  = setdiff(vars, (rvar,))
  pvars  = [ovars; rvar]
  ptable = table |> Select(pvars)

  # transformation
  X = Tables.matrix(ptable)
  Y = applymatrix(transform, X)

  # new variable names
  n = newvars(transform, pvars)

  # return same table type
  𝒯 = (; zip(n, eachcol(Y))...)
  newtable = 𝒯 |> Tables.materializer(table)

  newtable, (rvar, rind)
end

function revert(transform::LogRatio, table, cache)
  # retrieve cache
  rvar, rind = cache

  # trasformation
  Y = Tables.matrix(table)
  X = revertmatrix(transform, Y)

  # original variable names
  vars = Tables.columnnames(table)
  n = oldvars(transform, vars, rvar)

  # permute reference variable
  n[[rind,end]]   .= n[[end,rind]]
  X[:,[rind,end]] .= X[:,[end,rind]]

  # return same table type
  𝒯 = (; zip(n, eachcol(X))...)
  𝒯 |> Tables.materializer(table)
end

# to be implemented by log-ratio transforms
function refvar end
function newvars end
function oldvars end
function applymatrix end
function revertmatrix end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("transforms/alr.jl")
include("transforms/clr.jl")
include("transforms/ilr.jl")


text:
"""
Julia-fjåge gateway, standalone container, and slave container.

Notes:
- This implementation is not thread-safe.
"""
module Fjage

using Sockets, Distributed, Base64, UUIDs, Dates, Logging, Random
using JSON, MacroTools, PkgVersion
using InteractiveUtils: subtypes

export BLOCKING

const BLOCKING = -1
const MAX_QUEUE_LEN = 256
const VERSION = @PkgVersion.Version

const timeout = Ref(1000)

"""
    default_timeout(millis)

Set default timeout for requests. This timeout is used when not explicitly
specified. The timeout is given in milliseconds.
"""
default_timeout(t) = (timeout[] = t)

include("aid.jl")
include("msg.jl")
include("const.jl")
include("gw.jl")
include("container.jl")

end


text:
using ZXCalculus
using ZXCalculus: SpiderType
st = Dict(
    0+1 => SpiderType.In,
    1+1 => SpiderType.In,
    2+1 => SpiderType.In,
    3+1 => SpiderType.In,
    4+1 => SpiderType.In,
    5+1 => SpiderType.Z,
    6+1 => SpiderType.Z,
    7+1 => SpiderType.Z,
    8+1 => SpiderType.Z,
    9+1 => SpiderType.Z,
    10+1 => SpiderType.Z,
    11+1 => SpiderType.Z,
    12+1 => SpiderType.Z,
    13+1 => SpiderType.Z,
    14+1 => SpiderType.Z,
    15+1 => SpiderType.Z,
    16+1 => SpiderType.Z,
    17+1 => SpiderType.Z,
    18+1 => SpiderType.Z,
    19+1 => SpiderType.Z,
    20+1 => SpiderType.Z,
    21+1 => SpiderType.Z,
    22+1 => SpiderType.Z,
    23+1 => SpiderType.Out,
    24+1 => SpiderType.Out,
    25+1 => SpiderType.Out,
    26+1 => SpiderType.Out,
    27+1 => SpiderType.Out,
    28+1 => SpiderType.Z,
    29+1 => SpiderType.Z,
    30+1 => SpiderType.Z,
    31+1 => SpiderType.Z,
    32+1 => SpiderType.Z,
    33+1 => SpiderType.Z,
    34+1 => SpiderType.Z,
    35+1 => SpiderType.Z,
    36+1 => SpiderType.Z,
    37+1 => SpiderType.Z,
    38+1 => SpiderType.Z,
    39+1 => SpiderType.Z,
    40+1 => SpiderType.Z,
    41+1 => SpiderType.Z,
    42+1 => SpiderType.Z,
    43+1 => SpiderType.Z,
    44+1 => SpiderType.Z,
    45+1 => SpiderType.Z,
    46+1 => SpiderType.Z,
    47+1 => SpiderType.Z,
    48+1 => SpiderType.Z,
    49+1 => SpiderType.Z,
    50+1 => SpiderType.Z,
    51+1 => SpiderType.Z,
)

ps = Dict(
    0+1 => 0//1,
    1+1 => 0//1,
    2+1 => 0//1,
    3+1 => 0//1,
    4+1 => 0//1,
    5+1 => 0//1,
    6+1 => 5//4,
    7+1 => 0//1,
    8+1 => 5//4,
    9+1 => 3//2,
    10+1 => 0//1,
    11+1 => 0//1,
    12+1 => 0//1,
    13+1 => 1//4,
    14+1 => 0//1,
    15+1 => 0//1,
    16+1 => 1//4,
    17+1 => 0//1,
    18+1 => 7//4,
    19+1 => 0//1,
    20+1 => 0//1,
    21+1 => 0//1,
    22+1 => 3//2,
    23+1 => 0//1,
    24+1 => 0//1,
    25+1 => 0//1,
    26+1 => 0//1,
    27+1 => 0//1,
    28+1 => 1//2,
    29+1 => 7//4,
    30+1 => 1//2,
    31+1 => 1//4,
    32+1 => 1//2,
    33+1 => 1//4,
    34+1 => 1//2,
    35+1 => 1//2,
    36+1 => 7//4,
    37+1 => 1//2,
    38+1 => 0//1,
    39+1 => 1//4,
    40+1 => 0//1,
    41+1 => 3//4,
    42+1 => 0//1,
    43+1 => 1//4,
    44+1 => 0//1,
    45+1 => 3//4,
    46+1 => 0//1,
    47+1 => 7//4,
    48+1 => 0//1,
    49+1 => 0//1,
    50+1 => 0//1,
    51+1 => 0//1
)

es = Dict(
    (0+1, 8+1) => EdgeType.HAD,
    (1+1, 6+1) => EdgeType.HAD,
    (2+1, 13+1) => EdgeType.HAD,
    (3+1, 16+1) => EdgeType.HAD,
    (4+1, 5+1) => EdgeType.HAD,
    (5+1, 6+1) => EdgeType.HAD,
    (5+1, 8+1) => EdgeType.HAD,
    (5+1, 9+1) => EdgeType.HAD,
    (6+1, 8+1) => EdgeType.HAD,
    (6+1, 10+1) => EdgeType.HAD,
    (6+1, 11+1) => EdgeType.HAD,
    (6+1, 13+1) => EdgeType.HAD,
    (6+1, 18+1) => EdgeType.HAD,
    (6+1, 19+1) => EdgeType.HAD,
    (6+1, 20+1) => EdgeType.HAD,
    (6+1, 40+1) => EdgeType.HAD,
    (6+1, 42+1) => EdgeType.HAD,
    (6+1, 46+1) => EdgeType.HAD,
    (6+1, 49+1) => EdgeType.HAD,
    (7+1, 8+1) => EdgeType.HAD,
    (7+1, 9+1) => EdgeType.HAD,
    (7+1, 30+1) => EdgeType.HAD,
    (8+1, 10+1) => EdgeType.HAD,
    (8+1, 13+1) => EdgeType.HAD,
    (8+1, 18+1) => EdgeType.HAD,
    (8+1, 19+1) => EdgeType.HAD,
    (8+1, 21+1) => EdgeType.HAD,
    (8+1, 38+1) => EdgeType.HAD,
    (8+1, 40+1) => EdgeType.HAD,
    (8+1, 46+1) => EdgeType.HAD,
    (8+1, 48+1) => EdgeType.HAD,
    (9+1, 10+1) => EdgeType.HAD,
    (9+1, 11+1) => EdgeType.HAD,
    (9+1, 13+1) => EdgeType.HAD,
    (9+1, 18+1) => EdgeType.HAD,
    (9+1, 38+1) => EdgeType.HAD,
    (9+1, 42+1) => EdgeType.HAD,
    (9+1, 44+1) => EdgeType.HAD,
    (9+1, 46+1) => EdgeType.HAD,
    (10+1, 35+1) => EdgeType.HAD,
    (11+1, 32+1) => EdgeType.HAD,
    (12+1, 16+1) => EdgeType.HAD,
    (12+1, 18+1) => EdgeType.HAD,
    (12+1, 33+1) => EdgeType.HAD,
    (13+1, 14+1) => EdgeType.HAD,
    (13+1, 15+1) => EdgeType.HAD,
    (13+1, 17+1) => EdgeType.HAD,
    (13+1, 22+1) => EdgeType.HAD,
    (13+1, 50+1) => EdgeType.HAD,
    (14+1, 16+1) => EdgeType.HAD,
    (14+1, 18+1) => EdgeType.HAD,
    (14+1, 36+1) => EdgeType.HAD,
    (15+1, 16+1) => EdgeType.HAD,
    (15+1, 29+1) => EdgeType.HAD,
    (16+1, 51+1) => EdgeType.HAD,
    (17+1, 18+1) => EdgeType.HAD,
    (17+1, 31+1) => EdgeType.HAD,
    (18+1, 22+1) => EdgeType.HAD,
    (19+1, 22+1) => EdgeType.HAD,
    (19+1, 28+1) => EdgeType.HAD,
    (20+1, 22+1) => EdgeType.HAD,
    (20+1, 37+1) => EdgeType.HAD,
    (21+1, 22+1) => EdgeType.HAD,
    (21+1, 34+1) => EdgeType.HAD,
    (22+1, 27+1) => EdgeType.HAD,
    (22+1, 38+1) => EdgeType.HAD,
    (22+1, 42+1) => EdgeType.HAD,
    (22+1, 44+1) => EdgeType.HAD,
    (22+1, 46+1) => EdgeType.HAD,
    (23+1, 48+1) => EdgeType.HAD,
    (24+1, 49+1) => EdgeType.HAD,
    (25+1, 50+1) => EdgeType.HAD,
    (26+1, 51+1) => EdgeType.HAD,
    (38+1, 39+1) => EdgeType.HAD,
    (40+1, 41+1) => EdgeType.HAD,
    (42+1, 43+1) => EdgeType.HAD,
    (44+1, 45+1) => EdgeType.HAD,
    (46+1, 47+1) => EdgeType.HAD,
)

using Graphs
zxg = ZXGraph(ZXDiagram(0))
vs = 1:52
for v in vs
    ZXCalculus.add_spider!(zxg, st[v], ZXCalculus.Phase(ps[v]))
end
for (e, _) in es
    Graphs.add_edge!(zxg, e[1], e[2])
end
for i = 1:5
    push!(zxg.inputs, i)
    push!(zxg.outputs, i+23)
end

ZXCalculus.ancilla_extraction(zxg)


text:
module ParticleSimulation

using StaticArrays
using Printf
using LinearAlgebra: norm

export Vec2D, random_point
export md, force_pair, wrap

function wrap(x,side)
    x = mod(x,side)
    if x >= side/2
        x -= side
    elseif x < -side/2
        x += side
    end
    return x
end

function force_pair(x::T,y::T,cutoff,side) where T
    Δv = wrap.(y - x, side)
    d = norm(Δv)
    if d > cutoff
        return zero(T)
    else
        return 2*(d - cutoff)*(Δv/d)
    end
end

function forces!(f::Vector{T},x::Vector{T},force_pair::F) where {T,F}
    fill!(f,zero(T))
    n = length(x)
    for i in 1:n-1
        @inbounds for j in i+1:n
            fpair = force_pair(i,j,x[i],x[j])
            f[i] += fpair
            f[j] -= fpair
        end
    end
    return f
end

function md(x0,v0,mass,dt,nsteps,isave,force_pair::F) where F
    x = copy(x0)
    v = copy(v0)
    f = similar(x0)
    a = similar(x0)
    trajectory = [ copy(x0) ] # will store the trajectory
    for step in 1:nsteps
        # Compute forces
        forces!(f,x,force_pair)
        # Accelerations
        @. a = f / mass
        # Update positions
        @. x = x + v*dt + a*dt^2/2
        # Update velocities
        @. v = v + a*dt
        # Save if required
        if mod(step,isave) == 0
            push!(trajectory,copy(x))
        end
    end
    return trajectory
end

struct Vec2D{T} <: FieldVector{2,T}
    x::T
    y::T
end

function random_point(::Type{Vec2D{T}},range) where T
    p = Vec2D(
        range[begin] + rand(T)*(range[end]-range[begin]),
        range[begin] + rand(T)*(range[end]-range[begin])
    )
    return p
end

end # module

using .ParticleSimulation

function main()
    n = 100
    cutoff = 5.
    side = 100.
    nsteps = 200_000
    trajectory = md((
        x0 = [random_point(Vec2D{Float64},(-50,50)) for _ in 1:n ],
        v0 = [random_point(Vec2D{Float64},(-1,1)) for _ in 1:n ],
        mass = [ 10.0 for _ in 1:100 ],
        dt = 0.001,
        nsteps = nsteps,
        isave = 1000,
        force_pair = (i,j,p1,p2) -> force_pair(p1,p2,cutoff,side)
    )...)
    file = open("traj_julia.xyz","w")
    for (step,x) in pairs(trajectory)
        println(file,n)
        println(file," step = ", step)
        for p in x
            println(file,"He $(wrap(p.x,side)) $(wrap(p.y,side)) 0.")
        end
    end
    close(file)
end

main()


text:
using JSON

default_endpoints = Dict(
  "base" => "https://plot.ly",
  "api" => "https://api.plot.ly/v2")

type PlotlyCredentials
    username::String
    api_key::String
end

type PlotlyConfig
    plotly_domain::String
    plotly_api_domain::String
end

function signin(username::String, api_key::String, endpoints=nothing)
# Define session credentials/endpoint configuration, where endpoint is a Dict

    global plotlycredentials = PlotlyCredentials(username, api_key)

    # if endpoints are specified both the base and api domains must be specified
    if endpoints != nothing
        try
            base_domain = endpoints["plotly_domain"]
            api_domain = endpoints["plotly_api_domain"]
            global plotlyconfig = PlotlyConfig(base_domain, api_domain)
        catch
            error("You must specify both the base and api endpoints.")
        end
    end
end

function get_credentials()
# Return the session credentials if defined --> otherwise use .credentials specs

    if !isdefined(Plotly, :plotlycredentials)

        creds = get_credentials_file()

        try
            username = creds["username"]
            api_key = creds["api_key"]

            global plotlycredentials = PlotlyCredentials(username, api_key)

        catch
            error("Please 'signin(username, api_key)' before proceeding. See
            http://plot.ly/API for help!")
        end
    end

    # will persist for the remainder of the session
    return plotlycredentials
end

function get_config()
# Return the session configuration if defined --> otherwise use .config specs

    if !isdefined(Plotly,:plotlyconfig)

        config = get_config_file()

        if isempty(config)
            base_domain = default_endpoints["base"]
            api_domain = default_endpoints["api"]
        else
            base_domain = get(config, "plotly_domain", default_endpoints["base"])
            api_domain = get(config, "plotly_api_domain", default_endpoints["api"])
        end

        global plotlyconfig = PlotlyConfig(base_domain, api_domain)

    end

    # will persist for the remainder of the session
    return plotlyconfig
end

function set_credentials_file(input_creds::Dict)
# Save Plotly endpoint configuration as JSON key-value pairs in
# userhome/.plotly/.credentials. This includes username and api_key.

    # plotly credentials file
    userhome = homedir()
    plotly_credentials_folder = joinpath(userhome, ".plotly")
    plotly_credentials_file = joinpath(plotly_credentials_folder, ".credentials")

    #check to see if dir/file exists --> if not, create it
    try
        mkdir(plotly_credentials_folder)
    catch err
        isa(err, SystemError) || rethrow(err)
    end

    prev_creds = get_credentials_file()

    #merge input creds with prev creds
    if !isempty(prev_creds)
        creds = merge(prev_creds, input_creds)
    else
        creds = input_creds
    end

    #write the json strings to the cred file
    creds_file = open(plotly_credentials_file, "w")
    write(creds_file, JSON.json(creds))
    close(creds_file)
end

function set_config_file(input_config::Dict)
# Save Plotly endpoint configuration as JSON key-value pairs in
# userhome/.plotly/.config. This includes the plotly_domain, and plotly_api_domain.

    # plotly configuration file
    userhome = homedir()
    plotly_config_folder = joinpath(userhome, ".plotly")
    plotly_config_file = joinpath(plotly_config_folder, ".config")

    #check to see if dir/file exists --> if not create it
    try
        mkdir(plotly_config_folder)
    catch err
        isa(err, SystemError) || rethrow(err)
    end

    prev_config = get_config_file()

    #merge input config with prev config
    if !isempty(prev_config)
        config = merge(prev_config, input_config)
    else
        config = input_config
    end

    #write the json strings to the config file
    config_file = open(plotly_config_file, "w")
    write(config_file, JSON.json(config))
    close(config_file)
end

function get_credentials_file()
# Load user credentials as a Dict

    # plotly credentials file
    userhome = homedir()
    plotly_credentials_folder = joinpath(userhome, ".plotly")
    plotly_credentials_file = joinpath(plotly_credentials_folder, ".credentials")

    if !isfile(plotly_credentials_file)
        creds = Dict()
    else
        creds_file = open(plotly_credentials_file)
        creds = JSON.parse(creds_file)

        if creds == nothing
            creds = Dict()
        end

    end

    return creds

end

function get_config_file()
# Load endpoint configuration as a Dict

    # plotly configuration file
    userhome = homedir()
    plotly_config_folder = joinpath(userhome, ".plotly")
    plotly_config_file = joinpath(plotly_config_folder, ".config")

    if !isfile(plotly_config_file)
        config = Dict()
    else
        config_file = open(plotly_config_file)
        config = JSON.parse(config_file)

        if config == nothing
            config = Dict()
        end

    end

    return config
end