text:
abstract type SCM end

"""
    SynthControlModel

Takes in a `TreatmentPanel` and holds the results of calling the `fit!` method, optimal
weights for comparator units and the resulting predicted outcome of the synthetic control
unit as well as the estimated treatment effect from comparing synthetic control and observed
outcome of the treated unit.

"""

mutable struct SynthControlModel <: SCM
    treatment_panel::TreatmentPanel
    w::Vector{Float64} # Vector of weights
    ŷ₁::Vector{Float64} # Vector of post-treatment predicted outcomes for synthetic unit
    τ̂::Vector{Float64} # Vector of estimated treatmet impacts
    p_test_res::Array{Float64, 2} # Matrix of outcomes for untreated under placebo
end

function SynthControlModel(tp::TreatmentPanel)
    @unpack y₁₀, y₁₁, J, T₁ = tp

    ŷ₁ = zeros(length(y₁₀) + length(y₁₁))
    w = zeros(J)
    τ̂ = zeros(T₁)
    p_test_res = zeros(length(ŷ₁), J)

    SynthControlModel(tp, w, ŷ₁, τ̂, p_test_res)
end

"""
isfitted(s::SynthControlModel)

Check whether the `SynthControlModel` object `s` has been fitted.

# Examples
```julia-repl
julia> df = load_brexit();

julia> s = SynthControlModel(df, :country, :dateid, :realgdp, 86, "United Kingdom");

julia> isfitted(s)
false

julia> fit!(s);

julia> isfitted(s)
true

```
"""
isfitted(s::SynthControlModel) = sum(s.w) ≈ 0.0 ? false : true

"""
    fit!(s::SynthControlModel; placebo_test = false)

Fit the `SynthControlModel` `s` by finding the weights that minimize the distance between
the pre-treatment outcomes for the observational unit of interest and the weighted average
pre-treatment outcomes for unweighted units.

If `placebo_test = true` is supplied, additional placebo tests will be performed by using
every non-treated unit in the data set as the treated unit in turn and estimating the
treatment impact on this unit. Results are stored in the `p_test_res` field and can be
used as the basis for inference.
"""
function fit!(s::SynthControlModel; placebo_test = false)

    isnothing(s.treatment_panel.predictors) || "Predictors currently not implemented"

    @unpack y₁₀, y₁₁, yⱼ₀, yⱼ₁, J, T₀ = s.treatment_panel

    # Objective function
    # Hat tip to Mathieu Tanneau who suggested switching to a quadratic formulation
    obj(w) = dot(y₁₀ .- vec(w'*yⱼ₀), y₁₀ .- vec(w'*yⱼ₀)) + 1e6*(1.0 - sum(w))^2

    # Initial condition and bounds
    initial = [1/J for _ in 1:J]
    lower = zeros(J)
    upper = ones(J)

    # Get weights through optimization
    res = optimize(obj, lower, upper, initial)
    s.w = res.minimizer

    # Get estimates
    s.ŷ₁ = vec(s.w'*[yⱼ₀ yⱼ₁])
    s.τ̂ = ([y₁₀; y₁₁] .- s.ŷ₁)[T₀+1:end]

    if placebo_test
        placebos = unique(s.treatment_panel.data[(s.treatment_panel.data[!,s.pid] .!= s.treat_id), s.pid])
        p = deepcopy(s)

        for n ∈ 1:p.J
            # change treatment ID
            p.treat_id = placebos[n]
            # change outcome data
            p.y = p.data[p.data[!, p.pid] .== p.treat_id, p.outcome]
            p.comps = collect(reshape(p.data[(p.data[!,p.pid] .!= p.treat_id) .& (p.data[!,p.tid].<p.treat_t), p.outcome],
                    p.J, p.no_pretreat_t)')
            fit!(p)
            s.p_test_res[:, n] = p.y .- p.ŷ₁
        end
    end
    return s
end

# Pretty printing
import Base.show

function show(io::IO, ::MIME"text/plain", s::SynthControlModel)

    println(io, "\nSynthetic Control Model")
    println(io, "\tOutcome variable: ", s.treatment_panel.outcome)
    println(io, "\tTime dimension: ",string(s.treatment_panel.t_var)," with ",
            string(length(unique(s.treatment_panel.data[!,s.treatment_panel.t_var]))), " unique values")
    println(io, "\tTreatment period: ", string(s.treatment_panel.treatment[2]))
    println(io, "\tID variable: ",string(s.treatment_panel.id_var), " with ",
            string(length(unique(s.treatment_panel.data[!,s.treatment_panel.id_var]))), " unique values")
    println(io, "\tTreatment ID: ",string(s.treatment_panel.treatment[1]))

    if isfitted(s)
      println(io, "\tModel is fitted")
      println(io, "\tImpact estimates: ",round.(s.τ̂, digits=3))
    else
      println(io, "\n\tModel is not fitted")
    end
end


text:
mutable struct Scalar{T}
    type::Type
    pointer::Ptr

    function Scalar{T}(p::Ptr) where T
        if !haskey(TYPE_MAP, T)
            error("Type $T is not supported.")
        end
        ret = new(T, p)
        finalizer(ret) do s
            ccall((:scalar_destroy, libtorch_capi),
                  Cvoid, (Ptr{Cvoid},),
                  s.pointer)
        end
        return ret
    end

end

function Scalar{T}(s::U) where {T<:TorchNumber, U<:TorchNumber}
    if !haskey(TYPE_MAP, T)
        error("Type $T is not supported.")
    end
    data = T[convert(T, s)]
    ptr = ccall((:scalar_from_data, libtorch_capi),
                Ptr{Cvoid}, (Ptr{Cvoid}, Cchar),
                data, TYPE_MAP[T])
    Scalar{T}(ptr)
end

Scalar(s::T) where T = Scalar{T}(s)

function value(s::Scalar{T}) where T
    data = T[zero(T)]
    ccall((:scalar_value, libtorch_capi),
          Cvoid, (Ptr{Cvoid}, Cchar, Ptr{Cvoid}),
          s.pointer, TYPE_MAP[T], data)
    return data[1]
end

Base.getindex(s::Scalar) = value(s)

function Base.show(io::IO, s::Scalar{T}) where {T}
    write(io, "PyTorch.Scalar{$T} = $(value(s))\n")
end

function Base.display(s::Scalar)
    show(stdout, s)
end

ext: jl


text:
"""
prototyping tools to help with boxy & friends
"""

# using Images, FileIO           # load elevation from PNG
using FileIO
using Interpolations
using Gadfly, Colors
using Optim
using Test
using JSON2
import Caesar

prjPath = dirname(dirname(pathof(Caesar)))
include( joinpath(prjPath, "examples","dev","scalar","Sequences.jl") )

# some Gadfly stuff
# Gadfly.set_default_plot_size(35cm,25cm)
# st=style( major_label_font="CMU Serif",
#           minor_label_font="CMU Serif",
#           major_label_font_size=20pt,
#           minor_label_font_size=20pt)

## Load unique test data

datapath = joinpath(prjPath, "examples", "dev", "scalar", "data")

fid = open(joinpath(datapath, "seq_data_a.json"), "r")
myData_a = JSON2.read(fid, Dict{Symbol, Vector{Float64}}); close(fid)
fid = open(joinpath(datapath, "seq_data_b.json"), "r")
myData_b = JSON2.read(fid, Dict{Symbol, Vector{Float64}}); close(fid)

xseq_a = myData_a[:x_seq]
zseq_a = myData_a[:z_seq]

xseq_b = myData_b[:x_seq]
zseq_b = myData_b[:z_seq]

s_a  = Sequences.MeasurementSequence(myData_a[:x_seq], myData_a[:z_seq])
s_b  = Sequences.MeasurementSequence(myData_b[:x_seq], myData_b[:z_seq])

## Test 1. Globally-referenced sequences on correct grid (zero displacement)

@test Sequences.ssd(s_a,s_b) < 1e-3 #FIXME replace with test value

# find relative position through minimizer
f(x) = Sequences.ssd(s_a, Sequences.displace(s_b, x))
r = Optim.optimize(f, -10.0, 10.0)
@test abs(r.minimizer) < 0.1 # valid for low noise

## Test 2. Locally-referenced
s_a  = Sequences.MeasurementSequence(myData_a[:x_seq], myData_a[:z_seq])
s_b  = Sequences.MeasurementSequence(myData_a[:x_seq], myData_b[:z_seq])
f(x) = Sequences.ssd(s_a, Sequences.displace(s_b, x))
r = Optim.optimize(f, -10.0, 10.0)
@test r.minimizer - 5  < 0.1 # valid for low noise

Gadfly.set_default_plot_size(35cm, 25cm) # make smaller for bigger font
st=style( major_label_font="CMU Serif",
          minor_label_font="CMU Serif",
          major_label_font_size=20pt,
          minor_label_font_size=20pt)
# plot input
p = plot(st)
push!(p, layer(x=s_a.x, y=s_a.y, Geom.line,color=[colorant"orange"]))
push!(p, layer(x=s_b.x, y=s_b.y, Geom.line,color=[colorant"hotpink"]))

# plot output/aligned sequences
p = plot(st)
push!(p, layer(x=s_a.x, y=s_a.y,Geom.line,color=[colorant"orange"]))
s3 = Sequences.displace(s_b, r.minimizer)
push!(p, layer(x=s3.x, y=s3.y,Geom.line,color=[colorant"hotpink"]))

# plot intensity (correlator output)
qr = collect(LinRange(-10.0,10.0,1001))  # displacement query range
v = f.(qr)
k = abs(maximum(v))
pv = v./k
pv = exp.(-pv)
pv ./= sum(pv)
plot(x=qr,y=pv,Geom.line) #,Scale.y_log)
qr[argmax(pv)] # most likely displacement

Sequences.ssd(s1, s2) # should be low

p = plot(st)
push!(p, layer(x=s1.x, y=s1.y,Geom.line,color=[colorant"orange"]))
push!(p, layer(x=s2.x, y=s2.y,Geom.line,color=[colorant"hotpink"]))

f(x) = Sequences.ssd(s1, Sequences.displace(s2, x))
r = optimize(f, -10.0, 10.0)

@test abs(r.minimizer - 5.0)<0.1

# show results
p = plot(st)
push!(p, layer(x=s1.x, y=s1.y,Geom.line,color=[colorant"orange"]))
s3 = Sequences.displace(s2, r.minimizer)
push!(p, layer(x=s3.x, y=s3.y,Geom.line,color=[colorant"hotpink"]))

## 3. slightly more practical use - map ssd to probability

s1 = Sequences.MeasurementSequence(collect(LinRange(0,15,n)), terrain[1:n]+sigma_y*randn(n))
s2 = Sequences.MeasurementSequence(collect(LinRange(0,15,n)), terrain[1+Int(n/3):1+Int(n/3)+n-1]+sigma_y*randn(n))

p = plot(st)
push!(p, layer(x=s1.x, y=s1.y, Geom.line,color=[colorant"orange"]))
push!(p, layer(x=s2.x, y=s2.y, Geom.line,color=[colorant"hotpink"]))

f(x) = Sequences.ssd(s1, Sequences.displace(s2, x))
qr = collect(LinRange(-10.0,10.0,10001))
v = f.(qr)

k = abs(maximum(v))
pv = v./k
pv = exp.(-pv)
pv ./= sum(pv)
plot(x=qr,y=pv,Geom.line,Scale.y_log)
qr[argmax(pv)]


text:
using Documenter
using SVDD

makedocs(
    sitename = "SVDD Documentation",
    pages = [
        "Home" => "index.md",
        "Getting Started" => "start.md",
        "SVDD" => "svdd.md",
        "Custom Solvers" => [
            "SMO" => "smo.md"
            ]
        ],
    format = Documenter.HTML(),
    modules = [SVDD]
)
deploydocs(
    repo = "github.com/englhardt/SVDD.jl.git"
)


text:
# Copyright (c) 2021
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using LinearAlgebra;

"""
    L0EM(A::Matrix{Float64}, b::Vector{Float64}; maxiter=100, epsilon=0.01)

Find the solution to Ax=b using an efficient EM algoritm that directly solves
the L0 optimization problem.

### Input

- `A`       -- Matrix: Ax=b
- `b`       -- Vector: Ax=b
- `maxiter` -- (optional) number of optmization iterations
- `epsilon` -- (optional) threshold to stop optimizing
- `lambda`  -- (optional) regularization

### Output

Solution to Ax=b (Vector{Float64})

### Example

julia> A = randn(10, 100);
julia> b = randn(10);
julia> x = L0EM(A, b, maxiter=5, epsilon=.01);
julia> println(x)

### Algorithm

This function implements Liu and Li's L0EM algorithms in https://arxiv.org/pdf/1407.7508v1.pdf.

"""
function L0EM(A::Matrix{Float64},
              b::Vector{Float64};
              lambda=.001,
              epsilon=.001,
              maxiter=50)
    local n, p, theta, A_eta, eta;
    eps_stop = .01;
    eps_zero = .01;
    n, p = size(A);
    eye = Matrix{Float64}(I, p, p);

    # get the initial solution
    theta = inv(A'*A + lambda*eye)*A'*b;

    # continue to optimize theta
    for i = 1:maxiter
        eta = theta;

        # copy the squared eta terms into a matrix so we can use the hammard product.
        A_eta = repeat(eta.^2, 1, n)';
        A_eta = A_eta.*A;

        # update theta
        theta = inv(A_eta'*A + lambda*eye)*A_eta'*b;

        # check to break the loop due to a small difference in norm
        if norm(theta-eta, 2) <= epsilon
            break;
        end
    end
    x = theta;
    i = abs.(x) .< epsilon;
    x[i] = zeros(sum(i));
    return x;
end

# A = randn(10, 50); b = randn(10);
# x = L0EM(A, b);
# println(x)


text:
using Base.Test
using AperturePhotometry

data = ones(10, 10)

r = 3.0
@test_approx_eq sum_circle(data, 5.0, 5.0, r) pi*r^2

println("tests passed")


text:
"""
`HyperRank.jl`
============

This file contains implementations of [motif-based PageRank](http://www.cse.ust.hk/~yqsong/papers/2018-AAAI-MotifPageRank.pdf)
and the novel Higher-Order HyperRank. Higher-Order HyperRank is a generalization
of motif-based PageRank to hypergraphs, where hypergraph motifs are defined as in
https://arxiv.org/pdf/2003.01853.pdf.

The high-level idea is to perform a hypergraph PageRank in which the random walker
is more likely to traverse hyperedges that take part in more hypergraph motif
instances. However, in practice, the problem is equivalent to performing PageRank
on the projected graph of the hypergraph, where edges between vertices sharing
hyperedges that take part in hypergraph motifs are more likely to be traversed.
In other words, the probability that the random walker goes to a node with which
it shares many hypergraph motifs is increased. Thus, Higher-Order HyperRank can
be computed as a motif-based PageRank with a slightly different motif co-occurrence
matrix definition.
"""

include("Hypergraph.jl") # Hypergraph data structur
include("h_motifs.jl") # Hypergraph-motif finding, classification

using MatrixNetworks # Used for pagerank, graph data structure
using SparseArrays # Used for converting MatrixNetworks and matrices for PageRanking
using Profile # Used for performance measurement

"""
`motif_pagerank`
================

Perform [motif-based PageRank](http://www.cse.ust.hk/~yqsong/papers/2018-AAAI-MotifPageRank.pdf)
on a network.

Arguments
---------
    - `G::Matrix{T}`: The adjacency matrix of the network
    - `M::Matrix{T}`: The motif co-occurrence matrix of the network
    - `d::Float64(=0.5)`: The balance of the weighted average between `G` and `M`.
                          Must be in [0,1]. `d=1.0` means the PageRanked matrix
                          is G, and `d=0.0` means the PageRanked matrix is M.
    - `α::Float64(=0.85)`: The PageRank damping term. Must be in [0,1).
    - `ϵ::Float64(=1.0e-6)`: Tolerance in calulating the PageRank vector

Returns
-------
A vector `x`, where `x[i]` is the motif-based PageRank of node `i`
"""
function motif_pagerank(G::Matrix{T}, M::Matrix{V}; d::Float64=0.5, α::Float64=0.85, ϵ::Float64=1.0e-6) where T <: Real where V <: Real
    if d < 0 || d > 1 throw(ArgumentError("Expected d ∈ [0,1], got $d.")) end # Invalid d
    if α < 0 || α >= 1 throw(ArgumentError("Expected α ∈ [0,1), got $d.")) end # Invalid α

    average = d*(G/maximum(G)) + (1-d)*(M/maximum(M))
    return pagerank(MatrixNetwork(sparse(average)), α, ϵ) # PageRank the weighted average of G and M
end

"""
`hyper_rank`
============

Perform Higher-Order HyperRank on a hypergraph.

Arguments
---------
    - `M::MatrixHypergraph`: The hypergraph to analyze
    - `m::Int8`: The hypergraph motif id, must be in [1,26]. The motifs corresponding
                 to each id can be found on page 3 of the [paper](http://www.cse.ust.hk/~yqsong/papers/2018-AAAI-MotifPageRank.pdf).
    - `d::Float64(=0.5)`: The balance of the weighted average between `G` and `M`.
                          Must be in [0,1]. `d=1.0` means the PageRanked matrix
                          is G, and `d=0.0` means the PageRanked matrix is M.
    - `α::Float64(=0.85)`: The PageRank damping term. Must be in [0,1).
    - `ϵ::Float64(=1.0e-6)`: Tolerance in calculating the PageRank vector.

Returns
-------
A vector `x`, where `x[i]` is the Higher-Order HyperRank of node `i`
"""
function hyper_rank(M::MatrixHypergraph, m::Int8; d::Float64=0.5, α::Float64=0.85, ϵ::Float64=1.0e-6)
    if m < 1 || m > 26 throw(ArgumentError("Expected m ∈ [1,26], got $m.")) end

    G = dyadic_projection(M) # Projection of the hypergraph
    motifs = all_motifs(M) # Find and classify motif instances
    W = motif_cooccurence(M,motifs,m) # Construct the motif co-occurrence matrix
    return motif_pagerank(Matrix(sparse(G)),Matrix(W); d=d, α=α, ϵ=ϵ)
end

"""
`hyper_rank`
============

Perform Higher-Order HyperRank on a hypergraph.

Arguments
---------
    - `M::MatrixHypergraph`: The hypergraph to analyze
    - `motifs::Vector{Vector{Tuple{Int64,Int64,Int64}}}`: The categorization of all triples of hyperedges into hypergraph motifs,
                                                          where motifs[i] contains all triples that are instances of motif `i`
    - `m::Int8`: The hypergraph motif id, must be in [1,26]. The motifs corresponding
                 to each id can be found on page 3 of the [paper](http://www.cse.ust.hk/~yqsong/papers/2018-AAAI-MotifPageRank.pdf).
    - `d::Float64(=0.5)`: The balance of the weighted average between `G` and `M`.
                          Must be in [0,1]. `d=1.0` means the PageRanked matrix
                          is G, and `d=0.0` means the PageRanked matrix is M.
    - `α::Float64(=0.85)`: The PageRank damping term. Must be in [0,1).
    - `ϵ::Float64(=1.0e-6)`: Tolerance in calculating the PageRank vector.

Returns
-------
A vector `x`, where `x[i]` is the Higher-Order HyperRank of node `i`
"""
function hyper_rank(M::MatrixHypergraph, motifs::Vector{Vector{Tuple{Int64,Int64,Int64}}}, m::Int8; d::Float64=0.5, α::Float64=0.85, ϵ::Float64=1.0e-6)
    if m < 1 || m > 26 throw(ArgumentError("Expected m ∈ [1,26], got $m.")) end

    G = dyadic_projection(M) # Projection of the hypergraph
    W = motif_cooccurence(M,motifs,m) # Construct the motif co-occurrence matrix
    return motif_pagerank(Matrix(sparse(G)),Matrix(W); d=d, α=α, ϵ=ϵ)
end

max_stars_repo_name: joshua-matt/HO-HyperRank

max_issues_repo_name: joshua-matt/HO-HyperRank

max_forks_repo_name: joshua-matt/HO-HyperRank


text:
### Fit a full rank model with sketched Frank Wolfe

# todo: make it work for mpca

import FirstOrderOptimization: frank_wolfe_sketched, FrankWolfeParams, DecreasingStepSize,
                               AbstractSketch, AsymmetricSketch,
                               IndexingOperator, LowRankOperator, mysvds
import LowRankModels: fit!, ConvergenceHistory, get_yidxs, grad, evaluate
import Base: axpy!, scale!

export fit_sketch!, FrankWolfeParams

### FITTING
function fit_sketch!(gfrm::GFRM, params::FrankWolfeParams = FrankWolfeParams();
			  ch::ConvergenceHistory=ConvergenceHistory("FrankWolfeGFRM"),
              z::AbstractVector = zeros(sum(map(length, gfrm.observed_examples))),
              sketch::AbstractSketch = AsymmetricSketch(size(gfrm.A)..., gfrm.k),
			  verbose=true,
			  kwargs...)

    if !isa(gfrm.r, TraceNormConstraint)
        error("Frank Wolfe fitting is only implemented for trace norm constrained problems")
    end

    # the functions below close over all this problem data
    yidxs = get_yidxs(gfrm.losses)
    d = maximum(yidxs[end])
    m,n = size(gfrm.A)
    alpha = gfrm.r.scale

    nobsj = map(length, gfrm.observed_examples)
    startobsj = cumsum(vcat(0, nobsj))
    nobs = sum(nobsj)
    # start_obsj = append!([0], cumsum(nobsj)[1:end-1])
    # end_obsj = start_obsj+nobsj
    obs = Array(Int, nobs)
    iobs = 1
    for j=1:n
        for i=gfrm.observed_examples[j]
            obs[iobs] = m*(j-1) + i
            iobs += 1
        end
    end
    indexing_operator = IndexingOperator(m, n, obs)
    @assert size(indexing_operator, 1) == nobs
    @assert length(z) == nobs

    function f(z)
        obj = 0
        for j=1:n
            for (iobs,i) in enumerate(gfrm.observed_examples[j])
                obj += evaluate(gfrm.losses[j], z[startobsj[j]+iobs], gfrm.A[i,j])
            end
        end
        return obj
    end

    ## Grad of f
    function grad_f(z;
                    g = Array(Float64, size(z)))
        iobs = 1
        for j=1:n
            lj = gfrm.losses[j]
            for i=gfrm.observed_examples[j]
                g[iobs] = grad(lj, z[iobs], gfrm.A[i,j])
                iobs += 1
            end
        end
        # return G = A'*g
        return IndexedLowRankOperator(indexing_operator, g)
    end

    const_nucnorm(z) = alpha # we'll always saturate the constraint, don't bother computing it
    # returns solution, optval of min <G, Delta> st ||Delta||_* \leq alpha
    function min_lin_st_nucnorm_sketched(G, alpha)
        ga = Array(G)
        u,s,v = mysvds(ga, nsv=1)
        return LowRankOperator(-alpha*u, v')
    end

    # recover
    t = time()
    X_sketched = frank_wolfe_sketched(
        z,
        f, grad_f,
        const_nucnorm,
        alpha,
        min_lin_st_nucnorm_sketched,
        sketch,
        params,
        ch,
        verbose=verbose
        )

    t = time() - t
    update_ch!(ch, t, f(z))

    gfrm.W = X_sketched

    return X_sketched, ch
end


text:
function test_fft(AT)
    for n = 1:3
        @testset "FFT with ND = $n" begin
            dims = ntuple(i -> 40, n)
            @test compare(fft!, AT, rand(ComplexF32, dims))
            @test compare(ifft!, AT, rand(ComplexF32, dims))

            @test compare(fft, AT, rand(ComplexF32, dims))
            @test compare(ifft, AT, rand(ComplexF32, dims))
        end
    end
end


text:
module LiquidsStructure

### Imports
using SpecialFunctions

### Exports
export AbstractLiquid, AttractiveHardSpheres, CompositePotential,
       DipolarHardSpheres, DipoleDipole, HardCore, HardDisks, HardSpheres, MSA,
       NonLinearSharmaSharma, PercusYevick, RosenfeldFMT, SharmaSharma, SquareWell,
       StructureFactor, VerletWeis, Yukawa, dimensionality, temperature, volume_fraction

### Implementation
include("potentials.jl")
include("liquids.jl")
include("schemes.jl")
include("structure.jl")
include("utils.jl")

end # module


text:
using DrWatson
@quickactivate "EPExperiment"

using EPExperiments

using Random
using Logging
using LoggingExtras
using JLD2
using FileIO
using CSV
using DataFrames
using LinearAlgebra: I
using StatsFuns: logsumexp
using StatsPlots

const LOG_DIR_NAME = "bananas_estimate"
const LOG_FILENAME = "out.log"

@with_kw struct BananasEstimateConfig{T<:AnISConfig}
    experiment_name::String = "test"
    num_replications::Int = 1
    num_samples::Int = 10
    num_betas::Int = 101
    anis_params::T = AnISMHConfig(
        transition_kernel_cfg=MHFunction(
            kernel_fn=kernel_fn
        )
    )
    seed::Int = 1234
end

function banana_reprocess_results(experiment_name)
    result_folder = datadir(LOG_DIR_NAME, experiment_name)
    results = load(joinpath(result_folder, "results.jld2"))

    logger = TeeLogger(
        ConsoleLogger(),
        FileLogger(joinpath(result_folder, LOG_FILENAME))
    )

    banana_process_results(
        results["exp_estimate"],
        results["tabi_results"],
        results["experiment_config"],
        result_folder,
        logger
    )
end

function banana_process_results(
    exp_estimate,
    tabi_results,
    experiment_config,
    result_folder,
    logger
)
    Z1p_chns = tabi_results[:Z1_positive_info]
    Z1n_chns = tabi_results[:Z1_negative_info]
    Z2_chns = tabi_results[:Z2_info]

    gt_file = load(datadir("bananas_ground_truth", "ground_truth.jld2"))
    ground_truth = gt_file["expectation_estimate"]
    @show ground_truth

    ep, _ = convergence_plot(
        "banana",
        Z1p_chns,
        Z2_chns,
        banana_f,
        ground_truth;
        Z1_negative_chn=Z1n_chns,
        ylims=(10^(-6.5),5),
        yticks=[1, 0.01, 0.0001, 0.000001],
        num_xlogspaced=100,
        legend=:bottomleft
    )
    savefig(ep, joinpath(result_folder, "error_plot.pdf"))

    for (name, chn) in pairs(tabi_results)
        with_logger(logger) do
            @info "$name"
            lws = Array(chn[:log_weight])
            for chain_ix in 1:size(lws, 2)
                ess = compute_ess(lws[:,chain_ix])
                @info "Chain $chain_ix: $ess"
            end
        end
    end
end

function bananas_estimate(experiment_config)
    @unpack seed, num_samples, num_replications = experiment_config
    @unpack num_betas, anis_params = experiment_config

    result_folder = mkpath(joinpath(
        datadir(LOG_DIR_NAME),
        savename(experiment_config)
    ))
    logger = TeeLogger(
        ConsoleLogger(),
        FileLogger(joinpath(result_folder, LOG_FILENAME))
    )

    Random.seed!(seed)

    bananas_exp = banana()

    # Make AnIS algorithm
    betas = (geomspace(1, 1001, num_betas) .- 1) ./ 1000
    anis_alg = get_anis_alg(anis_params, betas)
    tabi = TABI(
        TuringAlgorithm(anis_alg, num_samples),
        TuringAlgorithm(anis_alg, num_samples),
        TuringAlgorithm(anis_alg, num_samples)
    )

    # Replicate results N times
    Z1p_chns = Chains[]
    Z1n_chns = Chains[]
    Z2_chns = Chains[]
    exp_ests = []
    for n in 1:num_replications
        println("Run $n:")
        @time begin
        exp_estimate, tabi_results = estimate_expectation(
            bananas_exp,
            tabi;
            store_intermediate_samples=true
        )
        end
        push!(Z1p_chns, tabi_results[:Z1_positive_info])
        push!(Z1n_chns, tabi_results[:Z1_negative_info])
        push!(Z2_chns, tabi_results[:Z2_info])
        push!(exp_ests, exp_estimate)
    end

    Z1p_chns = chainscat(Z1p_chns...)
    Z1n_chns = chainscat(Z1n_chns...)
    Z2_chns = chainscat(Z2_chns...)
    tabi_results = (
        Z1_positive_info = Z1p_chns,
        Z1_negative_info = Z1n_chns,
        Z2_info = Z2_chns
    )
    exp_estimate = mean(exp_ests)

    # Save results and experiment config
    @tagsave(joinpath(result_folder, "results.jld2"), Dict(
        "exp_estimate" => exp_estimate,
        "tabi_results" => tabi_results,
        "experiment_config" => experiment_config
    ))

    banana_process_results(
        exp_estimate,
        tabi_results,
        experiment_config,
        result_folder,
        logger
    )
end

kernel_fn(prior_sample::T, i) where {T<:Real} = Normal(0, 1)
kernel_fn(prior_sample::NamedTuple, i) = map(x -> kernel_fn(x, i), prior_sample)

bananas_estimate(BananasEstimateConfig(
    experiment_name="banana",
    num_samples=1000,
    num_betas=201,
    num_replications=5,
    anis_params=AnISMHConfig(
        transition_kernel_cfg=MHFunction(
            kernel_fn=kernel_fn
        )
    )
))


text:
export Zernike,
       zernike_norm,
       zernike_nm_to_fringe,
       zernike_nm_to_ansi_j,
       zernike_ansi_j_to_nm,
       zernike_noll_to_nm,
       zernike_fringe_to_nm,
       zernike_zero_separation

"""
    Zernike([T, ]n, m, norm=true)

Zernike polynomial of orders `n` and `m`. If `norm` is true the output will be normalized. The output type can be specified with `T`, which will default to `Float64`.

# Examples

```jldoctest
julia> z = Zernike(10, 3)
Zernike{Float64}(n=10, m=3)

julia> z(0.8, 0.2)
-1.6860955965619917

julia> Zernike(BigFloat, 10, 3)(0.8, 0.2)
-1.68609559656199170518675600760616362094879150390625
```
"""
struct Zernike{T<:AbstractFloat,F1,F2}
    n::Int
    m::Int
    basis::F1
    trig_func::F2
    Zernike{T}(n, m, basis, trig_func) where {T} = new{T,typeof(basis),typeof(trig_func)}(n, m, basis, trig_func)
end

function Zernike(T, n, m, norm::Bool=true)
    norm_val = norm ? zernike_norm(n, m) : 1.0
    trig_func = m < 0 ? sin : cos
    basis = x -> norm_val * jacobi((n - m) ÷ 2, 0, abs(m), x)
    Zernike{T}(n, m, basis, trig_func)
end

Zernike(n, m, norm::Bool=true) = Zernike(Float64, n, m, norm)

function Base.show(io::IO, z::Zernike{T}) where T
    print(io, "Zernike{$T}(n=$(z.n), m=$(z.m))")
end

function (z::Zernike{T})(ρ, θ) where T
    x = ρ^2 - 1
    return T(ρ^abs(z.m) * z.basis(x) * z.trig_func(z.m * θ))
end

"""
    kronecker(i,j)

1 if i==j, else 0; mathematical kronecker function
"""
kronecker(i, j) = Int(i == j)

"""
    zernike_norm(n, m)

Norm of Zernike polynomial of radial order n, azimuthal order m.

The norm is the average squared distance to zero.  By multiplying a zernike
value by the norm, the term is given unit stdev or RMS.
"""
function zernike_norm(n, m)
    num = √(2 * (n+1)) / (1 + kronecker(m, 0))
end

"""
    zernike_nm_to_fringe(n, m)

Map (n,m) ANSI indices to a single fringe index.
"""
function zernike_nm_to_fringe(n, m)
    term1 = (1 + (n + abs(m))/2)^2
    term2 = 2*abs(m)
    term3 = (1 + sign(m)) / 2
    return int(term2 - term2 - term3) + 1
end

"""
    zernike_nm_to_ansi_j(n, m)

Map (n,m) ANSI indices to a single ANSI j index.

See also:
    - [`zernike_ansi_j_to_nm`](@ref) (reciprocal of this function)
"""
function zernike_nm_to_ansi_j(n, m)
    return int((n * (n + 2) + m) / 2)
end

"""
    zernike_ansi_to_ansi_j(n, m)

Map (n,m) ANSI indices to a single ANSI j index.

See also:
    - [`zernike_nm_to_ansi_j`](@ref) (reciprocal of this function)
"""
function zernike_ansi_j_to_nm(j)
    n = int(ceil((-3 + √(9 + 8j))/2))
    m = 2j - n * (n + 2)
    return n, m
end

"""
    zernike_noll_to_nm(j)

Map j Noll index to ANSI (n,m) indices.
"""
function zernike_noll_to_nm(j)
    n = int(ceil((-1 + √(1 + 8j))/2) - 1)
    if n == 0
        m = 0
    else
        nseries = int((n+1) * (n+2) / 2)
        residual = j - nseries - 1

        if isodd(j)
            sign = -1
        else
            sign = 1
        end

        if isodd(n)
            ms = [1,1]
        else
            ms = [0]
        end

        for i=0:n÷2
            push!(ms, ms[end]+2)
            push!(ms, ms[end])
        end

        m = ms[residual] * sign
    end
    return n, m
end

"""
    zernike_fringe_to_nm(j)

Map j Fringe index to ANSI (n,m) indices.
"""
function zernike_fringe_to_nm(j)
    m_n = 2 * ceil(√j - 1)
    g_s = (m_n / 2)^2 + 1
    n = m_n / 2 + floor((j-g_s)/2)
    m = m_n - n * (1 - mod(j-g_s, 2) * 2)
    return int(n), int(m)
end

"""
    zernike_zero_separation(n)

Minimum zero separation of Zernike polynomial of radial order n.  Useful for
computing sample count requirements.
"""
function zernike_zero_separation(n)
    return 1 / n^2
end

"""
    zernike(n, m, ρ, θ[; norm])

Zernike polynomial of radial order n and azimuthal order m, evaluated at the
point (ρ, θ).  No normalization is required of (ρ, θ), though the polynomials
are orthogonal only over the unit disk.

norm is a boolean flag indicating whether the result should be orthonormalized
(scaled to unit RMS) or not.
"""
function zernike(n, m, ρ, θ; norm::Bool=true)
    x = ρ^2 - 1
    n_j = (n - m) / 2
    am = abs(m)
    # α=0, β=|m|
    # there is a second syntax where you have x reversed, 1 - ρ^2,
    # in which ase you swap α and β.  It makes absolutely no difference
    out = jacobi(n_j, 0, am, x)
    if m != 0
        if sign(m) == -1
            f = sin
        else
            f = cos
        end
        out *= (ρ^am * f(m*θ))
    end
	if norm
		out *= zernike_norm(n,m)
	end
    return out
end


text:
export propensity_forwarddiff

function propensity_forwarddiff(f::StandardTimeInvariantPropensity, parameter_count)
    pardiffs = []
    for ip in 1:parameter_count
        push!(pardiffs,
            (x, p) -> ForwardDiff.derivative(z -> f(x,[p[1:ip-1];z;p[ip+1:end]]), p[ip])
        )
    end
    return StandardTimeInvariantPropensityGradient(pardiffs)
end

function propensity_forwarddiff(f::JointTimeVaryingPropensity, parameter_count)
    pardiffs = []
    for ip in 1:parameter_count
        push!(pardiffs,
            (t, x, p) -> ForwardDiff.derivative(z -> f(t, x,[p[1:ip-1];z;p[ip+1:end]]), p[ip])
        )
    end
    return JointTimeVaryingPropensityGradient(pardiffs)
end

function propensity_forwarddiff(f::SeparableTimeVaryingPropensity, parameter_count)
    tfactor_pardiffs = []
    for ip in 1:parameter_count
        push!(tfactor_pardiffs,
            (t, p) -> ForwardDiff.derivative(z -> f.tfactor(t,[p[1:ip-1];z;p[ip+1:end]]), p[ip])
        )
    end
    statefactor_pardiffs = []
    for ip in 1:parameter_count
        push!(statefactor_pardiffs,
            (x, p) -> ForwardDiff.derivative(z -> f.statefactor(x,[p[1:ip-1];z;p[ip+1:end]]), p[ip])
        )
    end
    return SeparableTimeVaryingPropensityGradient(f.tfactor, f.statefactor, tfactor_pardiffs, statefactor_pardiffs)
end


text:
# General case - applies to PrescribedHeatModel or PrescribedWaterModel.
# We can't dispach on boundaries at this point because Prescribed models
# do not have that field.
function soil_boundary_state!(
    nf,
    bctype,
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)

end

function soil_boundary_state!(
    nf,
    bctype,
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)

end

# SoilWaterModel methods

"""
    function soil_boundary_state!(
        nf,
        land::LandModel,
        soil::SoilModel,
        water::SoilWaterModel,
        state⁺::Vars,
        aux⁺::Vars,
        nM,
        state⁻::Vars,
        aux⁻::Vars,
        bctype,
        t,
    )

The Dirichlet-type method for `soil_boundary_state!` for the
`SoilWaterModel`.
"""
function soil_boundary_state!(
    nf,
    bctype,
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    execute_based_on_boundaries!(
        nf,
        land,
        soil,
        water,
        water.boundaries,
        state⁺,
        aux⁺,
        nM,
        state⁻,
        aux⁻,
        bctype,
        t,
    )
end

"""
    function soil_boundary_state!(
        nf,
        land::LandModel,
        soil::SoilModel,
        water::SoilWaterModel,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        bctype,
        t,
    )

The Neumann-type method for `soil_boundary_state!` for the
`SoilWaterModel`.
"""
function soil_boundary_state!(
    nf,
    bctype,
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    execute_based_on_boundaries!(
        nf,
        land,
        soil,
        water,
        water.boundaries,
        state⁺,
        diff⁺,
        aux⁺,
        n̂,
        state⁻,
        diff⁻,
        aux⁻,
        bctype,
        t,
    )
end

"""
    function execute_based_on_boundaries!(
        nf,
        land::LandModel,
        soil::SoilModel,
        water::SoilWaterModel,
        boundaries::SurfaceDrivenWaterBoundaryConditions,
        state⁺::Vars,
        aux⁺::Vars,
        nM,
        state⁻::Vars,
        aux⁻::Vars,
        bctype,
        t,
    )

The Dirichlet-type boundary method for `execute_based_on_boundaries!` for the
`SoilWaterModel` when used with `SurfaceDrivenWaterBoundaryConditions`. As
`SurfaceDrivenWaterBoundaryConditions` are Neumann, this does nothing.
"""
function execute_based_on_boundaries!(
    nf,
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    boundaries::SurfaceDrivenWaterBoundaryConditions,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
)
    nothing
end

"""
    function execute_based_on_boundaries!(
        nf,
        land::LandModel,
        soil::SoilModel,
        water::SoilWaterModel,
        boundaries::SurfaceDrivenWaterBoundaryConditions,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        bctype,
        t,
    )

The Neumann-type boundary method for `execute_based_on_boundaries!` for the
`SoilWaterModel` when used with `SurfaceDrivenWaterBoundaryConditions`. This applies
zero flux at the bottom of the domain, and a physical flux based on evaporation,
precipitation, and runoff at the top.
"""
function execute_based_on_boundaries!(
    nf,
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    boundaries::SurfaceDrivenWaterBoundaryConditions,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
)
    precip_model = boundaries.precip_model
    runoff_model = boundaries.runoff_model
    #compute surface flux
    net_surface_flux =
        compute_surface_flux(runoff_model, precip_model, state⁻, t)
    #top
    if bctype == 2
        diff⁺.soil.water.K∇h = n̂ * (-net_surface_flux)
        #bottom
    elseif bctype == 1
        diff⁺.soil.water.K∇h = n̂ * eltype(state⁻)(0.0)
    end

end

"""
    function execute_based_on_boundaries!(
        nf,
        land::LandModel,
        soil::SoilModel,
        water::SoilWaterModel,
        boundaries::GeneralBoundaryConditions,
        state⁺::Vars,
        aux⁺::Vars,
        nM,
        state⁻::Vars,
        aux⁻::Vars,
        bctype,
        t,
    )

The Dirichlet-type boundary method for `execute_based_on_boundaries!` for the
`SoilWaterModel` when used with `GeneralBoundaryConditions`. This applies the
user-supplied functions of space and time as Dirichlet conditions on `ϑ_l`.

Note that not supplying a function (so that it is `nothing`) results in no
boundary condition of this type applied.
"""
function execute_based_on_boundaries!(
    nf,
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    boundaries::GeneralBoundaryConditions,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
)

    bc = boundaries.dirichlet_bc
    if bctype == 2
        if bc.surface_state != nothing
            state⁺.soil.water.ϑ_l = bc.surface_state(aux⁻, t)
        else
            nothing
        end
    elseif bctype == 1
        if bc.bottom_state != nothing
            state⁺.soil.water.ϑ_l = bc.bottom_state(aux⁻, t)
        else
            nothing
        end
    end
end

"""
    function execute_based_on_boundaries!(
        nf,
        land::LandModel,
        soil::SoilModel,
        water::SoilWaterModel,
        boundaries::GeneralBoundaryConditions,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        bctype,
        t,
    )

The Neumann-type boundary method for `execute_based_on_boundaries!` for the
`SoilWaterModel` when used with `GeneralBoundaryConditions`. This applies the
user-supplied functions of space and time as Neumann conditions on the flux,
equal to -K∇h.

Note that not supplying a function (so that it is `nothing`) results in no
boundary condition of this type applied.
"""
function execute_based_on_boundaries!(
    nf,
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    boundaries::GeneralBoundaryConditions,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
)

    bc = boundaries.neumann_bc
    if bctype == 2
        if bc.surface_flux != nothing
            # Note that -K∇h is the flux, so we need a minus sign here.
            diff⁺.soil.water.K∇h = n̂ * (-bc.surface_flux(aux⁻, t))
        else
            nothing
        end
    elseif bctype == 1
        if bc.bottom_flux != nothing
            # Note that -K∇h is the flux, so we add in a minus sign
            # we add a second minus sign because we choose to specify a BC
            # in terms of ẑ, not n̂, as then it doesnt change directions.
            diff⁺.soil.water.K∇h = -n̂ * (-bc.bottom_flux(aux⁻, t))
        else
            nothing
        end
    end
end

# SoilHeatModel - GeneralBoundaryConditions only so far

"""
    function soil_boundary_state!(
        nf,
        bctype,
        land::LandModel,
        soil::SoilModel,
        heat::SoilHeatModel,
        state⁺::Vars,
        aux⁺::Vars,
        nM,
        state⁻::Vars,
        aux⁻::Vars,
        t,
    )

The Dirichlet-type method for `soil_boundary_state!` for the
`SoilHeatModel`.
"""
function soil_boundary_state!(
    nf,
    bctype,
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    execute_based_on_boundaries!(
        nf,
        land,
        soil,
        heat,
        heat.boundaries,
        state⁺,
        aux⁺,
        nM,
        state⁻,
        aux⁻,
        bctype,
        t,
    )
end

"""
    function soil_boundary_state!(
        nf,
        bctype,
        land::LandModel,
        soil::SoilModel,
        heat::SoilHeatModel,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        t,
    )

The Neumann-type method for `soil_boundary_state!` for the
`SoilHeatModel`.
"""
function soil_boundary_state!(
    nf,
    bctype,
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    execute_based_on_boundaries!(
        nf,
        land,
        soil,
        heat,
        heat.boundaries,
        state⁺,
        diff⁺,
        aux⁺,
        n̂,
        state⁻,
        diff⁻,
        aux⁻,
        bctype,
        t,
    )
end

"""
    function execute_based_on_boundaries!(
        nf,
        land::LandModel,
        soil::SoilModel,
        heat::SoilHeatModel,
        boundaries::GeneralBoundaryConditions,
        state⁺::Vars,
        aux⁺::Vars,
        nM,
        state⁻::Vars,
        aux⁻::Vars,
        bctype,
        t,
    )

The Dirichlet-type boundary method for `execute_based_on_boundaries!` for the
`SoilHeatModel` when used with `GeneralBoundaryConditions`. This applies the
user-supplied functions of space and time as Dirichlet conditions on `ρe_int`.

Note that not supplying a function (so that it is `nothing`) results in no
boundary condition of this type applied. Additionally, note that the user
supplies a Dirichlet condition for temperature, and this is converted into
a boundary condition on volumetric. internal energy.
"""
function execute_based_on_boundaries!(
    nf,
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    boundaries::GeneralBoundaryConditions,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
)
    bc = boundaries.dirichlet_bc
    if bctype == 2
        if bc.surface_state != nothing
            ϑ_l, θ_i = get_water_content(land.soil.water, aux⁻, state⁻, t)
            θ_l = volumetric_liquid_fraction(
                ϑ_l,
                land.soil.param_functions.porosity,
            )
            ρc_s = volumetric_heat_capacity(
                θ_l,
                θ_i,
                land.soil.param_functions.ρc_ds,
                land.param_set,
            )

            ρe_int_bc = volumetric_internal_energy(
                θ_i,
                ρc_s,
                bc.surface_state(aux⁻, t),
                land.param_set,
            )

            state⁺.soil.heat.ρe_int = ρe_int_bc
        else
            nothing
        end
    elseif bctype == 1
        if bc.bottom_state != nothing
            ϑ_l, θ_i = get_water_content(land.soil.water, aux⁻, state⁻, t)
            θ_l = volumetric_liquid_fraction(
                ϑ_l,
                land.soil.param_functions.porosity,
            )
            ρc_s = volumetric_heat_capacity(
                θ_l,
                θ_i,
                land.soil.param_functions.ρc_ds,
                land.param_set,
            )

            ρe_int_bc = volumetric_internal_energy(
                θ_i,
                ρc_s,
                bc.bottom_state(aux⁻, t),
                land.param_set,
            )

            state⁺.soil.heat.ρe_int = ρe_int_bc
        else
            nothing
        end
    end
end

"""
    function execute_based_on_boundaries!(
        nf,
        land::LandModel,
        soil::SoilModel,
        heat::SoilHeatModel,
        boundaries::GeneralBoundaryConditions,
        state⁺::Vars,
        diff⁺::Vars,
        aux⁺::Vars,
        n̂,
        state⁻::Vars,
        diff⁻::Vars,
        aux⁻::Vars,
        bctype,
        t,
    )

The Neumann-type boundary method for `execute_based_on_boundaries!` for the
`SoilHeatModel` when used with `GeneralBoundaryConditions`. This applies the
user-supplied functions of space and time as Neumann conditions on the flux,
equal to -κ∇T.

Note that not supplying a function (so that it is `nothing`) results in no
boundary condition of this type applied.
"""
function execute_based_on_boundaries!(
    nf,
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    boundaries::GeneralBoundaryConditions,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
)
    bc = boundaries.neumann_bc
    if bctype == 2
        if bc.surface_flux != nothing
            # Again, the flux is -κ∇T
            diff⁺.soil.heat.κ∇T = n̂ * (-bc.surface_flux(aux⁻, t))
        else
            nothing
        end
    elseif bctype == 1
        if bc.bottom_flux != nothing
            # two minus signs - the flux is minus κ∇T, and n̂ is -ẑ at the
            # bottom of the domain.
            diff⁺.soil.heat.κ∇T = -n̂ * (-bc.bottom_flux(aux⁻, t))
        else
            nothing
        end
    end
end

max_forks_repo_name: jm-c/CLIMA


text:
import LinearAlgebra.norm
export norm_inf, norm, norm_1

# deprecate these soon
norm_inf(x::AbstractExpr) = maximum(abs(x))
norm_1(x::AbstractExpr) = sum(abs(x))
norm_fro(x::AbstractExpr) = norm2(vec(x))

# behavior of norm should be consistent with julia:
# * vector norms for vectors
# * operator norms for matrices
function norm(x::AbstractExpr, p::Real=2)
    if length(size(x)) <= 1 || minimum(size(x))==1
        # x is a vector
        if p == 1
            return norm_1(x)
        elseif p == 2
            return norm2(x)
        elseif p == Inf
            return norm_inf(x)
        elseif p > 1
            # TODO: allow tolerance in the rationalize step
            return rationalnorm(x, rationalize(Int, float(p)))
        else
            error("vector p-norms not defined for p < 1")
        end
    else
        # TODO: After the deprecation period, allow this again but make it consistent with
        # LinearAlgebra, i.e. make norm(x, p) for x a matrix the same as norm(vec(x), p).
        Base.depwarn("`norm(x, p)` for matrices will in the future be equivalent to " *
                     "`norm(vec(x), p)`. Use `opnorm(x, p)` for the Julia 0.6 behavior of " *
                     "computing the operator norm for matrices.", :norm)
        return opnorm(x, p)
    end
end

if isdefined(LinearAlgebra, :vecnorm) # deprecated but defined
    import LinearAlgebra: vecnorm
end
Base.@deprecate vecnorm(x::AbstractExpr, p::Real=2) norm(vec(x), p)


text:
using ConvexBodyProximityQueries
using CurveProximityQueries
using Plots
import Random: seed!

pgfplots()

include("EllipseBounds.jl")
using .EllipseBounds

function drawellipses!(a, b::Curve, maxnum; atol=1e-8)
    leaves, glb, gub, soln = CurveProximityQueries._init(a, b)
    i = 0
    while (gub - glb) > atol && i < maxnum
        glb, gub, soln = CurveProximityQueries._loop(a, b, leaves, glb, gub)
        i += 1
    end
    for interval in leaves
        β, lb = interval
        if diam(β) > 1e-3
            bell = Ellipse(b, β)
            plot!(bell,
                 linecolor=RGBA(102/255,130/225,223/255,0.5),
                 fillcolor=RGBA(102/255,130/225,223/255,0.5),
                 )
        end
    end
end

seed!(22)
poly = randpoly([-0.5, 0.4]; n=20, scale=0.5)

curve = Bernstein([[0.944223, 0.177264], [0.46193, 0.105887], [0.791169, 0.172386], [0.256546, 0.463735], [0.70966, 0.3173], [0.664713, 0.830956], [0.952688, 0.736242], [0.00763912, 0.97189], [0.270595, 0.827442], [0.16636, 0.996615], [0.381073, 0.543077], [0.225381, 0.754049], [0.570043, 0.638543], [0.254696, 0.988563]])

pts = closest_points(poly, curve)

cases = [nothing, 0, 1, 3, 7, 100]
alph = 'a':'z'
for (i, c) in enumerate(cases)
    plot(poly,
         fillalpha=1.0,
         fillcolor=Gray(0.7),
        )
    isnothing(c) || drawellipses!(poly, curve, c)
    plot!(curve,
          linewidth=2.0,
          linecolor=:black,
         )
    savefig("out/fig6$(alph[i]).pdf")
end
plot!(pts,
      linestyle=:dash,
      markersize=8.0,
      linewidth=2.0,
     )
savefig("out/fig6$(alph[length(cases)]).pdf")


text:
# -- PRIVATE METHODS ---------------------------------------------------------------------------------------------------- #
# ----------------------------------------------------------------------------------------------------------------------- #

# -- PUBLIC METHODS ----------------------------------------------------------------------------------------------------- #

"""
    compute_binary_price_tree(initialPrice::Float64, treeHeight::Int64, deltaFunction::Function)

``\\frac{n!}{k!(n - k)!} = \\binom{n}{k}``

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function compute_binary_price_tree(initialPrice::Float64, treeHeight::Int64, deltaFunction::Function)::PSResult

    # TODO: check - are the args legit?

    # initialize -
    number_of_nodes = Int(2^(treeHeight + 1) - 1)
    priceTree = Array{Float64,1}(undef, number_of_nodes)

    try

        # set the current price -
        priceTree[1] = initialPrice

        # main loop -
        for node_index = 1:number_of_nodes

            # for this node, what are the delta values?
            (u, d) = deltaFunction(treeHeight)

            # what is the current price -
            basePrice = priceTree[node_index]

            # build my kid's prices -
            upValue = u * basePrice
            downValue = d * basePrice

            # left -
            leftIndex = 2 * node_index
            if (leftIndex <= number_of_nodes)
                priceTree[leftIndex] = upValue
            end

            # right
            rightIndex = 2 * node_index + 1
            if (rightIndex <= number_of_nodes)
                priceTree[rightIndex] = downValue
            end
        end

        # return -
        return PSResult{Array{Float64,1}}(priceTree)
    catch error
        return PSResult(error)
    end
end

"""
    compute_ternary_price_tree(initialPrice::Float64, treeHeight::Int64, deltaFunction::Function)

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function compute_ternary_price_tree(initialPrice::Float64, treeHeight::Int64, deltaFunction::Function)::PSResult

    # TODO: check - are the args legit?

    # initialize -
    number_of_nodes = Int((3^(treeHeight + 1) - 1) / 2)
    priceTree = Array{Float64,1}(undef, number_of_nodes)

    try

        # set the current price -
        priceTree[1] = initialPrice

        # main loop -
        for node_index = 1:number_of_nodes

            # for this node, what are the delta values?
            (u, m, d) = deltaFunction(treeHeight)

            # TODO: check, is the u,m,d values ok ...?

            # what is the current price -
            basePrice = priceTree[node_index]

            # build my kid's prices -
            upValue = u * basePrice
            midValue = m * basePrice
            downValue = d * basePrice

            # set my kids -

            # left -
            leftIndex = 3 * node_index - 1
            if (leftIndex <= number_of_nodes)
                priceTree[leftIndex] = upValue
            end

            # mid -
            midIndex = 3 * node_index
            if (midIndex <= number_of_nodes)
                priceTree[midIndex] = midValue
            end

            # right
            rightIndex = 3 * node_index + 1
            if (rightIndex <= number_of_nodes)
                priceTree[rightIndex] = downValue
            end
        end

        # return -
        return PSResult{Array{Float64,1}}(priceTree)
    catch error
        return PSResult(error)
    end
end
# ----------------------------------------------------------------------------------------------------------------------- #


text:
# File holding the random forest object

using Random
using Base.Threads
include("DecisionTree.jl")

mutable struct RFR <: AbstractRegressor

    # internal variables
    n_features::Union{Int, Nothing}
    trees::Vector{DTRegressor}

    # external parameters
    n_trees::Int
    max_depth::Union{Int,Nothing}
    max_features::Union{Int, Nothing}
    min_samples_leaf::Int
    # random_state::Union{AbstractRNG, Int}
    random_state::Int
    bootstrap::Bool
    α::Float64
    param_dict::Dict

    # Initialize random forest object with default parameters
    RFR(;
        n_trees=100,
        max_depth=nothing,
        max_features=nothing,
        min_samples_leaf=1,
        random_state=0,
        bootstrap=true,
        α=0.0,
        param_dict=Dict()) = new(
            nothing, [],
            haskey(param_dict, :n_trees) ? param_dict[:n_trees][1] : n_trees,
            haskey(param_dict, :max_depth) ? param_dict[:max_depth][1] : max_depth,
            haskey(param_dict, :max_features) ? param_dict[:max_features][1] : max_features,
            haskey(param_dict, :min_samples_leaf) ? param_dict[:min_samples_leaf][1] : min_samples_leaf,
            # check_random_state(random_state),
            haskey(param_dict, :random_state) ? param_dict[:random_state][1] : random_state,
            haskey(param_dict, :bootstrap) ? param_dict[:bootstrap][1] : bootstrap,
            haskey(param_dict, :α) ? param_dict[:α][1] : α,
            param_dict)
end

# Fit random forest to data
function fit!(forest::RFR, X::Matrix, Y::Matrix)
    @assert size(Y, 2) == 1 "Y must be 1d"

    # set internal variables
    forest.n_features = size(X, 2)

    # Allocate space
    # trees is a list containing all regression trees
    forest.trees = Array{DTRegressor}(undef, forest.n_trees)

    # Grow the trees
    # Using multi-threading to run over multiple processor cores
    @inbounds @threads for i in 1:forest.n_trees
        # rng_states[i] = copy(forest.random_state)
        forest.trees[i] = create_tree(forest, X, Y)
    end

    return
end

# Create a tree within a forest
function create_tree(forest::RFR, X::Matrix, Y::Matrix)
    n_samples = size(X, 1)

    # Tree random state is random state + ThreadID
    tree_rs = MersenneTwister(Threads.threadid() + forest.random_state)

    # In case of bootstrapping allocate space and choose observations randomly
    inds = Array{Int}(undef, n_samples)
    if forest.bootstrap
        inds = [rand(tree_rs, 1:n_samples) for i in 1:n_samples]
        unique!(inds)
        X_ = copy(X[inds, :])
        Y_ = copy(Y[inds, :])
    else
        X_ = copy(X)
        Y_ = copy(Y)
    end

    # Make tree by invoking DTRegressor object with the parameters of the random forest
    new_tree = DTRegressor(max_depth= forest.max_depth,
                            max_features = forest.max_features,
                            min_samples_leaf = forest.min_samples_leaf,
                            random_state = tree_rs,
                            α = forest.α)

    # This function calls the DTRegressor version of fit!
    fit!(new_tree, X_, Y_)

    return new_tree
end

# Prediction of random forest for a given data set X
function predict(forest::RFR, X::Matrix)
    @assert forest.n_features == size(X, 2) "# of features are not the same"

    # Throw error if tree is not fitted
    if length(forest.trees) == 0
        error("Forest is not fitted")
    end

    # Allocate space for predictions
    prediction = Array{Float64}(undef, (size(X, 1), forest.n_trees))

    # Get predictions from each tree within the forest
    for (ind, t) in enumerate(forest.trees)
        prediction[:,ind] = predict(t, X)
    end

    # Random forest prediction is then mean over all tree predictions
    out_arr = mean(prediction, dims=2)
    return out_arr
end

# Getting the strong variable selection frequency from a random forest
# Same as in a tree but averaged over all trees in a random forest
function strong_selection_freq(forest::RFR, var_index::Int)
    res_arr = Array{Float64}(undef, forest.n_trees)

    for (ind, tree) in enumerate(forest.trees)
        tmp = tree.split_dimensions[tree.split_dimensions .!= nothing]
        res_arr[ind] = sum(tmp .<= var_index)/length(tmp)
    end
    return mean(res_arr)
end

# Return the average depth of trees in a random forest
function average_depth(forest::RFR)

    # Allocate space
    depth_result = Array{Int}(undef, forest.n_trees)

    # Loop through all trees
    for (ind, tree) in enumerate(forest.trees)
        # Taking the maximum of the depth
        depth_result[ind] = maximum(tree.depth_list)
    end
    return mean(depth_result)
end

# Return all P(t_L) and depth - used for Figure 2
function lambda_depth(forest::RFR)

    # Allocate space
    return_arr = []
    depth_arr = []

    # Loop through trees
    for tree in forest.trees
        # Only add non-terminal nodes
        push!(return_arr, tree.pl[tree.pl .!= -2])
        push!(depth_arr, tree.depth_list[tree.pl .!= -2])
    end

    # Flatten array before returning
    return [collect(Iterators.flatten(return_arr)),
            collect(Iterators.flatten(depth_arr))]
end

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
module TestFilter

println("TestFilter: Tests various features of the symbolic handling.")

using ModiaMath: plot

using Modia
using Modia.Electric

@model LPfilter begin
    R = Resistor(R=100.0)
    C = Capacitor(C=2.5E-3, v=Float(start=0.0))
    V = ConstantVoltage(V=10.0)
    ground = Ground()
    @equations begin
        connect(V.n, ground.p)
        connect(V.p, R.p)
        connect(R.n, C.p)
        connect(C.n, V.n)
    end
end
result = simulate(LPfilter, 2)
plot(result, ("R.v", "C.v", "V.v"), figure=3)

checkSimulation(LPfilter, 2, "C.v", 9.996843043981416, logTranslation=true, logSimulation=true, storeEliminated=false)
simulate(LPfilter, 2, aliasElimination=true)
checkSimulation(LPfilter, 2, "C.v", 9.996843043981416, aliasElimination=true, logName="LPfilter aliasElimination")
checkSimulation(LPfilter, 2, "C.v", 9.996843043981416, aliasElimination=true, removeSingularities=true, logName="LPfilter aliasElimination removeSingularities")
simulate(LPfilter, 2.0, aliasElimination=true, removeSingularities=true, logName="LPfilter aliasElimination removeSingularities")
#checkSimulation(LPfilter, 2, "C.v", 9.996843043981416, aliasElimination=true, removeSingularities=true, expandArrayIncidence=true, logName="LPfilter aliasElimination removeSingularities expandArrayIncidence")
#checkSimulation(LPfilter, 2, "C.v", 9.996843043981416, aliasElimination=true, removeSingularities=true, expandArrayIncidence=true, useIncidenceMatrix=true, logName="LPfilter aliasElimination removeSingularities expandArrayIncidence useIncidenceMatrix")
#checkSimulation(LPfilter, 2, "C.v", 9.996843043981416, newStateSelection=true)

@model LPfilterWithoutGround begin
    R = Resistor(R=100.0)
    C = Capacitor(C=2.5E-3)
    V = ConstantVoltage(V=10.0)
    @equations begin
        connect(V.p, R.p)
        connect(R.n, C.p)
        connect(C.n, V.n)
    end
end

#res = simulate(LPfilterWithoutGround , 2)
checkSimulation(LPfilterWithoutGround, 2, "C.v", 9.996843043981416, removeSingularities=true)

@model LPfilterWithoutGround1 begin
    R = Resistor(R=100.0Ohm) # , n=Pin(v=Float(start=0.0)))
    C = Capacitor(C=2.5Milli * Farad, v=Voltage(start=0.0Volt))
    V = ConstantVoltage(V=10.0Volt)
    @equations begin
        connect(V.p, R.p)
        connect(R.n, C.p)
        connect(C.n, V.n)
    end
end
#checkSimulation(LPfilterWithoutGround1, 2, "C.v", 9.996843043981416)

@model LPfilterWithoutGround2 begin
    R = Resistor(R=100) # , n=Pin(v=Float(start=0.0)))
    C = Capacitor(C=2.5E-3, v=Float(start=0.0))
    V = ConstantVoltage(V=10)
    @equations begin
        connect(V.p, R.p)
        connect(R.n, C.p)
        connect(C.n, V.n)
    end
end

@model LPfilterAndSineSource begin
    R = Resistor(R=1)
    C = Capacitor(C=1)
    V = SineVoltage(V=1)
    @equations begin
        connect(V.p, R.p)
        connect(R.n, C.p)
        connect(C.n, V.n)
    end
end
result = simulate(LPfilterAndSineSource, 10)
plot(result, ("C.v", "V.v"), heading="LPfilterAndSineSource", figure=4)

@model LPfilterModel begin
    p = Pin()
    n = Pin()
    R = Resistor(R=1)
    C = Capacitor(C=1)
    ground = Ground()
    @equations begin
        connect(p, R.p)
        connect(n, R.n)
        connect(R.n, C.p)
        connect(C.n, ground.p)
    end
end

# Redeclaration of component models
@model HPfilter begin
    @extends LPfilter(R=Capacitor(C=2), C=Resistor(R=3))
end
checkSimulation(HPfilter, 2, "C.v", 7.16540372163548, removeSingularities=true)

# Redeclaration of models
@model GenericFilter begin
    m = Resistor
    @extends LPfilter(C=m(C=2E-3))
end

@model NewFilter begin
    @extends GenericFilter(m=Capacitor)
end
checkSimulation(NewFilter, 2, "C.v", 9.999596486928553, removeSingularities=true)

#high = true
# Conditional model constructor
@model CondFilter begin
    high = true  # Does not work without this.high.
    @extends LPfilter(R=if high Capacitor(C=2) else Resistor(R=3) end,
    C=if high Resistor(R=1) else Capacitor(C=1) end)
end
checkSimulation(CondFilter, 2, "C.v", 3.678778096777374, removeSingularities=true)

# Conditional model
@model CondFilter2 begin
    high = true
    @extends LPfilter(R=high ? Capacitor(C=2) : Resistor(R=1),
    C=high ? Resistor(R=1) : Capacitor(C=2))
end
checkSimulation(CondFilter2, 2, "C.v", 3.678778096777374, removeSingularities=true)

# Model arrays
Models = [Resistor, Capacitor]
high = true
@model FilterModels begin
    selectModel = 2
    @extends LPfilter(R=high ? Models[selectModel](C=2) : Resistor(R=2), C=high ? Resistor(R=1) : Capacitor(C=2))
    @equations begin
  end
end
checkSimulation(FilterModels, 2, "C.v", 3.678778096777374, removeSingularities=true)

# Model component arrays
ModelComponents = [Resistor(R=5), Capacitor(C=1)]

@model FilterComponents begin
    selectModel = 2
    high = true
    @extends LPfilter(R=high ? ModelComponents[selectModel] : Resistor(),
    C=high ? Resistor(R=2) : Capacitor())
end
checkSimulation(FilterComponents, 2, "C.v", 3.678778096777374)

@model TenCoupledFilters begin
    V = ConstantVoltage(V=1)
    ground = Ground()
    F1 = LPfilterModel()
    F2 = LPfilterModel()
    F3 = LPfilterModel()
    F4 = LPfilterModel()
    F5 = LPfilterModel()
    F6 = LPfilterModel()
    F7 = LPfilterModel()
    F8 = LPfilterModel()
    F9 = LPfilterModel()
    F10 = LPfilterModel()
    @equations begin
        connect(V.n, ground.p)
        connect(V.p, F1.p)
        connect(F1.n, F2.p)
        connect(F2.n, F3.p)
        connect(F3.n, F4.p)
        connect(F4.n, F5.p)
        connect(F5.n, F6.p)
        connect(F6.n, F7.p)
        connect(F7.n, F8.p)
        connect(F8.n, F9.p)
        connect(F9.n, F10.p)
    end
end
checkSimulation(TenCoupledFilters, 2, "F10.C.v", 1.232726022885833e-5)
checkSimulation(TenCoupledFilters, 2, "F10.C.v", 1.232726022885833e-5, aliasElimination=true)

end


text:
using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools

f(u,p,t) = (1.01) * u
f1(u,p,t) = (1.01)/2 * u
f2(u,p,t) = (1.01)/2 * u
σ(u,p,t) = 0.87u
f_split_analytic(u0,p,t,W) = @. u0*exp(0.63155t+0.87W)

prob = SplitSDEProblem{false}(f1,f2,σ,1/2,(0.0,1.0))
sol = solve(prob,SplitEM(),dt=1/10,save_noise=true)

prob = SDEProblem{false}(f,σ,1/2,(0.0,1.0),noise = NoiseWrapper(sol.W))
sol2 = solve(prob,EM(),dt=1/10)

@test sol[:] ≈ sol2[:]

u0 = rand(4)

ff_split = SplitSDEFunction(f1,f2,σ,analytic=f_split_analytic)
prob = SplitSDEProblem(ff_split,σ,u0,(0.0,1.0))

sol = solve(prob,SplitEM(),dt=1/10,save_noise=true)

prob = SDEProblem{false}(f,σ,u0,(0.0,1.0),noise = NoiseWrapper(sol.W))

sol2 = solve(prob,EM(),dt=1/10)

@test sol[end][:] ≈ sol2[end][:]

################################################################################

### Only first

println("Only First")

α = 0.1
β = 0.5
ff1 = (u,p,t) -> β./sqrt.(1+t) - u./(2*(1+t))
ff2 = (u,p,t) -> 0.0
σ2 = (u,p,t) -> α*β./sqrt.(1+t)
ff1_analytic(u0,p,t,W) = @. u0/sqrt(1+t) + β*(t+α*W)/sqrt(1+t)
f_ff1 = SplitSDEFunction(ff1,ff2,σ2,analytic=ff1_analytic)
prob = SplitSDEProblem(f_ff1,σ2,1.,(0.0,1.0))

sol = solve(prob,EM(),dt=1/10)
sol2 = solve(prob,SKenCarp(),dt=1/10)

dts = (1/2) .^(10:-1:2) #14->7 good plot
sim10 = test_convergence(dts,prob,SKenCarp(),trajectories=Int(1e1))
@test abs(sim10.𝒪est[:final]-2) < 0.3

### Only second

println("Only Second")

α = 0.1
β = 0.5
ff1 = (u,p,t) -> 0.0
ff2 = (u,p,t) -> β./sqrt.(1+t) - u./(2*(1+t))
σ2 = (u,p,t) -> α*β./sqrt.(1+t)
f_ff1 = SplitSDEFunction(ff1,ff2,σ2,analytic=ff1_analytic)
prob = SplitSDEProblem(f_ff1,σ2,1.,(0.0,1.0))

sol = solve(prob,EM(),dt=1/10)
sol2 = solve(prob,SKenCarp(),dt=1/10,seed=1)

dts = (1/2) .^ (10:-1:2) #14->7 good plot
sim10 = test_convergence(dts,prob,SKenCarp(),trajectories=Int(1e1))
@test abs(sim10.𝒪est[:final]-2) < 0.3

### Both

println("Both")

α = 0.1
β = 0.5
ff1 = (u,p,t) -> β./sqrt.(1+t)
ff2 = (u,p,t) -> - u./(2*(1+t))
σ2 = (u,p,t) -> α*β./sqrt.(1+t)
f_ff1 = SplitSDEFunction(ff1,ff2,σ2,analytic=ff1_analytic)
prob = SplitSDEProblem(f_ff1,σ2,1.,(0.0,1.0))

sol = solve(prob,EM(),dt=1/10)
sol2 = solve(prob,SKenCarp(),dt=1/10)

dts = (1/2) .^ (10:-1:2) #14->7 good plot
sim10 = test_convergence(dts,prob,SKenCarp(),trajectories=Int(1e1))
@test abs(sim10.𝒪est[:final]-2) < 0.3
sim10 = test_convergence(dts,prob,SKenCarp(nlsolve=StochasticDiffEq.NLFunctional()),trajectories=Int(1e1))
@test abs(sim10.𝒪est[:final]-2) < 0.3

################################################################################

### Only first

println("Only First")

α = 0.1
β = 0.5
ff1 = (du,u,p,t) -> @. du = β/sqrt(1+t) - u/(2*(1+t))
ff2 = (du,u,p,t) -> @. du = 0.0
σ2 = (du,u,p,t) -> @. du = α*β/sqrt(1+t)
f_ff1 = SplitSDEFunction(ff1,ff2,σ2,analytic=ff1_analytic)
prob = SplitSDEProblem(f_ff1,σ2,[1.],(0.0,1.0))

sol = solve(prob,EM(),dt=1/10)
sol2 = solve(prob,SKenCarp(),dt=1/10)

dts = (1/2) .^ (10:-1:2) #14->7 good plot
sim10 = test_convergence(dts,prob,SKenCarp(),trajectories=Int(1e1))
@test abs(sim10.𝒪est[:final]-2) < 0.3

### Only second

println("Only Second")

α = 0.1
β = 0.5
ff1 = (du,u,p,t) -> @. du = 0.0
ff2 = (du,u,p,t) -> @. du = β/sqrt(1+t) - u/(2*(1+t))
σ2 = (du,u,p,t) -> @. du = α*β/sqrt(1+t)
f_ff1 = SplitSDEFunction(ff1,ff2,σ2,analytic=ff1_analytic)
prob = SplitSDEProblem(f_ff1,σ2,[1.],(0.0,1.0))

sol = solve(prob,EM(),dt=1/10)
sol2 = solve(prob,SKenCarp(),dt=1/10)

dts = (1/2) .^ (10:-1:2) #14->7 good plot
sim10 = test_convergence(dts,prob,SKenCarp(),trajectories=Int(1e1))
@test abs(sim10.𝒪est[:final]-2) < 0.3

### Both

println("Both")

α = 0.1
β = 0.5
ff1 = (du,u,p,t) -> @. du = β/sqrt(1+t)
ff2 = (du,u,p,t) -> @. du = - u/(2*(1+t))
σ2 = (du,u,p,t) -> @. du = α*β/sqrt(1+t)
f_ff1 = SplitSDEFunction(ff1,ff2,σ2,analytic=ff1_analytic)
prob = SplitSDEProblem(f_ff1,σ2,[1.],(0.0,1.0))

sol = solve(prob,EM(),dt=1/10)
sol2 = solve(prob,SKenCarp(),dt=1/10)

dts = (1/2) .^ (10:-1:2) #14->7 good plot
sim10 = test_convergence(dts,prob,SKenCarp(),trajectories=Int(1e1))
@test abs(sim10.𝒪est[:final]-2) < 0.3
sim10 = test_convergence(dts,prob,SKenCarp(nlsolve=StochasticDiffEq.NLFunctional()),trajectories=Int(1e1))
@test abs(sim10.𝒪est[:final]-2) < 0.3


text:
#=
Applying SWAP Test on a MPS-reusable circuit to get states' overlap between 2 wave functions.
=#
export MSCpar, MScircuit, MStestRes, MStest, MSTtestRes, MSTtest, nMcheck

"""
    MSCpar(MSCircuit::ChainBlock)
    ->
    MSCpar{vBit::Int64, rBit::Int64, nBitT::Int64, nBitA::Int64, depth::Int64}
Get the set of nBit parameters from a circuit built by function `MScircuit`.
Fields:
\n`vBit::Int64`: Number of virtual qubits.
\n`rBit::Int64`: Number of reusable qubits.
\n`nBitT::Int64`: Number of qubits in target register.
\n`nBitA::Int64`: Number of qubits (lines) in MSCircuit.
"""
struct MSCpar
    vBit::Int64  # Number of virtual qubits.
    rBit::Int64  # Number of reusable qubits.
    nBitT::Int64 # Number of qubits in target register.
    nBitA::Int64 # Number of qubits (lines) in MSCircuit.
    depth::Int64 # Depth of each MPS block in MSCircuit.

    function MSCpar(MSCircuit::ChainBlock)
        nBitA = nqubits(MSCircuit)
        vBit = length(MSCircuit[end-1])
        rBit = length(MSCircuit[3][2])
        nBitT = nBitA - 1 - vBit - rBit
        if length(collect_blocks(QMPS.QDiff, MSCircuit)) == 0
            depth = 0
        else
            depth = length(content(MSCircuit[3][1])[1])
        end
        new(vBit, rBit, nBitT, nBitA, depth)
    end
end

"""
    MScircuit(nBitT::Int64, vBit::Int64, rBit::Int64, MPSblocks::Array{CompositeBlock,1}; ϕ::Float64=0.0)
    ->
    MSCircuit::ChainBlock
Function to build quantum circuit for MPS-Swap Test.
"""
function MScircuit(nBitT::Int64, vBit::Int64, rBit::Int64, MPSblocks::Array{CompositeBlock,1}; ϕ::Float64=0.0)
    par2nd = MPSpar(nBitT, vBit, rBit)
    nBlock = par2nd.nBlock
    nBitG = par2nd.nBit
    nBitA = 1 + nBitG + nBitT
    Cblocks = AbstractBlock[]
    push!(Cblocks, put(nBitA, nBitA=>H))
    push!(Cblocks, put(nBitA, nBitA=>shift(ϕ)))
    for i = 1:nBlock
        # MPSblock = put(nBitA, Tuple( (nBitT+1):(nBitA-1), )=>MPSblocks[i] ) #Use `subroutine` instead of `put` for CuYao compatibility and efficiency.
        MPSblock = subroutine( nBitA, MPSblocks[i], (nBitT+1):(nBitA-1) )
        SWAPblock = chain(nBitA, [control(nBitA, nBitA, ( (nBitA-1-irBit), (nBitT-(i-1)*rBit-irBit) )=>SWAP) for irBit = 0:rBit-1])
        OpBlock = chain(nBitA, MPSblock, SWAPblock)
        push!(Cblocks, OpBlock)
        MeasureBlock = Measure(nBitA, locs=(nBitT+vBit+1):(nBitA-1), resetto=0)
        push!(Cblocks, MeasureBlock)
    end
    SWAPvBit = chain(nBitA, [control(nBitA, nBitA, ( (nBitT+i),i )=>SWAP) for i=vBit:-1:1])
    push!(Cblocks, SWAPvBit)
    push!(Cblocks, put(nBitA, nBitA=>H))
    MSCircuit = chain(nBitA,Cblocks)
    MSCircuit
end

"""
    MStestRes{reg::ArrayReg, witnessOp::PutBlock, overlap::Float64}
Fields:
\n`reg::ArrayReg`: The register goes through `MStest`. nqubits(MStest.reg) = 1 + vBit + rBit + nBitT.
\n`witnessOp::PutBlock`: The witness(measure) operator of the `MStest` circuit.
\n`overlap::Float64`: Overlap bewteen register generated by `MScircuit` and target register.
"""
struct MStestRes
    reg::ArrayReg # The register goes through `MStest`. nqubits(MStest.reg) = 1 + vBit + rBit + nBitT.
    witnessOp::PutBlock  # The witness(measure) operator of the `MStest` circuit.
    overlap::Float64     # Overlap bewteen register generated by MSCircuit and target register.
end

"""
    MStest(regT::ArrayReg, MSCircuit::ChainBlock; nMeasure::Int64=1, useCuYao::Bool=CUDA_ON)
    ->
    MStestRes{reg::ArrayReg, witnessOp::PutBlock, overlap::Float64}
Return the main structure of MPS-Swap Test algorithm.
\n `regT::ArrayReg`: Target quantum state.
"""
function MStest(regT::ArrayReg, MSCircuit::ChainBlock; nMeasure::Int64=1, useCuYao::Bool=CUDA_ON)
    nMcheck(nMeasure, regT)
    nBitA = MSCpar(MSCircuit).nBitA
    nBitT = nqubits(regT)
    regA = join(zero_state(nBitA-nBitT, nbatch=nMeasure), repeat(regT,nMeasure) ) #This operation changes regA in the parent scope.
    useCuYao && (regA = regA |> CuYao.cu)
    resS = MStest(MSCircuit, regAll=regA)
    witnessOp = resS.witnessOp
    ActualOverlap = resS.overlap
    regA = resS.reg
    useCuYao && (regA = regA |> CuYao.cpu)
    MStestRes(regA, witnessOp, ActualOverlap)
end
"""
    MStest(MSCircuit::ChainBlock; regAll::ArrayReg)
    ->
    MStestRes{reg::ArrayReg, witnessOp::PutBlock, overlap::Float64}
Method 2 of `MStest`:
\n `regAll::ArrayReg`: Input register(state) for MSCircuit(MSQR circuit) built by `MScircuit`. `regAll = join(zero_state(1+nBitV+nBitR), regT)`
"""
function MStest(MSCircuit::ChainBlock; regAll::ArrayReg)
    regA = copy(regAll)
    nBitA = nqubits(regA)
    witnessOp = put(nBitA, nBitA=>Z)
    regA |> MSCircuit
    Overlaps = expect(witnessOp, regA)
    ActualOverlap = mean(Overlaps) |> real #Take the real part for simplicity since the imaginary part is 0.
    MStestRes(regA, witnessOp, ActualOverlap)
end

"""
    MSTtestRes{Aoverlap::Float64, Eoverlap::Float64, error::Float64}
Fields:
\n`Aoverlap::Float64`: Actual overlap generated by MPS-Swap Test.
\n`Eoverlap::Float64`: Expected overlap calculated through formula.
\n`error::Float64`:    Aoverlap/Eoverlap - 1
"""
struct MSTtestRes
    Aoverlap::Float64 # Actual overlap generated by MPS-Swap Test.
    Eoverlap::Float64 # Expected overlap calculated through formula.
    error::Float64    # Aoverlap/Eoverlap - 1
end

"""
    MSTtest(regT::ArrayReg, MSCircuit::ChainBlock, cExtend::ChainBlock; nMeasure::Int64=1, useCuYao::Bool=CUDA_ON)
    ->
    MSTtestRes{Aoverlap::Float64, Eoverlap::Float64, error::Float64}
Test function that verify the validity of MPS-Swap Test algorithm.
"""
function MSTtest(regT::ArrayReg, MSCircuit::ChainBlock, cExtend::ChainBlock; nMeasure::Int64=1, useCuYao::Bool=CUDA_ON)
    MSTres = MStest(regT, MSCircuit, nMeasure=nMeasure, useCuYao=useCuYao)
    ActualOverlap = MSTres.overlap
    regG = zero_state(nqubits(regT))
    regG |> cExtend
    ExpectOverlap = ((regT.state'*regG.state)[1] |> abs)^2 #[1]: Convert array to complex number.
    error = ActualOverlap/ExpectOverlap - 1
    MSTtestRes(ActualOverlap, ExpectOverlap, error)
end

"""
    nMcheck(nMeasure::Int64, regs::ArrayReg...)
Function that check if the user not only define nMeasure but also define nbatch.
"""
function nMcheck(nMeasure::Int64, regs::ArrayReg...)
    nbMean = nbatch.(regs) |> mean
    if nbMean!=1 && nMeasure!=1
        print("WARNING: Measure Times has already been defined in nbatch(reg)!!!
              \nNow the actual Measure Times = nMeasure*nbatch(reg).\n")
    end
end

# Default setting of CuYao(CUDA) support is off.
CUDA_ON = false


text:
require("$(dirname(@__FILE__))/MyModel.jl")
using MyModel
using Base.Test

init()
@test typeof(model) === MyModel.ModelType
@test MyModel.start_year === 2014
@test curr_year === MyModel.start_year
@test length(model.results_over_time) === 1
@test Epicenter.fetch_records() != {}

advance()
@test MyModel.curr_year === 2015
@test length(model.results_over_time) === 2

advance()
@test MyModel.curr_year === 2016
@test length(model.results_over_time) === 3

reset()
@test typeof(model) === MyModel.ModelType
@test MyModel.start_year === 2014
@test MyModel.curr_year === MyModel.start_year
@test length(model.results_over_time) === 1
@test Epicenter.take_records() != {}

advance()
@test MyModel.curr_year === 2015

# Test setparam
@test MyModel.model.price == MyModel.start_price
Epicenter.setparam([:(model.price), 15.0])
@test MyModel.model.price == 15.0

Epicenter.setparam([:curr_year, 20])
@test MyModel.curr_year == 20

Epicenter.setparam([:(dummy.x), 84])
@test MyModel.dummy.x == 84


text:
# fallback text/plain representation of any type:
writemime(io, ::MIME"text/plain", x) = showlimited(io, x)

function writemime(io::IO, ::MIME"text/plain", f::Function)
    if isgeneric(f)
        n = length(f.env)
        m = n==1 ? "method" : "methods"
        print(io, "$(f.env.name) (generic function with $n $m)")
    else
        show(io, f)
    end
end

function writemime(io::IO, ::MIME"text/plain", v::AbstractVector)
    if isa(v, Range)
        show(io, v)
    else
        print(io, summary(v))
        if !isempty(v)
            println(io, ":")
            with_output_limit(()->print_matrix(io, v))
        end
    end
end

writemime(io::IO, ::MIME"text/plain", v::AbstractArray) =
    with_output_limit(()->showarray(io, v, header=true, repr=false))

function writemime(io::IO, ::MIME"text/plain", v::DataType)
    show(io, v)
    methods(v) # force constructor creation
    if isgeneric(v)
        n = length(v.env)
        m = n==1 ? "method" : "methods"
        print(io, " (constructor with $n $m)")
    end
end

writemime(io::IO, ::MIME"text/plain", t::Associative) =
    showdict(io, t, limit=true)
writemime(io::IO, ::MIME"text/plain", t::Union(KeyIterator,ValueIterator)) =
    showkv(io, t, limit=true)

# showing exception objects as descriptive error messages

showerror(io::IO, e) = show(io, e)

function showerror(io::IO, e::TypeError)
    ctx = isempty(e.context) ? "" : "in $(e.context), "
    if e.expected === Bool
        print(io, "type: non-boolean ($(typeof(e.got))) ",
                  "used in boolean context")
    else
        if isa(e.got,Type)
            tstr = "Type{$(e.got)}"
        else
            tstr = string(typeof(e.got))
        end
        print(io, "type: $(e.func): ",
                  "$(ctx)expected $(e.expected), ",
                  "got $tstr")
        if e.func === :apply && e.expected <: Function && isa(e.got,AbstractArray)
            println(io)
            print(io, "Use square brackets [] for indexing.")
        end
    end
end

function showerror(io::IO, e, bt)
    try
        showerror(io, e)
    finally
        show_backtrace(io, bt)
    end
end

showerror(io::IO, e::LoadError) = showerror(io, e, {})
function showerror(io::IO, e::LoadError, bt)
    showerror(io, e.error, bt)
    print(io, "\nwhile loading $(e.file), in expression starting on line $(e.line)")
end

function showerror(io::IO, e::DomainError, bt)
    print(io, "DomainError")
    for b in bt
        code = ccall(:jl_lookup_code_address, Any, (Ptr{Void},Cint), b, true)
        if length(code) == 4
            if code[1] in (:log, :log2, :log10, :sqrt) # TODO add :besselj, :besseli, :bessely, :besselk
                print(io, "\n", code[1],
                      " will only return a complex result if called with a complex argument.",
                      "\ntry ", code[1], "(complex(x))")
            end
            break
        end
    end
    show_backtrace(io, bt)
end

showerror(io::IO, e::SystemError) = print(io, "$(e.prefix): $(strerror(e.errnum))")
showerror(io::IO, ::DivideError) = print(io, "integer division error")
showerror(io::IO, ::StackOverflowError) = print(io, "stack overflow")
showerror(io::IO, ::UndefRefError) = print(io, "access to undefined reference")
showerror(io::IO, e::UndefVarError) = print(io, e.var, " not defined")
showerror(io::IO, ::EOFError) = print(io, "read: end of file")
showerror(io::IO, e::ErrorException) = print(io, e.msg)
showerror(io::IO, e::KeyError) = (print(io, "key not found: "); show(io, e.key))
showerror(io::IO, e::InterruptException) = print(io, "interrupt")

function showerror(io::IO, e::MethodError)
    name = isgeneric(e.f) ? e.f.env.name : :anonymous
    if isa(e.f, DataType)
        print(io, "no method $(e.f)(")
    else
        print(io, "no method $(name)(")
    end
    for (i, arg) in enumerate(e.args)
        if isa(arg,Type) && arg != typeof(arg)
            print(io, "Type{$(arg)}")
        else
            print(io, typeof(arg))
        end
        i == length(e.args) || print(io,", ")
    end
    print(io, ")")
    if isdefined(Base,name)
        f = eval(Base,name)
        if f !== e.f && isgeneric(f) && applicable(f,e.args...)
            println(io)
            print(io, "you may have intended to import Base.$(name)")
        end
    end
end

function show_trace_entry(io, fname, file, line, n)
    print(io, "\n")
    print(io, " in ", fname, " at ", file)
    if line >= 1
        try
            print(io, ":", line)
        catch
            print(io, '?') #for when dec is not yet defined
        end
    end
    if n > 1
        print(io, " (repeats ", n, " times)")
    end
end

function show_backtrace(io::IO, t, set=1:typemax(Int))
    # we may not declare :eval_user_input
    # directly so that we get a compile error
    # in case its name changes in the future
    show_backtrace(io, try
                         symbol(string(eval_user_input))
                       catch
                         :(:) #for when client.jl is not yet defined
                       end, t, set)
end

# show the backtrace, up to (but not including) top_function
function show_backtrace(io::IO, top_function::Symbol, t, set)
    n = 1
    lastfile = ""; lastline = -11; lastname = symbol("#")
    local fname, file, line
    count = 0
    for i = 1:length(t)
        lkup = ccall(:jl_lookup_code_address, Any, (Ptr{Void},Cint), t[i], true)
        if lkup === ()
            continue
        end
        fname, file, line, fromC = lkup
        if fromC; continue; end
        if i == 1 && fname == :error; continue; end
        if fname == top_function; break; end
        count += 1
        if !in(count, set); continue; end
        if file != lastfile || line != lastline || fname != lastname
            if lastline != -11
                show_trace_entry(io, lastname, lastfile, lastline, n)
            end
            n = 1
            lastfile = file; lastline = line; lastname = fname
        else
            n += 1
        end
    end
    if n > 1 || lastline != -11
        show_trace_entry(io, lastname, lastfile, lastline, n)
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
module iterative

	import Base.BLAS

	include("cg.jl")
	include("gmres.jl")

end


text:
struct CarEvaluation <: Tabular end
url(::CarEvaluation) = "http://archive.ics.uci.edu/ml/machine-learning-databases/car/car.data"
checksum(::CarEvaluation) = "b703a9ac69f11e64ce8c223c0a40de4d2e9d769f7fb20be5f8f2e8a619893d83"
prep(::CarEvaluation) = path -> preprocess(path, CarEvaluation())
target(::CarEvaluation) = 7
categorical(::CarEvaluation) = 1:6
size(::CarEvaluation) = (1728, 0, 0)
problem(::CarEvaluation) = Classification

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
# ---
# title: 262. Trips and Users
# id: problem262
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories:
# link: <https://leetcode.com/problems/trips-and-users/description/>
# hidden: true
# ---
#
# The `Trips` table holds all taxi trips. Each trip has a unique Id, while
# Client_Id and Driver_Id are both foreign keys to the Users_Id at the `Users`
# table. Status is an ENUM type of ('completed', 'cancelled_by_driver',
# 'cancelled_by_client').
#
#
#
#     +----+-----------+-----------+---------+--------------------+----------+
#     | Id | Client_Id | Driver_Id | City_Id |        Status      |Request_at|
#     +----+-----------+-----------+---------+--------------------+----------+
#     | 1  |     1     |    10     |    1    |     completed      |2013-10-01|
#     | 2  |     2     |    11     |    1    | cancelled_by_driver|2013-10-01|
#     | 3  |     3     |    12     |    6    |     completed      |2013-10-01|
#     | 4  |     4     |    13     |    6    | cancelled_by_client|2013-10-01|
#     | 5  |     1     |    10     |    1    |     completed      |2013-10-02|
#     | 6  |     2     |    11     |    6    |     completed      |2013-10-02|
#     | 7  |     3     |    12     |    6    |     completed      |2013-10-02|
#     | 8  |     2     |    12     |    12   |     completed      |2013-10-03|
#     | 9  |     3     |    10     |    12   |     completed      |2013-10-03|
#     | 10 |     4     |    13     |    12   | cancelled_by_driver|2013-10-03|
#     +----+-----------+-----------+---------+--------------------+----------+
#
#
# The `Users` table holds all users. Each user has an unique Users_Id, and Role
# is an ENUM type of ('client', 'driver', 'partner').
#
#
#
#     +----------+--------+--------+
#     | Users_Id | Banned |  Role  |
#     +----------+--------+--------+
#     |    1     |   No   | client |
#     |    2     |   Yes  | client |
#     |    3     |   No   | client |
#     |    4     |   No   | client |
#     |    10    |   No   | driver |
#     |    11    |   No   | driver |
#     |    12    |   No   | driver |
#     |    13    |   No   | driver |
#     +----------+--------+--------+
#
#
# Write a SQL query to find the cancellation rate of requests made by unbanned
# users (both client and driver must be unbanned) between **Oct 1, 2013** and
# **Oct 3, 2013**. The cancellation rate is computed by dividing the number of
# canceled (by client or driver) requests made by unbanned users by the total
# number of requests made by unbanned users.
#
# For the above tables, your SQL query should return the following rows with the
# cancellation rate being rounded to _two_ decimal places.
#
#
#
#     +------------+-------------------+
#     |     Day    | Cancellation Rate |
#     +------------+-------------------+
#     | 2013-10-01 |       0.33        |
#     | 2013-10-02 |       0.00        |
#     | 2013-10-03 |       0.50        |
#     +------------+-------------------+
#
#
# **Credits:**
# Special thanks to
# [@cak1erlizhou](https://leetcode.com/discuss/user/cak1erlizhou) for
# contributing this question, writing the problem description and adding part of
# the test cases.
#
#
## @lc code=start
using LeetCode

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


text:
include("intcode.jl")

# Part 1: run the TEST program with input 1
prog = open(f -> parse.(Int32, split(readline(f), ",")), "day05_input.txt")
println.(run_intcode!(copy(prog), [1]))

# Part 2: run the TEST program with input 5
println.(run_intcode!(copy(prog), [5]))

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
module REPLComboShell

export setup_repl

using PEG

function process_pipeline(x)
    cmds = collect(Iterators.flatten(x))
    pipeline(cmds...)
end

@PEG.rule script = command[1] & pipe_recv[*] |> process_pipeline
@PEG.rule pipe_recv = r"\|"p & command |> x->x[2]
@PEG.rule command = argument[+] |> (x)->Base.cmd_gen(x)
@PEG.rule argument = word, single_quote, double_quote
@PEG.rule single_quote = "'" & r"(\\.|[^'])+" & "'"  |> x->x[2]
@PEG.rule double_quote = "\"" & r"(\\.|[^\"])+" & "\""  |> x->x[2]
@PEG.rule word = r"[^|\s\"']+"p

function is_shell_command(s)
    args = split(s, ' ')
    !isnothing(Sys.which(args[1])) && args[1] != "import" || args[1] == "cd"
end

using REPL
import REPL.LineEdit

struct JulishCompletionProvider <: REPL.LineEdit.CompletionProvider end

def_completions = REPL.REPLCompletionProvider()

function REPL.complete_line(c::JulishCompletionProvider, s::REPL.LineEdit.PromptState)
    partial = REPL.beforecursor(s.input_buffer)
    if is_shell_command(partial)
        completion_command = Base.shell_escape(`complete -C "$partial"`)
        fish_complete = readlines(`fish -c "$completion_command"`)
        fish_complete = map(x->split(x, '\t')[1], fish_complete)
        ret = REPL.REPLCompletions.Completion[REPL.REPLCompletions.PathCompletion(x) for x in fish_complete]
        range = findlast(' ', partial)+1:length(partial)
        ret = unique!(map(REPL.REPLCompletions.completion_text, ret))
        should_complete = length(ret)==1
        return ret, partial[range], should_complete
    else
        return REPL.complete_line(def_completions, s)
    end
end

function setup_repl(repl, default_username)
    julsh = LineEdit.Prompt("test>";
        prompt_prefix = "",
        prompt_suffix = "",
        keymap_dict   = REPL.LineEdit.default_keymap_dict,
        on_enter      = (s -> true),
        complete      = JulishCompletionProvider(),
        sticky        = true,
    )

    function set_prompt(julsh)
        location = if haskey(ENV, "SSH_CONNECTION")
            hostname = gethostname()
            "$hostname "
        else
            ""
        end
        path_full = replace(pwd(), homedir() => "~")
        path_dirname = dirname(path_full)
        path = if length(path_dirname) > 1
            path_parts = split(path_dirname, "/") .|> x->if length(x) > 0 x[1:1] else x end
            push!(path_parts, basename(path_full))
            join(path_parts, "/")
        else
            path_full
        end
        branch = if isfile(".git/HEAD")
            head = readchomp(".git/HEAD")
            matching = match(r"ref: refs/heads/(.+)", head)
            label = if isnothing(matching)
                "detached"
            else
                matching.captures[1]
            end
            " $label"
        else
            ""
        end
        julsh.prompt = repeat(" ", length(location) + length(path) + length(branch) + 2)

        username = ENV["USER"]
        julsh.prompt_suffix = if username == default_username
            "\x1b[G$location\x1b[32m$path\x1b[36m$branch\x1b[39m> "
        elseif username == "root"
            "\x1b[G$location\x1b[31m$path\x1b[33m$branch\x1b[39m# "
        else
            "\x1b[G$location\x1b[33m$path\x1b[32m$branch\x1b[39m% "
        end
    end

    function parse_to_expr(s)
        args = split(s, ' ')
        if !isnothing(Sys.which(args[1])) && args[1] != "import"
            cmd = nothing
            try
                cmd = parse_whole(script, s)
            catch e
                if isa(e, Base.Meta.ParseError)
                    print("\x1b[31m", e.msg, "\x1b[m")
                    return
                else
                    rethrow(e)
                    return
                end
            end
            open(pipeline(stdin, ignorestatus(cmd), stdout), "r") do io
                while true
                    try
                        sleep(0) # Required? Yes. Cursed? Yep!
                        success(io)
                        break
                    catch
                    end
                end
            end
            set_prompt(julsh);
            nothing
        elseif args[1] == "cd"
            if length(args) != 2
                cd()
            elseif !isdir(args[2])
                println("$(args[2]): directory not found")
            else
                cd(expanduser(args[2]));
            end
            set_prompt(julsh);
            nothing
        else
            Meta.parse(s)
        end
    end

    set_prompt(julsh)

    if !isdefined(repl, :interface)
        repl.interface = REPL.setup_interface(repl)
    end

    julsh.on_done = REPL.respond(parse_to_expr, repl, julsh)

    main_mode = repl.interface.modes[1]
    help_mode = repl.interface.modes[3]
    history_mode = repl.interface.modes[4]
    pkg_mode = repl.interface.modes[end]

    history_mode.hp.mode_mapping[:julsh] = julsh
    julsh.hist = history_mode.hp

    search_prompt, skeymap = LineEdit.setup_search_keymap(history_mode.hp)

    julsh.keymap_dict = LineEdit.keymap(Dict{Any,Any}[
        skeymap,
        REPL.mode_keymap(main_mode),
        LineEdit.history_keymap,
        LineEdit.default_keymap,
        LineEdit.escape_defaults,
    ])

    julsh.keymap_dict = LineEdit.keymap_merge(julsh.keymap_dict, Dict{Any, Any}(
        '\b' => function (s, args...)
            LineEdit.edit_backspace(s)
        end,
        ']' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                LineEdit.transition(s, pkg_mode)
            else
                LineEdit.edit_insert(s, ']')
            end
        end,
        '?' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, help_mode) do
                    LineEdit.state(s, help_mode).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, '?')
            end
        end,
        ',' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, main_mode) do
                    LineEdit.state(s, main_mode).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, ',')
            end
        end,
        "\x03" => function (s, args...)
            print(LineEdit.terminal(s), "^C\n\n")
            LineEdit.transition(s, julsh)
            LineEdit.transition(s, :reset)
            LineEdit.refresh_line(s)
        end,
    ))

    main_mode.keymap_dict['\b'] = function (s, args...)
        if isempty(s) || position(LineEdit.buffer(s)) == 0
            buf = copy(LineEdit.buffer(s))
            LineEdit.transition(s, julsh) do
                LineEdit.state(s, julsh).input_buffer = buf
            end
        else
            LineEdit.edit_backspace(s)
        end
    end

    pkg_mode.keymap_dict = LineEdit.keymap_merge(pkg_mode.keymap_dict, Dict{Any, Any}(
        '\b' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, julsh) do
                    LineEdit.state(s, julsh).input_buffer = buf
                end
            else
                LineEdit.edit_backspace(s)
            end
        end,
        "\x03" => function (s, args...)
            print(LineEdit.terminal(s), "^C\n\n")
            LineEdit.transition(s, julsh)
            LineEdit.transition(s, :reset)
            LineEdit.refresh_line(s)
        end,
    ))

    help_mode.keymap_dict = LineEdit.keymap_merge(help_mode.keymap_dict, Dict{Any, Any}(
        '\b' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, julsh) do
                    LineEdit.state(s, julsh).input_buffer = buf
                end
            else
                LineEdit.edit_backspace(s)
            end
        end,
        "\x03" => function (s, args...)
            print(LineEdit.terminal(s), "^C\n\n")
            LineEdit.transition(s, julsh)
            LineEdit.transition(s, :reset)
            LineEdit.refresh_line(s)
        end,
    ))

    help_mode.on_done = function (s, buf, ok::Bool)
        if !ok
            return LineEdit.transition(s, :abort)
        end
        line = String(take!(buf)::Vector{UInt8})
        f = line->REPL.helpmode(REPL.outstream(repl), line)
        ast = Base.invokelatest(f, line)
        response = REPL.eval_with_backend(ast, REPL.backend(repl))
        REPL.print_response(repl, response, true, REPL.hascolor(repl))
        println()
        return LineEdit.transition(s, julsh)
    end

    repl.interface.modes[1] = julsh;
end

end # module

max_stars_repo_name: aaronjanse/REPLJulish

max_issues_repo_name: aaronjanse/REPLJulish

max_forks_repo_name: aaronjanse/REPLJulish


text:
include("Kernals.jl")
include("gaussian_process.jl")

function uniform(a,b,N)
    rand(N)*(b-a)+a
end

using Distributions
n=50 #number of test points
N=10;#Number of training points
Xtest=linspace(-5,5,n); #Xtest are all of available points to check on the axis
X=uniform(-5,5,N); #X is the values we will sample at and train
noise_dist=Normal(0,10.0^(-6)) #Distribution of Noise
Y=sin(X)+rand(noise_dist) #Y is our function values with noise
K=cov_gen(std_exp_square_ker,Xtest,Xtest)+eye(length(Xtest))*1e-6
L=ctranspose(chol(K));
dist=MvNormal(zeros(n),eye(K)) # We draw f from N(0,I)L which is the same as N(0,K)
fprior=L*rand(dist,3); #Each column is a function

using PyPlot
fig = figure("pyplot_plot",figsize=(5,5))
ax = axes()
fill_between(Xtest,-2*diag(K),2*diag(K),facecolor="#a6a6a6")#This fills confindence interval for two standard
#deviations, currently the variance is one (as we take off diagonal)
plot(Xtest,fprior,alpha=0.75)
title("N=3 Multivariate Gaussians")
ylabel("f(x)")
xlabel("x")
grid("off")
show()

D=[(X[i],Y[i]) for i=1:length(X)];
mu,sigma,D = gaussian_process(std_exp_square_ker,D,1e-6,Xtest);
mu=reshape(mu,length(mu));
sigma=reshape(sigma,length(sigma));
y=map(x->x[2],D);   # these are our y noisy function
x=map(x->x[1],D); #These are our x training points
using PyPlot
fig = figure("pyplot_plot",figsize=(5,5))
ax = axes()
fill_between(Xtest,mu-2*sigma,mu+2*sigma,facecolor="#a6a6a6")#This fills confindence interval for two standard
#deviations, currently the variance is one (as we take off diagonal)
plot(x,y,linewidth=0,marker="o")

title("N=3 Multivariate Gaussians")
ylabel("f(x)")
xlabel("x")
grid("off")
show()


text:
"""
    Uniqueness
"""
abstract type Uniqueness end

struct AllUniqueTrait <: Uniqueness end
const AllUnique = AllUniqueTrait()

struct NotUniqueTrait <: Uniqueness end
const NotUnique = NotUniqueTrait()

struct UnkownUniqueTrait <: Uniqueness end
const UnkownUnique = UnkownUniqueTrait()

Uniqueness(::T) where {T} = Uniqueness(T)
Uniqueness(::Type{T}) where {T} = UnkownUnique
Uniqueness(::Type{T}) where {T<:AbstractRange} = AllUnique

"""
    all_unique(x, u::Uniqueness=UnkownUnique)

Functions similarly to `allunique` but provides the option to bypass checking
each element if this has already been checked or if it's already known at
compile time (e.g., ranges). Note that the default `UnkownUniqueTrait` is
treated like a wild card, so that whatever `Uniqueness(x)` returns determines
whether each element needs to be checked.
"""
function all_unique(x::T, u::Uniqueness=UnkownUnique) where {T}
    return _catch_all_unique(x, _all_unique(Uniqueness(T), u))
end
_all_unique(tu::U, u::U) where {U<:Uniqueness} = tu
_all_unique(::AllUniqueTrait, ::UnkownUniqueTrait) = AllUnique
_all_unique(::UnkownUniqueTrait, ::AllUniqueTrait) = AllUnique
_all_unique(::NotUniqueTrait, ::UnkownUniqueTrait) = NotUnique
_all_unique(::UnkownUniqueTrait, ::NotUniqueTrait) = NotUnique
_all_unique(::AllUniqueTrait, ::NotUniqueTrait) = error
_all_unique(::NotUniqueTrait, ::AllUniqueTrait) = error

_catch_all_unique(x, ::NotUniqueTrait) = false
_catch_all_unique(x, ::AllUniqueTrait) = true
_catch_all_unique(x, ::UnkownUniqueTrait) = allunique(x)
function _catch_all_unique(x::T, ::typeof(error)) where {T}
    error("$T cannot have it's uniqueness specified differently than what is
           determined at compile time. Consider not specifying  `u` in
           `all_unique`.")
end


text:
using Documenter, HomotopyContinuation

makedocs(
    sitename = "HomotopyContinuation.jl",
    pages = [
        "Introduction" => "index.md",
        "Problem formulation" => [
            "ModelKit" => "model_kit.md",
            "Systems" => "systems.md",
            "Homotopies" => "homotopies.md",
            "Linear Subspaces" => "linear_subspaces.md",
        ],
        "Solving Systems" => [
            "Solve (finitely many solutions)" => "solve.md",
            "Witness Sets (infinitely many solutions)" => "witness_sets.md",
            "Monodromy" => "monodromy.md",
            "Results" => "result.md",
            "Certification" => "certification.md",
            "Examples" => "solve_examples.md",
            "Start systems" => "start_systems.md",
        ],
        "Trackers" => ["Tracker" => "tracker.md", "EndgameTracker" => "endgame_tracker.md"],
        "Miscellaneous" => "misc.md",
    ],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    strict = false,
)

deploydocs(
    repo = "github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl.git",
    push_preview = false,
)