text:
#################################################################################################
# AIM 1: nanopolish's LRT vs Accuracy
#################################################################################################
## Deps
using Dierckx
using CodecZlib
using StatsPlots
using Distributions
using DelimitedFiles
using Plots.PlotMeasures

## Constants
const noise_levels = ["$(s)" for s in collect(2.0:0.5:3.5)]
const lr_ints = [-1000:0.01:-10,-10:0.01:-7.5,-7.5:0.01:-5,-5:0.01:-4,-4:0.01:-3,-3:0.01:-2,
                 -2:0.01:-1,-1:0.01:-0.5,-0.5:0.01:0.5,0.5:0.01:1,1:0.01:2,2:0.01:3,3:0.01:4,
                 4:0.01:5,5:0.01:7.5,7.5:0.01:10,10:0.01:1000]
const lr_labs = ["$(minimum(x)):$(maximum(x))" for x in lr_ints]
const data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/Aim-1"

# Colors
const blnd_col = ["#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7"]
const s_color_cde = Dict("0.5" => blnd_col[1], "1.0" => blnd_col[2], "1.5" => blnd_col[3], "2.0" => blnd_col[4],
                         "2.5" => blnd_col[5], "3.0" => blnd_col[6], "3.5" => blnd_col[7], "4.0" => blnd_col[8])

## Default attributes
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

#################################################################################################
# Functions
#################################################################################################

# Function to compute accuracy
function comp_accu(truth, pred)

    # Return accuracy
    return sum(truth .== pred) / length(truth)

end

function map_noise_ex(s, lr_ints)

    # Read in
    stream = GzipDecompressorStream(open("$(data_dir)/nanopolish/gm12878_chr22_sigma_$(s)_nanopolish_tuples_sample.tsv.gz"))
    in_data = readdlm(stream)
    close(stream)

    # Remove NaN
    kp_ind = .! map(isequal("n/a"), in_data[:,2])
    in_data = in_data[kp_ind,:]

    # Get data
    accs = []
    n_pts = []
    total_pts = 0
    for lr in lr_ints
        # Get LRs
        obs_lrs = in_data[:,2]
        lr_ind =  [x in lr for x in obs_lrs]
        if sum(lr_ind) > 0
            # Get calls
            true_x = in_data[lr_ind,1]
            pred_x = in_data[lr_ind,3]
            # Update total calls
            total_pts += sum(lr_ind)
            # Push
            push!(n_pts, sum(lr_ind))
            push!(accs, comp_accu(true_x, pred_x))
        else
            # Push
            push!(n_pts, 0)
            push!(accs, 1.0)
        end
        println("$(lr)")
    end

    # Return accuracies
    return accs, n_pts / total_pts

end

#################################################################################################
# Calls
#################################################################################################

# Plot
n_bins_int = div(length(lr_ints), 2)
plt_crv = plot(xlabel="percentage of calls", ylabel="accuracy", xlim=(0, 1), ylim=(0.5, 1.0), legend=:bottomleft);
for s in noise_levels
    println("Sigma $(s)")

    # Get data
    accs, n_pts = map_noise_ex(s, lr_ints)

    # Plot bars
    p1 = plot(lr_labs,accs,seriestype=:bar,xrotation=45,ylabel="accuracy",title="\\sigma=$(s)",
        ylim=(0.45, 1),label="",top_margin=10px);
    p2 = plot(lr_labs,n_pts,seriestype=:bar,xrotation=45,xlabel="LR",ylabel="percentage of calls",
        ylim=(0.0, 0.51),label="");
    p = plot(p1, p2, layout=(2, 1), size=(700, 700))
    savefig("$(data_dir)/Nanopolish-LR-vs-Accuracy-Sigma-$(s).pdf")

    # Plot curve percent vs accuracy
    acc_curve = 0.5 * ( accs[1:(1 + n_bins_int)] + accs[end:-1:(1 + n_bins_int)])
    perc_curve = n_pts[1:n_bins_int] + n_pts[end:-1:(end - n_bins_int + 1)]
    push!(perc_curve, n_pts[n_bins_int + 1])
    perc_curve = [sum(perc_curve[1:i]) for i = 1:length(perc_curve)]
    acc_curve = [sum(acc_curve[1:i] .* perc_curve[1:i]) / sum(perc_curve[1:i]) for i = 1:length(perc_curve)]

    # Add (0,1) point
    pushfirst!(perc_curve, 0.0)
    pushfirst!(acc_curve, 1.0)

    # Fit spline curve
    spl = Spline1D(perc_curve, acc_curve, k=3)

    # Plot spline
    plot!(plt_crv, 0.0:0.01:1.0, spl(0.0:0.01:1.0), label="\\sigma=$(s)")

end
plot(plt_crv, size=(700, 700))
savefig("$(data_dir)/Nanopolish-Accuracy-vs-Percentage.pdf")


text:
"""
    LinearKernel(; c::Real=0.0)

Linear kernel with constant offset `c`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the linear kernel with constant offset
``c \\geq 0`` is defined as
```math
k(x, x'; c) = x^\\top x' + c.
```

See also: [`PolynomialKernel`](@ref)
"""
struct LinearKernel{Tc<:Real} <: SimpleKernel
    c::Vector{Tc}

    function LinearKernel(c::Real)
        @check_args(LinearKernel, c, c >= zero(c), "c ≥ 0")
        return new{typeof(c)}([c])
    end
end

LinearKernel(; c::Real=0.0) = LinearKernel(c)

@functor LinearKernel

kappa(κ::LinearKernel, xᵀy::Real) = xᵀy + only(κ.c)

metric(::LinearKernel) = DotProduct()

Base.show(io::IO, κ::LinearKernel) = print(io, "Linear Kernel (c = ", only(κ.c), ")")

"""
    PolynomialKernel(; degree::Int=2, c::Real=0.0)

Polynomial kernel of degree `degree` with constant offset `c`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the polynomial kernel of degree
``\\nu \\in \\mathbb{N}`` with constant offset ``c \\geq 0`` is defined as
```math
k(x, x'; c, \\nu) = (x^\\top x' + c)^\\nu.
```

See also: [`LinearKernel`](@ref)
"""
struct PolynomialKernel{Tc<:Real} <: SimpleKernel
    degree::Int
    c::Vector{Tc}

    function PolynomialKernel{Tc}(degree::Int, c::Vector{Tc}) where {Tc}
        @check_args(PolynomialKernel, degree, degree >= one(degree), "degree ≥ 1")
        @check_args(PolynomialKernel, c, only(c) >= zero(Tc), "c ≥ 0")
        return new{Tc}(degree, c)
    end
end

function PolynomialKernel(; degree::Int=2, c::Real=0.0)
    return PolynomialKernel{typeof(c)}(degree, [c])
end

# The degree of the polynomial kernel is a fixed discrete parameter
function Functors.functor(::Type{<:PolynomialKernel}, x)
    reconstruct_polynomialkernel(xs) = PolynomialKernel{typeof(xs.c)}(x.degree, xs.c)
    return (c=x.c,), reconstruct_polynomialkernel
end

kappa(κ::PolynomialKernel, xᵀy::Real) = (xᵀy + only(κ.c))^κ.degree

metric(::PolynomialKernel) = DotProduct()

function Base.show(io::IO, κ::PolynomialKernel)
    return print(io, "Polynomial Kernel (c = ", only(κ.c), ", degree = ", κ.degree, ")")
end

ext: jl


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# ––––––––
# Emphasis
# ––––––––

type Italic
    text
end

@trigger '*' ->
function asterisk_italic(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "*")
    return result === nothing ? nothing : Italic(parseinline(result, md))
end

type Bold
    text
end

@trigger '*' ->
function asterisk_bold(stream::IO, md::MD)
    result = parse_inline_wrapper(stream, "**")
    return result === nothing ? nothing : Bold(parseinline(result, md))
end

# ––––
# Code
# ––––

@trigger '`' ->
function inline_code(stream::IO, md::MD)
    withstream(stream) do
        ticks = startswith(stream, r"^(`+)")
        result = readuntil(stream, ticks)
        if result === nothing
            nothing
        else
            result = strip(result)
            # An odd number of backticks wrapping the text will produce a `Code` node, while
            # an even number will result in a `LaTeX` node. This allows for arbitary
            # backtick combinations to be embedded inside the resulting node, i.e.
            #
            # `a`, ``a``, `` `a` ``, ``` ``a`` ```, ``` `a` ```, etc.
            #  ^     ^        ^            ^             ^
            #  C     L        L            C             C       with C=Code and L=LaTeX.
            #
            isodd(length(ticks)) ? Code(result) : LaTeX(result)
        end
    end
end

# ––––––––––––––
# Images & Links
# ––––––––––––––

type Image
    url::String
    alt::String
end

@trigger '!' ->
function image(stream::IO, md::MD)
    withstream(stream) do
        startswith(stream, "![") || return
        alt = readuntil(stream, ']', match = '[')
        alt ≡ nothing && return
        skipwhitespace(stream)
        startswith(stream, '(') || return
        url = readuntil(stream, ')', match = '(')
        url ≡ nothing && return
        return Image(url, alt)
    end
end

type Link
    text
    url::String
end

@trigger '[' ->
function link(stream::IO, md::MD)
    withstream(stream) do
        startswith(stream, '[') || return
        text = readuntil(stream, ']', match = '[')
        text ≡ nothing && return
        skipwhitespace(stream)
        startswith(stream, '(') || return
        url = readuntil(stream, ')', match = '(')
        url ≡ nothing && return
        return Link(parseinline(text, md), url)
    end
end

@trigger '[' ->
function footnote_link(stream::IO, md::MD)
    withstream(stream) do
        regex = r"^\[\^(\w+)\]"
        str = startswith(stream, regex)
        if isempty(str)
            return
        else
            ref = match(regex, str).captures[1]
            return Footnote(ref, nothing)
        end
    end
end

# –––––––––––
# Punctuation
# –––––––––––

type LineBreak end

@trigger '\\' ->
function linebreak(stream::IO, md::MD)
    if startswith(stream, "\\\n")
        return LineBreak()
    end
end

@trigger '-' ->
function en_dash(stream::IO, md::MD)
    if startswith(stream, "--")
        return "–"
    end
end

const escape_chars = "\\`*_#+-.!{}[]()\$"

@trigger '\\' ->
function escapes(stream::IO, md::MD)
    withstream(stream) do
        if startswith(stream, "\\") && !eof(stream) && (c = read(stream, Char)) in escape_chars
            return string(c)
        end
    end
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using .MyExample
using Test

#@test my_f(2,1) == 7

@testset "MyExample.jl" begin
    @test my_f(2,1) == 7
    @test my_f(2,1) == 13
end


text:
using ZMQ

module zmqtest

# Write your package code here.

end

function ZMQ.recv(socket::ZMQ.Socket, timeout::Number)
    @info "try recv"
    t = @task (_ -> ZMQ.recv(socket))
    schedule(t)
    timedwait(_ -> istaskdone(t), timeout)
    return istaskdone(t) ? fetch(t) : nothing
end

REQUEST_TIMEOUT = 2500

max_stars_repo_name: jwdevantier/julia-zmqtest

max_issues_repo_name: jwdevantier/julia-zmqtest

max_forks_repo_name: jwdevantier/julia-zmqtest


text:
"""
Implementation of a fixed-lag smoother based on Wikipedia
(https://en.wikipedia.org/wiki/Kalman_filter#Fixed-lag_smoother) and
Roger Labbe Jr.'s filterpy package (https://filterpy.readthedocs.io/en/latest/kalman/FixedLagSmoother.html)

Author: Wouter Kouw
Last update: 06-01-2022
"""

using Distributions
using LinearAlgebra
using Random

function fixedlag_rts_smoother(y, A, C, Q, R, x0; N=1)
    """
    Fixed-lag Rauch-Tung-Striebel smoother (Th. 8.2)

    This filter is built for a linear Gaussian dynamical system with known
    transition coefficients, process and measurement noise.
    """

    # Dimensionality
    Dx = size(Q,1)
    Dy = size(R,1)

    # Define identity matrix
    Ix = Matrix{Float64}(I,Dx,Dx)

    # Recast process noise to matrix
    if Dx == 1
        if typeof(Q) != Array{Float64,2}
            Q = reshape([Q], 1, 1)
        end
    end
    if Dy == 1
        if typeof(R) != Array{Float64,2}
            R = reshape([R], 1, 1)
        end
    end

    # Time horizon
    T = length(y)

    # Initialize estimate arrays
    mk = zeros(Dx, T)
    Pk = zeros(Dx, Dx, T)

    # Start previous state variable
    m_kmin1 = x0[1]
    P_kmin1 = x0[2]

    for k in 1:T

        # Predict step
        m_pred = A*m_kmin1
        P_pred = A*P_kmin1*A' + Q

        # Update step
        v = y[k] .- C*m_pred
        S = C*P_pred*C' + R
        SI = inv(S)
        K = P_pred*C'*SI

        # Populate arrays
        m_updt = m_pred + K*v
        P_updt = (Ix - K*C)*P_pred
        mk[:,k] = m_pred
        Pk[:,:,k] = P_updt

        # Compute invariants
        CTSI = C'*SI
        A_KC = (A - K*C)'

        # Start hindsight loop
        if k >= N

            # Initialize smoothed covariance
            PS = P_updt
            for i =1:N

                # Smoothed gain
                K = PS*CTSI

                # Smoothed covariance
                PS = PS*A_KC

                # Hindsight index
                ix = k-i+1

                # Populate smoothed arrays
                mk[:,ix] = mk[:,ix] + K*v
                Pk[:,:,ix] = PS
            end
        else
            mk[:,k] = m_updt
        end

        # Update previous states
        m_kmin1 = mk[:,k]
        P_kmin1 = Pk[:,:,k]
    end
    return mk, Pk
end

max_stars_repo_name: wmkouw/Bayesian-smoothers

max_issues_repo_name: wmkouw/Bayesian-smoothers

max_forks_repo_name: wmkouw/Bayesian-smoothers


text:
module V

const NAMES = [
    "E"
    "H"
    "E1"
    "E2"
    "E3"
    "E4"
    "E_E1"
    "H_E3"
    "H_E4"
    "E11"
    "E12"
    "E23"
    "E34"
    "E24"
    "E44"
    "E11P"
    "E12P"
    "E23P"
    "E34P"
    "E24P"
    "E44P"
    "G"
    "S"
    "I"
    "R"
    "O"
    "A"
    "E11G"
    "E11S"
    "E11R"
    "E12G"
    "E12S"
    "E12R"
    "E23G"
    "E23S"
    "E23I"
    "E23R"
    "E34G"
    "E34S"
    "E34I"
    "E34R"
    "E24G"
    "E24S"
    "E24I"
    "E24R"
    "E44G"
    "E44S"
    "E44I"
    "E44R"
    "sigmaG"
    "sigmaS"
    "sigmaI"
    "sigmaR"
    "sigmaA"
    "sigmaSP"
    "sigmaAP"
    "sigmaG_O"
    "sigmaG_A"
    "sigmaSP_G"
    "sigmaAP_S"
    "sigmaAP_I"
    "sigmaAP_R"
    "P3_A"
    "P2"
    "P3"
    "Akt"
    "Aktstar"
    "RsD"
    "RsT"
    "sigmaRP"
    "Raf"
    "Rafstar"
    "MEK"
    "MEKstar"
    "ERK"
    "ERKstar"
    "OP"
    "AP"
    "A_sigmaG_O"
    "sigmaA_G"
    "sigmaA_G_O"
    "sigmaO"
    "E13"
    "E14"
    "E13P"
    "E14P"
    "E13G"
    "E13S"
    "E13I"
    "E13R"
    "E14G"
    "E14S"
    "E14I"
    "E14R"
    "fint"
    "T"
    "E11T"
    "E12T"
    "E23T"
    "E34T"
    "E24T"
    "E44T"
    "E13T"
    "E14T"
    "sigmaAP_T"
    "sigmaT"
    "E1_PT"
    "E2_PT"
    "E4_PT"
    "E_E1_PT"
    "H_E4_PT"
    "pERK"
    "ERK_MEKstar"
    "pERK_MEKstar"
    "ERKPpase"
    "ERKstar_ERKPpase"
    "pERK_ERKPpase"
]

for (index,value) in enumerate(NAMES)
    eval(Meta.parse("const $value=$index"))
end

const NUM = length(NAMES)

end  # module


text:
#
# TODO
#   - caching doesn't work. fix it with mutable structs
#   - fuse composed operators if possible
#   - avoid caching for AffineOp() created with +,-,λ* where
#     the second op is just Identity, or Null
#
###
# AbstractOperator interface
###

SciMLBase.has_adjoint(::AbstractOperator) = true
SciMLBase.has_mul(::AbstractOperator) = true
SciMLBase.has_mul!(::AbstractOperator) = true

function (A::AbstractOperator{<:Number,D})(u::AbstractField{<:Number,D}, p, t::Number) where{D}
    A * u
end

function (A::AbstractOperator{<:Number,D})(du::AbstractField{<:Number,D}, u, p, t::Number) where{D}
    mul!(du, A, u)
end

function LinearAlgebra.mul!(v::AbstractField{<:Number,D},A::AbstractOperator{<:Number,D},u::AbstractField{<:Number,D}) where{D}
    ArgumentError("LinearAlgebra.mul! not defined for $A")
end

function Base.:\(A::AbstractOperator{<:Number,D},u::AbstractField{<:Number,D}) where{D}
    ArgumentError("Operator inversion not defined for $A")
end

function Base.:*(A::AbstractOperator{<:Number,D}, B::AbstractOperator{<:Number,D}) where{D}
#   @warn "Operator fusion not defined for $A * $B. falling back to lazy composition, ∘"
    A ∘ B
end

# caching
function init_cache(A::AbstractOperator,u)
    @error "Caching behaviour not defined for $A"
end

function set_cache(A::AbstractOperator, cache)
    @set! A.cache = cache
    @set! A.isunset = false
    return A
end

Base.size(A::AbstractOperator, d::Integer) = d <= 2 ? size(A)[d] : 1

###
# NullOp
###

""" (Square) Null operator """
struct NullOp{D} <: AbstractOperator{Bool,D} end

Base.zero(::AbstractOperator{<:Number, D}) where{D} = NullOp{D}()
Base.zero(::Type{<:AbstractOperator{<:Number, D}}) where{D} = NullOp{D}()

Base.adjoint(Z::NullOp) = Z
issquare(::NullOp) = true

function LinearAlgebra.mul!(v::AbstractField{<:Number,D}, ::NullOp{D}, u::AbstractField{<:Number,D}) where{D}
    mul!(v,I, false)
end

# overloads
for op in (
           :*, :∘, :+,
          )
    @eval $op(::NullOp{D}, ::AbstractOperator{<:Number,D}) where{D} = NullOp{D}()
    @eval $op(::AbstractOperator{<:Number,D}, ::NullOp{D}) where{D} = NullOp{D}()
end

Base.:-(::NullOp{D}, A::AbstractOperator{<:Number,D}) where{D} = -A
Base.:-(A::AbstractOperator{<:Number,D}, ::NullOp{D}) where{D} = A

###
# IdentityOp
###

""" (Square) Identity operator """
struct IdentityOp{D} <: AbstractOperator{Bool,D} end

Base.one(::AbstractOperator{<:Number, D}) where{D} = IdentityOp{D}()
Base.one(::Type{<:AbstractOperator{<:Number, D}}) where{D} = IdentityOp{D}()

Base.adjoint(Id::IdentityOp) = Id
Base.inv(Id::IdentityOp) = Id

issquare(::IdentityOp) = true

SciMLBase.has_ldiv(::IdentityOp) = true
SciMLBase.has_ldiv!(::IdentityOp) = true

function LinearAlgebra.mul!(v::AbstractField{<:Number,D}, ::IdentityOp{D}, u::AbstractField{<:Number,D}) where{D}
    copy!(v, u)
end

function LinearAlgebra.ldiv!(v::AbstractField{<:Number,D}, ::IdentityOp{D}, u::AbstractField{<:Number,D}) where{D}
    copy!(v, u)
end

function LinearAlgebra.ldiv!(Id::IdentityOp{D}, u::AbstractField{<:Number,D}) where{D}
    u
end

# overloads
for op in (
           :*, :∘,
          )
    @eval $op(::IdentityOp{D}, A::AbstractOperator{<:Number,D}) where{D} = A
    @eval $op(A::AbstractOperator{<:Number,D}, ::IdentityOp{D}) where{D} = A
end

Base.:/(A::AbstractOperator{<:Number,D}, ::IdentityOp{D}) where{D} = A

for op in (
           :*, :∘,
          )
    @eval $op(::NullOp{D}, ::IdentityOp{D}) where{D} = NullOp{D}()
    @eval $op(::IdentityOp{D}, ::NullOp{D}) where{D} = NullOp{D}()
end

###
# AffineOp
###

""" Lazy affine operator combinations αA + βB """
struct AffineOp{T,D,
                Ta <: AbstractOperator{<:Number,D},
                Tb <: AbstractOperator{<:Number,D},
                Tα,Tβ,Tc
               } <: AbstractOperator{T,D}
    A::Ta
    B::Tb
    α::Tα
    β::Tβ

    cache::Tc
    isunset::Bool

    function AffineOp(A::AbstractOperator{Ta,D}, B::AbstractOperator{Tb,D}, α, β,
                      cache = nothing, isunset = cache === nothing) where{Ta,Tb,D}
        T = promote_type(Ta,Tb)
        new{T,D,typeof(A),typeof(B),typeof(α),typeof(β),typeof(cache)}(A, B, α, β, cache, isunset)
    end
end

issquare(A::AffineOp) = issquare(A.A) & issquare(A.B)
function Base.adjoint(A::AffineOp)
    if issquare(A)
        AffineOp(A.A',A.B',A.α, A.β, A.cache, A.isunset)
    else
        AffineOp(A.A',A.B',A.α', A.β')
    end
end

function init_cache(A::AffineOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    cache = A.B * u
end

function Base.:*(A::AffineOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    @unpack A, B, α, β = A
    α * (A * u) + β * (B * u)
end

function LinearAlgebra.mul!(v::AbstractField{<:Number,D}, A::AffineOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    mul!(v, A.A, u)
    lmul!(A.α, v)

    if isunset
        cache = init_cache(A,u)
        A = set_cache(A, cache)
    else
        mul!(A.cache, A.B, u)
    end

    lmul!(A.β, A.cache)
    axpy!(true, A.cache, v)
end

function Base.:+(A::AbstractOperator{<:Number,D}, B::AbstractOperator{<:Number,D},) where{D}
    AffineOp(A,B,true,true)
end

function Base.:-(A::AbstractOperator{<:Number,D}, B::AbstractOperator{<:Number,D}) where{D}
    AffineOp(A,B,true,-true)
end

function Base.:+(A::AbstractOperator{<:Number,D}, λ::Number) where{D}
    Id = IdentityOp{D}()
    AffineOp(A, Id, true, λ)
end

function Base.:-(A::AbstractOperator{<:Number,D}, λ::Number) where{D}
    Id = IdentityOp{D}()
    AffineOp(A, Id, -true, λ)
end

function Base.:*(A::AbstractOperator{<:Number,D}, λ::Number) where{D}
    Z = NullOp{D}()
    AffineOp(A, Z, true, λ)
end

function Base.:/(A::AbstractOperator{<:Number,D}, λ::Number) where{D}
    Z = NullOp{D}()
    AffineOp(A, Z, -true, λ)
end

###
# ComposedOp
###

""" Lazy Composition A ∘ B """
struct ComposedOp{T,D,Ti,To,Tc} <: AbstractOperator{T,D}
    inner::Ti
    outer::To

    cache::Tc
    isunset::Bool

    function ComposedOp(inner::AbstractOperator{Ti,D},
                        outer::AbstractOperator{To,D},
                        cache = nothing,
                        isunset::Bool = cache === nothing
                       ) where{Ti,To,D}
#       @assert size(outer, 1) == size(inner, 2)
        T = promote_type(Ti, To)
        isunset = cache === nothing
        new{T,D,typeof(inner),typeof(outer),typeof(cache)}(inner, outer, cache, isunset)
    end
end

function Base.:∘(outer::AbstractOperator{Number,D}, inner::AbstractOperator{Number,D}) where{D}
    ComposedOp(inner,outer)
end

Base.size(A::ComposedOp) = (size(A.outer, 1), size(A.inner, 2))
Base.adjoint(A::ComposedOp) = A.inner' ∘ A.outer'
Base.inv(A::ComposedOp) = inv(A.inner) ∘ inv(A.outer)

SciMLBase.has_ldiv(A::ComposedOp) = has_ldiv(A.inner) & has_ldiv(A.outer)
SciMLBase.has_ldiv!(A::ComposedOp) = has_ldiv!(A.inner) & has_ldiv!(A.outer)
issquare(A::ComposedOp) = issquare(A.inner) & issquare(A.outer)

function init_cache(A::ComposedOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    cache = A.inner * u
    return cache
end

function Base.:*(A::ComposedOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    @unpack inner, outer = A
    outer * inner * u
end

function Base.:\(A::ComposedOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    @assert has_ldiv(A)
    @unpack inner, outer = A

    outer \ (inner \ u)
end

function LinearAlgebra.mul!(v::AbstractField{<:Number,D}, A::ComposedOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    @unpack inner, outer = A

    if A.isunset
        cache = init_cache(A, u)
        A = set_cache(A, cache) # pass this A to calling context
    else
        mul!(cache, inner, u)
    end

    mul!(v, outer, cache)
end

function LinearAlgebra.ldiv!(A::ComposedOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    @assert has_ldiv!(A)
    @unpack inner, outer = A

    ldiv!(inner, u)
    ldiv!(outer, u)
end

function LinearAlgebra.ldiv!(v::AbstractField{<:Number,D}, A::ComposedOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    @assert has_ldiv!(A)
    @unpack inner, outer = A

    ldiv!(v, inner, u)
    ldiv!(outer, v)
end

###
# InverseOp
###

""" Lazy Inverse Operator """
struct InverseOp{T,D,Ta} <: AbstractOperator{T,D}
    A::Ta

    function InverseOp(A::AbstractOperator{T,D}) where{T,D}
        @assert issquare(A)
        new{T,D,typeof(A)}(A)
    end
end

function Base.:*(A::InverseOp{<:Number,D}, u::AbstractField{<:Number,D}) where{D}
    A.A \ u
end

SciMLBase.has_ldiv(::InverseOp) = true
SciMLBase.has_ldiv!(::InverseOp) = true
issquare(::InverseOp) = true

Base.inv(A::AbstractOperator) = InverseOp(A)
Base.size(A::InverseOp) = size(A.A)
Base.adjoint(A::InverseOp) = inv(A.A')
LinearAlgebra.ldiv!(y, A::InverseOp, x) = mul!(y, A.A, x)
LinearAlgebra.mul!(y, A::InverseOp, x) = ldiv!(y, A.A, x)

# overload inverse for matrices of operators
#function Base.inv(A::Matrix{<:AbstractOperator})
#end
#


text:
module SPTrackingToolkit
using DefaultArrays
using LapSolve
using LinearAlgebra

export track
export SPT

const ENDMARK=-1
const STARTMARK=1
const FAKEMARK=-100
const COMPLETEMARK=0
const MISSINGFRAME=-1

Base.@kwdef struct SPT
    maxtimegap::Int
    maxdist::Float64
    dims::Int
    rest::Int=0
    X::UnitRange{Int}=1:dims
    I::UnitRange{Int}=(dims+1):(dims+rest)
    verbose::Bool=false
end

function segments(config::SPT,Ma,Mb)
    pa=size(Ma,2)
    pb=size(Mb,2)
    C=DefaultArray(Inf,pa,pb)

    @inbounds for i in 1:pa, j in 1:pb
        cost=norm(Ma[config.X,i] .- Mb[config.X,j])
        cost <= config.maxdist || continue
        C[i,j] = cost
    end

    S=solve_stiff_lap(C,1.05)
    solution=zeros(Int,length(S),3)

    @inbounds for i in eachindex(S)
        ia, ib = S[i]
        solution[i,1]=ia
        solution[i,2]=ib
        solution[i,3]=COMPLETEMARK
        ia<0 && (solution[i,3]=STARTMARK)
        ib<0 && (solution[i,3]=ENDMARK)
        (ia<0 && ib<0) && (solution[i,3]=FAKEMARK)
    end

    solution
end

@inline function tC(v::Vector)
    v[1]+length(v)-2
end

function buildtracksegments(config::SPT,seg)
    tracks=Vector{Int}[]
    @inbounds for x in eachrow(seg[1])
        x[end]==FAKEMARK && continue
        push!(tracks,[1;x[1:2]])
    end
    @inbounds for t in 2:length(seg)
        @inbounds for x in eachrow(seg[t])
            x[end]==FAKEMARK && continue
            x[end]==STARTMARK && begin push!(tracks,[t;x[1:2]]); continue; end

            for iT in eachindex(tracks)
                (tC(tracks[iT])!=t) && continue;
                if tracks[iT][end]==x[1]
                    push!(tracks[iT],x[2])
                    break
                end
            end
        end
    end
    tracks
end

function buildcost_closegaps(config::SPT,frames,tracks)
    startseg=findall(x->x[2]==MISSINGFRAME,tracks)
    endseg=findall(x->x[end]==MISSINGFRAME,tracks)

    n=length(startseg)
    m=length(endseg)
    C=DefaultArray(Inf,m,n)
    @inbounds for i in 1:m, j in 1:n  ## GAP filling block
        ti=endseg[i]
        tj=startseg[j]

        if ti==tj
            continue
        end
        frameendidx = findlast(x->x != MISSINGFRAME,tracks[ti][2:end])+1
        framestartidx = findfirst(x->x != MISSINGFRAME,tracks[tj][2:end])+1
        Tend = frameendidx-2+tracks[ti][1]
        Tstart= framestartidx-2+tracks[tj][1]

        Δt=Tstart-Tend-1
        if 0< Δt <= config.maxtimegap
            pend=frames[Tend][config.X,tracks[ti][frameendidx]]
            pstart=frames[Tstart][config.X,tracks[tj][framestartidx]]
            cost=norm(pend.-pstart)/sqrt(Δt+1)
            cost <= config.maxdist || continue
            C[i,j] = cost
        end
    end
    return C
end

function closegaps(config::SPT,frames,origtracks)
    tracks=origtracks
    @label restartclosing
    startseg=findall(x->x[2]==MISSINGFRAME,tracks)
    endseg=findall(x->x[end]==MISSINGFRAME,tracks)
    n=length(startseg)
    m=length(endseg)
    m==0 && n==0 && return origtracks
    C=buildcost_closegaps(config,frames,tracks)
    S=solve_stiff_lap(C)
    newlinks=[(endseg[a],startseg[b]) for (a,b) in S if (0<a<=m && 0<b<=n)]
    config.verbose && println("Started GAP closing: ", length(newlinks)," gaps identified")
    linktr=Vector{Int}[]
    forbidden=Set{Int}()
    if length(newlinks)>0
        @inbounds for l in newlinks
            ti,tj=l
            if ti in forbidden || tj in forbidden
                continue
            end
            frameendidx = findlast(x->x != MISSINGFRAME,tracks[ti][2:end])+1
            framestartidx = findfirst(x->x != MISSINGFRAME,tracks[tj][2:end])+1
            Tend = frameendidx-2+tracks[ti][1]
            Tstart= framestartidx-2+tracks[tj][1]
            Δt=Tstart-Tend-1
            A=tracks[ti][1:frameendidx]
            B=tracks[tj][framestartidx:end]
            newtrack=[A;fill(MISSINGFRAME,Δt);B]
            push!(linktr,newtrack)
            push!(forbidden,ti)
            push!(forbidden,tj)
        end
    else
        config.verbose && println("GAP closing terminated")
        return tracks
    end
    for i in 1:length(tracks)
        if i in forbidden
            continue;
        end
        push!(linktr,tracks[i])
    end
    tracks = linktr
    config.verbose && println("GAP closing is being restarted: closed ",div(length(forbidden),2), " gaps")
    @goto restartclosing
end

function get_tracked(specs::SPT,frames,links)
    allP=Array{Float64,2}[]
    @inbounds for l in links
        t=l[1]
        idx=1
        n_particleinfo=length(specs.X)+length(specs.I)
        P=zeros(n_particleinfo+1,sum(x>0 for x in l)-1)
        for i in 2:length(l)
            f=l[i]
            if f>0
                P[:,idx].=[t;frames[t][:,f]]
                idx+=1
            end
            t+=1

        end
        push!(allP,P)
    end
    sort!(allP,by=x->-size(x,2))
    allP
end

function buildlinks(config::SPT,frames)
    l=length(frames)
    seg=[segments(config,frames[fa],frames[fb]) for (fa,fb) in zip(1:(l-1),2:l)]
    closegaps(config,frames,buildtracksegments(config,seg))
end

function track(config::SPT,frames)
    links=buildlinks(config,frames)
    get_tracked(config,frames,links)
end
end # module


text:
# Recombinations
# ==============
"""
    average(population)

Returns an *one* offspring individual of a multi-parent recombination by averaging `population`.
"""
function average(population::Vector{T}; kwargs...) where {T <: AbstractVector}
    obj = zeros(eltype(T), length(population[1]))
    l = length(population)
    for i in 1:l
        obj += population[i]
    end
    return obj./l
end

"""
    marriage(population)

Returns an *one* offspring individual of a multi-parent recombination by random copying from `population`.
"""
function marriage(population::Vector{T}; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    s = length(population)
    l = length(first(population))
    obj = zeros(eltype(T), l)
    idxs = rand(rng, 1:s, l)
    for i in 1:l
        obj[i] = population[idxs[i]][i]
    end
    return obj
end

# Strategy recombinations
# =======================

"""
    average(ss::Vector{<:AbstractStrategy})

Returns the average value of the mutation parameter ``\\sigma`` of strategies `ss`.
"""
function average(ss::Vector{<:AbstractStrategy}; kwargs...)
    s = copy(first(ss))
    l = length(ss)
    s.σ = mapreduce(s->s.σ/l, +, ss)
    return s
end

# ==================
# Genetic algorithms
# ==================

"""
    genop(v1, v2)

Returns the same parameter individuals `v1` and `v2` as an offspring pair (GENetic No OPeration).
"""
genop(v1::T, v2::T; kwargs...) where {T <: AbstractVector} = (v1,v2)

# Binary crossovers
# -----------------

"""
    SPX(v1, v2)

Single point crossover between `v1` and `v2` individuals.
"""
function SPX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    pos = rand(rng, 1:l)
    for i in pos:l
        vswap!(c1, c2, i)
    end
    return c1, c2
end

"""
    TPX(v1, v2)

Two point crossover between `v1` and `v2` individuals.
"""
function TPX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    from, to = rand(rng, 1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    for i in from:to
        vswap!(c1, c2, i)
    end
    return c1, c2
end

"""
    SHFX(v1, v2)

Shuffle crossover between the parents `v1` and `v2` that performs recombination similar to [`SPX`](@ref) preliminary shuffling these parents.
"""
function SHFX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    l = length(v1)
    # shuffle and perform 1-point XO
    prm = randperm(rng, l)
    tmp1, tmp2 = SPX(view(v1,prm), view(v2,prm), rng=rng)
    # unshuffle offspring
    c1, c2 = similar(tmp1), similar(tmp2)
    for (i,j) in enumerate(prm)
        c1[j] = tmp1[i]
        c2[j] = tmp2[i]
    end
    return c1, c2
end

"""
    UX(v1, v2)

Uniform crossover between `v1` and `v2` individuals.
"""
function UX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    xch = rand(rng, Bool, l)
    for i in 1:l
        if xch[i]
            vswap!(c1, c2, i)
        end
    end
    return c1, c2
end

"""
    BINX(Cr::Real=0.5)

Returns a uniform (binomial) crossover function, see [Recombination Interface](@ref), function with the probability `Cr` [^2].

The crossover probability value must be in unit interval, ``Cr \\in [0,1]``.
"""
function BINX(Cr::Real = 0.5)
    function binxvr(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
        l = length(v1)
        c1 = copy(v1)
        c2 = copy(v2)
        j = rand(rng, 1:l)
        for i in (((1:l).+j.-2).%l).+1
            if rand(rng) <= Cr
                vswap!(c1, c2, i)
            end
        end
        return c1, c2
    end
    return binxvr
end

"""
    EXPX(Cr::Real=0.5)

Returns an exponential crossover function, see [Recombination Interface](@ref), function with the probability `Cr` [^2].

The crossover probability value must be in unit interval, ``Cr \\in [0,1]``.
"""
function EXPX(Cr::Real = 0.5)
    function expxvr(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
        l = length(v1)
        c1 = copy(v1)
        c2 = copy(v2)
        j = rand(rng, 1:l)
        switch = true
        for i in (((1:l).+j.-2).%l).+1
            i == j && continue
            if switch && rand(rng) <= Cr
                c1[i] = v1[i]
            else
                switch = false
                c1[i] = v2[i]
            end
        end
        return c1, copy(c1)
    end
    return expxvr
end

"""
    BSX(k::Int)

Binary Subset Crossover[^7]. Produces an offspring by first pooling the unique
items of the two parents, and then creating each offspring by sampling without
replacement at most `k` elements from the pool of items.
"""
function BSX(k::Int)
    function BSX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector{Bool}}
        l = length(v1) # get number of available elements
        pooled = findall(v1 .| v2) # pool parents selections
        K = min(k,length(pooled)) # cannot sample more than the items in pool
        c1 = falses(l) # init child 1
        c2 = falses(l) # init child 2
        c1[shuffle(rng, pooled)[1:K]] .= true # fill child 1 with sample from pool w/o replacement
        c2[shuffle(rng, pooled)[1:K]] .= true # fill child 2 with sample from pool w/o replacement
        return c1, c2
    end
    return BSX
end

# Real valued crossovers
# ----------------------

"""
    DC(v1, v2)

Returns a randomly assembled offspring and its inverse from the elements of parents `v1` and `v2`.
"""
function DC(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    l = length(v1)
    c1 = similar(v1)
    c2 = similar(v2)
    sltc = rand(rng, Bool, 2, l)
    for i in 1:l
        c1[i] = sltc[1,i] ? v1[i] : v2[i]
        c2[i] = sltc[2,i] ? v2[i] : v1[i]
    end
    return c1, c2
end

"""
    WAX(w::Vector{<:Real})(v1, v2)

Returns a weighted average recombination operation, see [Recombination Interface](@ref), which generate an offspring as weighted average of the parents `v1` and `v2` with the weights `w`.
"""
function WAX(w::Vector{<:Real})
    function wavexvr(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
        c1 = (v1+v2)./w
        return c1, copy(c1)
    end
    return wavexvr
end

"""
    AX(v1, v2)

Average crossover generates an offspring by taking average of the parents `v1` and `v2`.
"""
function AX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    c1 = (v1+v2)./2
    return c1, copy(c1)
end

"""
    IC(d::Real=0.0)

Returns an extended intermediate recombination operation, see [Recombination Interface](@ref), which generates offspring `u` and `v` as

- ``u_i = x_i + \\alpha_i (y_i - x_i)``
- ``v_i = y_i + \\alpha_i (x_i - y_i)``

where ``\\alpha_i`` is chosen uniform y in the interval ``[-d;d+1]``.

"""
function IC(d::Real = 0.0)
    function intermxvr(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
        l = length(v1)
        α = (1.0+2d) * rand(rng, l) .- d
        c1 = v2 .+ α .* (v1 - v2)
        α = (1.0+2d) * rand(rng, l) .- d
        c2 = v1 .+ α .* (v2 - v1)
        return c1, c2
    end
    return intermxvr
end

"""
    LC(d::Real=0.0)

Returns a extended line recombination operation, see [Recombination Interface](@ref), which generates offspring `u` and `v` as

- ``u_i = x_i + \\alpha (y_i - x_i)``
- ``v_i = y_i + \\alpha (x_i - y_i)``

where ``\\alpha`` is chosen uniform randomly in the interval ``[-d;d+1]``.

"""
function LC(d::Real = 0.0)
    function linexvr(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
        α1, α2 = (1.0+2d) * rand(rng, 2) .- d
        c1 = v2 .+ α2 * (v1 - v2)
        c2 = v1 .+ α1 * (v2 - v1)
        return c1, c2
    end
    return linexvr
end

"""
    HX(x, y)

Heuristic crossover (HX) recombination operation[^3] generates offspring `u` and `v` as

- ``u = x + r (x - y)``
- ``v = y + r (y - x)``

where ``r`` is chosen uniform randomly in the interval ``[0;1)``.
"""
function HX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    c1 = v1 .+ rand(rng)*(v1 .- v2)
    c2 = v2 .+ rand(rng)*(v2 .- v1)
    return c1, c2
end

"""
    LX(μ::Real = 0.0, b::Real = 0.2)

Returns a Laplace crossover (LX) recombination operation[^4], see [Recombination Interface](@ref).
"""
function LX(μ::Real = 0.0, b::Real = 0.2) # location μ, scale b > 0
    function lxxvr(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
        u = rand(rng)
        β = u > 0.5 ? μ + b * log(u) : μ - b * log(u)
        S = β * abs.(v1 - v2)
        c1 = v1 + S
        c2 = v2 + S
        return c1, c2
    end
    return lxxvr
end

"""
    MILX(μ::Real = 0.0, b_real::Real = 0.15, b_int::Real = 0.35)

Returns a mixed integer Laplace crossover (MI-LX) recombination operation[^5], see [Recombination Interface](@ref).
"""
function MILX(μ::Real = 0.0, b_real::Real = 0.15, b_int::Real = 0.35) # location μ, scale b > 0
    function milxxvr(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
        @assert all([typeof(a) == typeof(b) for (a, b) in zip(v1, v2)]) "Types of variables in vectors do not match"
        l = length(v1)
        U, R = rand(rng, l), rand(rng, l)
        B = (isa(x, Integer) ? b_int : b_real for x in v1)
        βs = broadcast((u, r, b) -> r > 0.5 ? μ + b * log.(u) : μ - b * log.(u), U, R, B)
        S = βs .* abs.(v1 - v2)
        c1 = v1 + S
        c2 = v2 + S
        return c1, c2
    end
    return milxxvr
end

"""
    SBX(pm::Real = 0.5, η::Integer = 2)

Returns a Simulated Binary Crossover (SBX) recombination operation, see [Recombination Interface](@ref),
with the mutation probability `pm` of the recombinant component, and is the crossover distribution index `η`[^6].
"""
function SBX(pm::Real = 0.5, η::Integer = 2)
    function sbxv(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
        n = length(v1)
        u = rand(rng, n)
        mask = u .<= 0.5
        mask_neg = u .> 0.5
        β = similar(v1)
        β[mask] .= (2*u[mask]).^(1/(η+1))
        β[mask_neg] .= (2*(1 .- u[mask_neg])).^(-1/(η+1))
        μ = (v1 + v2)./2
        diff = v1 - v2
        c = β.*diff./2
        mask_set = n == 1 ? [1] : rand(rng, n) .<= pm
        c1 = copy(μ) # c2 = μ - c
        c1[mask_set] .-= c[mask_set]
        c2 = copy(μ) # c2 = μ + c
        c2[mask_set] .+= c[mask_set]
        return c1, c2
    end
    return sbxv
end

# Combinatorial crossovers
# ------------------------

"""
    PMX(v1, v2)

Partially mapped crossover which maps ordering and values information from
the parents `v1` and `v2` to  the offspring. A portion of one parent is mapped onto
a portion of the other parent string and the remaining information is exchanged.
"""
function PMX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    s = length(v1)
    from, to = rand(rng, 1:s, 2)
    from, to = from > to ? (to, from)  : (from, to)
    c1 = similar(v1)
    c2 = similar(v2)

    # Swap
    c1[from:to] = v2[from:to]
    c2[from:to] = v1[from:to]

    # Fill in from parents
    for i in vcat(1:from-1, to+1:s)
        # Check conflicting offspring
        in1 = inmap(v1[i], c1, from, to)
        if in1 == 0
            c1[i] = v1[i]
        else
            tmpin = in1
            while tmpin > 0
                tmpin = inmap(c2[in1], c1, from, to)
                in1 = tmpin > 0 ? tmpin : in1
            end
            c1[i] = v1[in1]
        end

        in2 = inmap(v2[i], c2, from, to)
        if in2 == 0
            c2[i] = v2[i]
        else
            tmpin = in2
            while tmpin > 0
                tmpin = inmap(c1[in2], c2, from, to)
                in2 = tmpin > 0 ? tmpin : in2
            end
            c2[i] = v2[in2]
        end
    end
    return c1, c2
end

"""
    OX1(v1, v2)

Order crossover constructs an offspring by choosing a substring of one parent
and preserving the relative order of the elements of the other parent.
"""
function OX1(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    s = length(v1)
    from, to = rand(rng, 1:s, 2)
    from, to = from > to ? (to, from)  : (from, to)
    c1 = zero(v1)
    c2 = zero(v2)
    # Swap
    c1[from:to] = v2[from:to]
    c2[from:to] = v1[from:to]
    # Fill in from parents
    k = to+1 > s ? 1 : to+1 #child1 index
    j = to+1 > s ? 1 : to+1 #child2 index
    for i in vcat(to+1:s,1:from-1)
        while in(v1[k],c1)
            k = k+1 > s ? 1 : k+1
        end
        c1[i] = v1[k]
        while in(v2[j],c2)
            j = j+1 > s ? 1 : j+1
        end
        c2[i] = v2[j]
    end
    return c1, c2
end

"""
    CX(v1, v2)

Cycle crossover creates an offspring from the parents `v1` and `v2` such that
every position is occupied by a corresponding element from one of the parents.
"""
function CX(v1::T, v2::T; kwargs...) where {T <: AbstractVector}
    s = length(v1)
    c1 = zero(v1)
    c2 = zero(v2)
    c1i = fill(true, length(v1))
    c2i = fill(true, length(v2))
    Z = zero(eltype(T))

    f1 = true #switch
    k = 1
    while k > 0
        idx = k
        if f1
            #cycle from v1
            while c1i[idx]
                c1[idx] = v1[idx]
                c2[idx] = v2[idx]
                c1i[idx] = false
                c2i[idx] = false
                idx = inmap(v2[idx],v1,1,s)
            end
        else
            #cycle from v2
            while c2i[idx]
                c1[idx] = v2[idx]
                c2[idx] = v1[idx]
                c1i[idx] = false
                c2i[idx] = false
                idx = inmap(v1[idx],v2,1,s)
            end
        end
        f1 ⊻= true
        k = inmap(true,c2i,1,s)
    end
    return c1,c2
end

"""
    OX2(v1, v2)

Order-based crossover selects at random several positions in the parent `v1`, and
the order of the elements in the selected positions of the parent `v1` is imposed on
the parent `v2`.
"""
function OX2(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    s = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    Z = zero(eltype(T))

    for i in 1:s
        if rand(rng, Bool)
            idx1 = inmap(v2[i],v1,1,s)
            idx2 = inmap(v1[i],v2,1,s)
            c1[idx1] = Z
            c2[idx2] = Z
        end
    end

    for i in 1:s
        if !in(v2[i],c1)
            tmpin = inmap(Z,c1,1,s)
            c1[tmpin] = v2[i]
        end
        if !in(v1[i],c2)
            tmpin = inmap(Z,c2,1,s)
            c2[tmpin] = v1[i]
        end
    end
    return c1,c2
end

"""
    POS(v1, v2)

Position-based crossover is a modification of the [`OX1`](@ref) operator.
It selects a random set of positions in the parents `v1` and `v2`, then imposes
the position of the selected elements of one parent on the corresponding elements
of the other parent.
"""
function POS(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    s = length(v1)
    c1 = zero(v1)
    c2 = zero(v2)
    Z = zero(eltype(T))

    for i in 1:s
        if rand(rng, Bool)
            c1[i] = v2[i]
            c2[i] = v1[i]
        end
end

    for i in 1:s
        if !in(v1[i],c1)
            tmpin = inmap(Z,c1,1,s)
            c1[tmpin] = v1[i]
        end
        if !in(v2[i],c2)
            tmpin = inmap(Z,c2,1,s)
            c2[tmpin] = v2[i]
        end
    end
    return c1,c2
end

"""
    SSX(v1, v2)

Subset crossover operator creates new offspring by pooling the unique indices of
the two parent vectors `v1` and `v2`, and then sampling a set of unique indices
from this pool, uniformly at random.
"""
function SSX(v1::T, v2::T; rng::AbstractRNG=default_rng()) where {T <: AbstractVector}
    l = length(v1)
    pool = unique([v1;v2])
    p = length(pool)
    c1, c2 = if l <= p
        pool[randperm(p)[1:l]], pool[randperm(p)[1:l]]
    else
        rand(pool, l), rand(pool, l)
    end
    return c1, c2
end

# ===================
# Genetic Programming
# ===================

"""
    crosstree(t1::Expr, t2::Expr)

Perform an arbitrary subtree swap between the expressions `t1` and `t2`.
"""
function crosstree(t1::Expr, t2::Expr; rng::AbstractRNG=default_rng())
    tt1, tt2 = copy(t1), copy(t2)
    i, j = rand(rng, 1:nodes(t1)-1), rand(rng, 1:nodes(t2)-1)
    ex1 = tt1[i]
    ex2 = tt2[j]
    tt1[i] = ex2
    tt2[j] = ex1
    tt1, tt2
end

# Utilities
# =========
function vswap!(v1::T, v2::T, idx::Int) where {T <: AbstractVector}
    val = v1[idx]
    v1[idx] = v2[idx]
    v2[idx] = val
end

function inmap(v::T, c::AbstractVector{T}, from::Int, to::Int) where {T}
    exists = 0
    for j in from:to
        if exists == 0 && v == c[j]
            exists = j
        end
    end
    return exists
end


text:
using SummationByParts

function nodecalc(sbp::TetSBP, isDG::Bool)
  vtx = [0.0 0 0; 1 0 0; 0 1 0; 0 0 1]
  r1 = vtx[1, :]
  r2 = vtx[2, :]
  r3 = vtx[3, :]
  r4 = vtx[4, :]

  T = zeros(3,3)
  T[:, 1] = r2 - r1
  T[:, 2] = r3 - r1
  T[:, 3] = r4 - r1

  # create operator
  if isDG
    coords = SummationByParts.SymCubatures.calcnodes(sbp.cub, vtx)
  else
    coords = calcnodes(sbp, vtx)
  end

  xi = zeros(coords)

  for i=1:size(coords,2)
    xi[:, i] = T\(coords[:, i] - r1.')
  end

  return xi, coords
end

function printCaseStatement(xi)
  nnodes = size(xi, 2)

  f = open("caseout.txt", "w")
  for i=1:nnodes
    xi1 = xi[1, i]
    xi2 = xi[2, i]
    if size(xi, 1) == 2
      xi3 = 0.0
    else
      xi3 = xi[3, i]
    end
    println(f, "case $(i-1):")
    println(f, "  {xi = Vector3($xi1, $xi2, $xi3); break; }")
  end

  println(f, "default:")
  println(f, "  {xi = Vector3(0, 0, 0); break; }")
  close(f)
end

function printBaryCoords(xi, coords)

  (tmp, nnodes) = size(xi)

  for i=1:nnodes
    println("node ", i, " coords: ", coords[:,i])
  end

  print("\n")

  for i=1:nnodes
    xi4 = 1 - sum(xi[:,i])
    println("node ", i, " coords: ", xi[:,i], " , xi4 = ", xi4)
  end
end

#=
sbp = TetSBP{Float64}(degree=4, reorder=false, internal=false)
xi, coords = nodecalc(sbp, true)
printBaryCoords(xi, coords)
#printCaseStatement(xi)
=#


text:
struct AffDyn{MT<:AbstractMatrix{Float64}, VT<:AbstractVector{Float64}}
    A::MT
    b::VT
end

function Base.:(-)(dyn::AffDyn)
    return AffDyn(-dyn.A, -dyn.b)
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
module PenaltiesLH

using ArgCheck, DocStringExtensions, Parameters

export Penalty, PenaltyVector
# Penalty
export name, compute_penalty, scalar_value, set_scalar_value, use, no_use, is_used
# PenaltyVector
export retrieve, penalty_dict, total_penalty, compute_all_penalties, report_penalties

include("types.jl")
include("penalty.jl")
include("penalty_vector.jl")

function format_number(x; ndigits = 2)
    xFmt = round(x, digits = ndigits);
    return string(xFmt);
end

end # module


text:
module Fermi
import Lints
import TBLIS
import DistributedArrays

include("Core/Options.jl")                             # Top level scope
include("Backend/Tensors.jl")                             # Top level scope
include("Backend/Environment.jl")
include("Backend/ComputeEnvironment.jl")
include("Backend/IO/Output.jl")
include("Backend/contract.jl")
include("Backend/DIIS.jl")
include("Core/PhysicalConstants.jl")
include("Core/AbstractWavefunctions.jl")
include("Core/Geometry.jl")
include("Core/Orbitals.jl")
include("Core/Integrals.jl")
include("Methods/HartreeFock/HartreeFock.jl")
include("Methods/MollerPlesset/MollerPlesset.jl")
include("Methods/ConfigurationInteraction/ConfigurationInteraction.jl")
include("Methods/CoupledCluster/CoupledCluster.jl")
#include("Methods/FocalPointAnalysis/FocalPointAnalysis.jl")

include("Core/SinglePointEnergy.jl")

end # module


text:
using COBS
using Test

# these are examples from [the wiki on COBS](https://en.wikipedia.org/wiki/Consistent_Overhead_Byte_Stuffing#Encoding_examples)
payloads = Vector{Vector{UInt8}}([[0],
            [0, 0],
            [1, 2, 0, 3],
            [1, 2, 3, 4],
            [1, 0, 0, 0],
            collect(1:254)])

messages = Vector{Vector{UInt8}}([[1, 1, 0],
            [1, 1, 1, 0],
            [3, 1, 2, 2, 3, 0],
            [5, 1, 2, 3, 4, 0],
            [2, 1, 1, 1, 1, 0],
            [255; 1:254; 0]])

@testset "encoding" begin
    for (payload, message) in zip(payloads, messages)
        @test encode(payload) == message
    end
end

@testset "encoding to port" begin
    for (payload, message) in zip(payloads, messages)
        io = IOBuffer()
        encode(io, payload)
        @test take!(io) == message
    end
end

@testset "decoding" begin
    for (message, payload) in zip(messages, payloads)
        @test decode(message) == payload
    end
end

@testset "decoding from port" begin
    for (message, payload) in zip(messages, payloads)
        io = IOBuffer(message)
        @test decode(io) == payload
    end
end


text:
include("Multi_HZ.jl")
include("Granular_Module.jl")
include("calculate_Pressure.jl")
include("Min_Energy.jl")
include("Deform_Config.jl")
using JuLIP
using Minimise_Energy
using Deform
using Pressure

# ===============================
#            new test
# ===============================

Grain_Num = 1024
α = 2.0
ϕ = 0.842

Jamming_Info = readdlm("jamming_config.txt")

x = Jamming_Info[:,1]'
y = Jamming_Info[:,2]'
z = zeros(Grain_Num)'
X = [x; y; z]
X = vecs(X)

grain_type = Jamming_Info[:,4] + 1
r = Jamming_Info[:,3]

Config = gran.setup_cell(Grain_Num, ϕ, α)
Config.X .= X
Config.M = r
Config.Z = grain_type

ϕ = π*sum(r.^2)

set_data!(Config,"Phi",ϕ)

#deformType = "grow"
#deformValue = 1.0
#Config= Deform.deform_Config(Config, deformType, deformValue)

# =======================
#   several grow steps
# =======================

N_Grow = 10
ϕ_Grow = get_data(Config, "Phi") + 10.^linspace(-6, -4, N_Grow)

p_Grow = zeros(N_Grow)'
Config_Grow = Array{AbstractAtoms, 1}(N_Grow)

deformType = "phi"
Config_Temp = deepcopy(Config)
for i = 1:N_Grow
    deformValue = ϕ_Grow[i]
    Config_Grow[i] = Deform.deform_Config(Config_Temp, deformType, deformValue)
    Config_Temp = deepcopy(Config_Grow[i])

    p_Grow[i] = trace(get_data(Config_Temp, "Pressure"))/get_data(Config_Temp, "Dimension")
end


text:
class C {
    public static void main(String[] args) {
	int{H} x = 7;
	int{L} y = 0;
	if (x > 0) {
	    y = 4;
	}
	int{L} z = y;
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
"""
    function eyes()

Evil-looking eyes

# Example

```jldoctest
julia> cowsay("Mwah, ha, ha!", cow=Cowsay.eyes)
 _______________
< Mwah, ha, ha! >
 ---------------
    \\
     \\
                                   .::!!!!!!!:.
  .!!!!!:.                        .:!!!!!!!!!!!!
  ~~~~!!!!!!.                 .:!!!!!!!!!UWWW\$\$\$
      :\$\$NWX!!:           .:!!!!!!XUWW\$\$\$\$\$\$\$\$\$P
      \$\$\$\$\$##WX!:      .<!!!!UW\$\$\$\$"  \$\$\$\$\$\$\$\$#
      \$\$\$\$\$  \$\$\$UX   :!!UW\$\$\$\$\$\$\$\$\$   4\$\$\$\$\$*
      ^\$\$\$B  \$\$\$\$\\     \$\$\$\$\$\$\$\$\$\$\$\$   d\$\$R"
        "*\$bd\$\$\$\$      '*\$\$\$\$\$\$\$\$\$\$\$o+#"
             \"""\"          \"""\"""\"
```
"""
function eyes(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
    $thoughts
     $thoughts
                                   .::!!!!!!!:.
  .!!!!!:.                        .:!!!!!!!!!!!!
  ~~~~!!!!!!.                 .:!!!!!!!!!UWWW\$\$\$
      :\$\$NWX!!:           .:!!!!!!XUWW\$\$\$\$\$\$\$\$\$P
      \$\$\$\$\$##WX!:      .<!!!!UW\$\$\$\$"  \$\$\$\$\$\$\$\$#
      \$\$\$\$\$  \$\$\$UX   :!!UW\$\$\$\$\$\$\$\$\$   4\$\$\$\$\$*
      ^\$\$\$B  \$\$\$\$\\     \$\$\$\$\$\$\$\$\$\$\$\$   d\$\$R"
        "*\$bd\$\$\$\$      '*\$\$\$\$\$\$\$\$\$\$\$o+#"
             \"""\"          \"""\"""\"
"""
return the_cow
end


text:
module LogNormalTest

using Test
using ForneyLab
using ForneyLab: prod!, unsafeMean, unsafeLogMean, unsafeVar, unsafeLogVar, unsafeCov, unsafeLogCov, outboundType, isApplicable, dims
using ForneyLab: SPLogNormalOutNPP, VBLogNormalOut

@testset "dims" begin
    @test dims(ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=1.0)) == 1
end

@testset "vague" begin
    @test vague(LogNormal) == ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=huge)
end

@testset "unsafe mean and variance" begin
    @test unsafeMean(ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0)) == exp(1.0 + 0.5*2.0)
    @test unsafeLogMean(ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0)) == 1.0
    @test unsafeVar(ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0)) == exp(2.0*1.0 + 2.0)*(exp(2.0) - 1.0)
    @test unsafeLogVar(ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0)) == 2.0
    @test unsafeCov(ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0)) == exp(2.0*1.0 + 2.0)*(exp(2.0) - 1.0)
    @test unsafeLogCov(ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0)) == 2.0
end

@testset "log pdf" begin
    @test isapprox(logPdf(ProbabilityDistribution(Univariate, LogNormal, m=1.2, s=0.5),2), -1.522411904058978)
end

@testset "Gamma approximatons to LogNormal" begin
    @test ForneyLab.laplace(Gamma, ProbabilityDistribution(Univariate, LogNormal, m=0.0, s=2.0)) == ProbabilityDistribution(Univariate, Gamma, a=0.5, b=0.5)
end

@testset "prod!" begin
    @test ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0) * ProbabilityDistribution(Univariate, PointMass, m=1.0) == ProbabilityDistribution(Univariate, PointMass, m=1.0)
    @test_throws Exception ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0) * ProbabilityDistribution(Univariate, PointMass, m=-1.0)
    @test ProbabilityDistribution(Univariate, LogNormal, m=0.0, s=2.0) * ProbabilityDistribution(Univariate, Gamma, a=3.0, b=4.0) == ProbabilityDistribution(Univariate, Gamma, a=2.5, b=4.5)
end

#-------------
# Update rules
#-------------

@testset "SPLogNormalOutNPP" begin
    @test SPLogNormalOutNPP <: SumProductRule{LogNormal}
    @test outboundType(SPLogNormalOutNPP) == Message{LogNormal}
    @test isApplicable(SPLogNormalOutNPP, [Nothing, Message{PointMass}, Message{PointMass}])

    @test ruleSPLogNormalOutNPP(nothing, Message(Univariate, PointMass, m=1.0), Message(Univariate, PointMass, m=2.0)) == Message(Univariate, LogNormal, m=1.0, s=2.0)
end

@testset "VBLogNormalOut" begin
    @test VBLogNormalOut <: NaiveVariationalRule{LogNormal}
    @test outboundType(VBLogNormalOut) == Message{LogNormal}
    @test isApplicable(VBLogNormalOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBLogNormalOut, [ProbabilityDistribution, ProbabilityDistribution, Nothing])

    @test ruleVBLogNormalOut(nothing, ProbabilityDistribution(Univariate, PointMass, m=1.5), ProbabilityDistribution(Univariate, PointMass, m=3.0)) == Message(Univariate, LogNormal, m=1.5, s=3.0)
end

@testset "averageEnergy and differentialEntropy" begin
    @test differentialEntropy(ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0)) == averageEnergy(LogNormal, ProbabilityDistribution(Univariate, LogNormal, m=1.0, s=2.0), ProbabilityDistribution(Univariate, PointMass, m=1.0), ProbabilityDistribution(Univariate, PointMass, m=2.0))
end

end #module


text:
# Copyright 2015 Conrad Scott

module Topos

using DataFrames

export words, documents

include("dictionary.jl")
include("corpus.jl")

include("vocab.jl")
include("uci-bag-of-words.jl")

include("tfidf.jl")

include("utility.jl")

end


text:
const TEMP_REF = 20.0
const PRESSURE_REF = 1.0
const TEMP_REF_UNITFUL = TEMP_REF * u"°C"
export TEMP_REF, PRESSURE_REF

const DISPFORM_NAMES = ["Schott", "Sellmeier1", "Herzberger", "Sellmeier2", "Conrady", "Sellmeier3", "HandbookOfOptics1", "HandbookOfOptics2", "Sellmeier4", "Extended1", "Sellmeier5", "Extended2", "Extended3"]
const STATUS = ["Standard", "Preferred", "Obsolete", "Special", "Melt"]

const MIL_GLASSES = Dict{Int,Glass}()
const MODEL_GLASSES = Vector{Glass}(undef, 0)

"""
Special glass to represent air. Refractive index is defined to always be 1.0 for any temperature and pressure (other indices are relative to this).
"""
const Air = AirType()


text:
@testset "Statistics" begin

    fname = joinpath(dirname(@__FILE__), "..", "data", "test_Hz19.5-testing.bdf")

    s = read_SSR(fname)

    @testset "Global Field Power" begin

        ~ = gfp(s.data)  # TODO: Check result
    end
end


text:
# Starting point for dot operator sample for Learning Julia

# The dot operator applies an operation to a collection
arr1 = [2,4,6,8,10]
println(arr1)

result1 = arr1 .* 2 # multiply each element in the array by 2
println(result1)

# functions can also be applied to collections this way
result2 = sqrt.(arr1)
println(result2)

# custom functions can also be used this way
f(x) = 3x
result3 = f.(arr1)
println(result3)

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
using Pkg;
using Flux;
using Zygote;
using LinearAlgebra;
using Optim;
using Plots;
using Serialization;
using PyCall;
using Statistics;
using Printf;
#Pkg.activate("./");
#Pkg.instantiate();

function calcSymmetricTensorBasis!(dest, M,rho)
    numBasisElements, numSamples = size(dest);
    @assert (numBasisElements == 42)
    @assert (numSamples == size(M)[end] == size(rho)[end]);

    Threads.@threads for i in 1:numSamples
        id = Matrix(I,3,3);
        S = 0.5*(M[:,:,i]+M[:,:,i]');
        W = 0.5*(M[:,:,i]-M[:,:,i]');
        g = rho[:,:,i]'*rho[:,:,i];

        dest[1,i]  = S
        dest[2,i]  = g
        dest[3,i]  = S^2
        dest[4,i]  = W^2
        dest[5,i]  = g^2
        dest[6,i]  = S*W-W*S
        dest[7,i]  = S*g + g*S
        dest[8,i]  = g*W - W*g
        dest[9,i]  = S*W^2 + W^2*S
        dest[10,i] = S*g^2 + g^2*S
        dest[11,i] = S^2*W - W*S^2
        dest[12,i] = g^2*W - W*g^2
        dest[13,i] = g*S^2 + S^2*g
        dest[14,i] = g*W^2 + W^2*g
        dest[15,i] = S*W*g - g*W*S
        dest[16,i] = S*g*W - W*g*S
        dest[17,i] = g*S*W - W*S*g
        dest[18,i] = S^2*W^2 + W^2*S^2
        dest[19,i] = S^2*g^2 + g^2*S^2
        dest[20,i] = W^2*g^2 + g^2*W^2
        dest[21,i] = S^2*W*g - g*W*S^2
        dest[22,i] = S^2*g*W - W*g*S^2
        dest[23,i] = W^2*g*S + S*g*W^2
        dest[24,i] = g^2*S*W - W*S*g^2
        dest[25,i] = g^2*W*S - S*W*g^2
        dest[26,i] = S*W^2*g + g*W^2*S
        dest[27,i] = S*g^2*W - W*g^2*S
        dest[28,i] = W*S^2*g - g*S^2*W
        dest[29,i] = g*S*W^2 + W^2*S*g
        dest[30,i] = S*g^2*S^2 + S^2*g^2*S
        dest[31,i] = W*S^2*W^2 - W^2*S^2*W
        dest[32,i] = W*g^2*W^2 - W^2*g^2*W
        dest[33,i] = g*S^2*g^2 + g^2*S^2*g
        dest[34,i] = g*W^2*g^2 + g^2*W^2*g
        dest[35,i] = S^2*W^2*g + g*W^2*S^2
        dest[36,i] = S^2*g^2*W - W*g^2*S^2
        dest[37,i] = W^2*g^2*S + S*g^2*W^2
        dest[38,i] = g^2*S^2*W - W*S^2*g^2
        dest[39,i] = S*W^2*g^2 + g^2*W^2*S
        dest[40,i] = g*S^2*W^2 + W^2*S^2*g
        dest[41,i] = id;
    end

    return;
end

function calcYifengSymmetricTensorBasis!(dest, M,rho)
    numBasisElements, numSamples = size(dest);
    @assert (numSamples == size(M)[end] == size(rho)[end]);

    Threads.@threads for i in 1:numSamples
        id = Matrix(I,3,3);
        S = 0.5*(M[:,:,i]+M[:,:,i]');
        W = 0.5*(M[:,:,i]-M[:,:,i]');
        g = rho[:,:,i]'*rho[:,:,i];

        dest[1,i]  = S # TB#1
        dest[2,i]  = S*W-W*S               # TB#2
        dest[3,i]  = S^2 - (tr(S^2)/3)*id; # TB#3
        dest[4,i]  = W^2 - (tr(W^2)/3)*id; # TB#4
        dest[5,i]  = S^2*W - W*S^2                       # TB#5
        dest[6,i]  = S*W^2 + W^2*S - (2*tr(S*W^2)/3)*id; # TB#6
        dest[7,i]  = W*S*W^2 - W^2*S*W # TB#7 -- justification needed!
        dest[8,i]  = S*W*S^2 - S^2*W*S # TB#8 -- justification needed!
        dest[9,i]  = S^2*W^2 + W^2*S^2 - (2*tr(S^2*W^2)/3)*id; # TB#9
        dest[10,i] = W*S^2*W^2 - W^2*S^2*W  # TB#10

    end

    return;
end

function calcSkewSymmetricTensorBasis!(dest, M,rho)
    numBasisElements, numSamples = size(dest);
    @assert (numBasisElements == 17)
    @assert (numSamples == size(M)[end] == size(rho)[end]);

    Threads.@threads for i in 1:numSamples
        id = Matrix(I,3,3);
        S = 0.5*(M[:,:,i]+M[:,:,i]');
        W = 0.5*(M[:,:,i]-M[:,:,i]');
        g = rho[:,:,i]'*rho[:,:,i];

        dest[1,i]  = g * W + W * g
        dest[2,i]  = S^2 * W + W * S^2
        dest[3,i]  = g^2 * W + W * g^2
        dest[4,i]  = S * W*g + g * W * S
        dest[5,i]  = S * g*W + W * g * S
        dest[6,i]  = g * S*W + W * S * g
        dest[7,i]  = S^2 * W * g + g * W*S^2
        dest[8,i]  = S^2 * g + W * g * S^2
        dest[9,i]  = S^2 * g * W + W * g * S^2
        dest[10,i] = g^2 * S * W + W * S * g^2
        dest[11,i] = g^2 * W * S + S * W * g^2
        dest[12,i] = S * g^2 * W + W * g^2 * S
        dest[13,i] = W * S^2 * g + g * S^2 * W
        dest[14,i] = W * S^2 * W^2 + W^2 * S^2 * W
        dest[15,i] = W * g^2 * W^2 + W^2 * g^2 * W
        dest[16,i] = S^2 * g^2 * W + W * g^2 * S^2
        dest[17,i] = g^2 * S^2 * W + W * S^2 * g^2
    end

    return;
end

function calcTensorBasis!(mij,roij,symTB,skewSymTB)
    calcSymmetricTensorBasis!(symTB, mij, roij);
    calcSkewSymmetricTensorBasis!(skewSymTB, mij, roij);
end

function calcInvariants!(dest,S,W,g)
    numSamples = size(S)[end];

    #TODO: compare with calcYifengInvariants! and merge changes
    Threads.@threads for i in 1:numSamples
        dest[1,i]  = tr(S[i])
        dest[2,i]  = tr(S[i]^2)
        dest[3,i]  = tr(g[i])
        dest[4,i]  = tr(S[i] * g[i])
        dest[5,i]  = tr(g[i]^2)
        dest[6,i]  = tr(S[i] * g[i]^2)
        dest[7,i]  = tr(W[i])
        dest[8,i]  = tr(S[i] * W[i])
        dest[9,i]  = tr(g[i] * W[i])
        dest[10,i] = tr(S[i] * g[i]*W[i])
        dest[11,i] = tr(W[i]^2)
        dest[12,i] = tr(S[i] * W[i]^2)
        dest[13,i] = tr(g[i] * W[i]^2)
        dest[14,i] = tr(S[i] * g[i]*W[i]^2)
    end
end

function calcYifengInvariants!(dest,mij,roij)
    numInvariants, numSamples = size(dest);

    Threads.@threads for i in 1:numSamples
        S = 0.5*(mij[:,:,i]+mij[:,:,i]');
        W = 0.5*(mij[:,:,i]-mij[:,:,i]');
        dest[1,i]  = tr(S^2)      #λ_1
        dest[2,i]  = tr(W^2)      #λ_2
        dest[3,i]  = tr(S^3);     #λ_3
        dest[4,i] = tr(S * W^2)   #λ_4
        dest[5,i] = tr(W^2 * S^2) #λ_5
    end
end

#function TBNN_LearningPressure(miniBatchSize, maxiters)
    #--------------------begin problem setup--------------------#
    basepath = "/groups/chertkov/data/tetradData/InertialRangeData/";
    vgtDataPath = string(basepath,"aij_1024_gaus8_tetrad.bin");
    rhoDataPath = string(basepath,"roij_1024_gaus8_tetrad.bin");
    pressureDataPath = string(basepath, "pij_1024_gaus8_tetrad.bin");
    viscousDataPath = string(basepath, "vis_1024_gaus8_tetrad.bin");

    numTrajs = 32768;
    numTsteps = 1500;
    numTstepsToConsider = 200;
    numTstepsToSkip = 10;
    dt = 3e-4;

    numSymBasis = 42;
    numSkewSymBasis = 17;
    numInvariants = 14;

    if (wantToGenerateData == true)
        py"""
        import numpy as np
        mij = np.fromfile(file=$vgtDataPath, dtype=float)
        mij = mij.reshape([$numTrajs, $numTsteps, 3, 3])

        roij = np.fromfile(file=$rhoDataPath, dtype=float)
        roij = roij.reshape([$numTrajs, $numTsteps, 3, 3])

        pij = np.fromfile(file=$pressureDataPath, dtype=float)
        pij = pij.reshape([$numTrajs, $numTsteps, 3,3]);

        vij = np.fromfile(file=$viscousDataPath, dtype=float)
        vij = vij.reshape([$numTrajs, $numTsteps, 3,3]);

        #dmdt = [[(mij[i,j,:,:] - mij[i,j-1,:,:])/$dt for i in range($numTrajs)] for j in range(1,$numTsteps)];
        #subgrid = np.array([[dmdt[j][i] - (-pij[i,j,:,:] - np.matmul(mij[i,j,:,:],mij[i,j,:,:]) + vij[i,j,:,:]) for i in range($numTrajs)] for j in range(0:($numTsteps-1))])

        tsteps = range(0,$numTstepsToConsider,$numTstepsToSkip);

        mij = mij[:, tsteps, :, :];
        roij = roij[:, tsteps, :,:];
        pij = pij[:, tsteps, :, :];
        vij = vij[:, tsteps, :, :];
        #subgrid = subgrid[:, tsteps, :, :];
        """

        roij = PyArray(py"roij"o);
        mij = PyArray(py"mij"o);
        pij = PyArray(py"pij"o);
        visij = PyArray(py"vij"o);
        #subgridij = PyArray(py"subgrid"o); # Π_ij = dm/dt-(-p_{ij} - m_{ik}m_{kj} + v_{ij})

        roij = permutedims(roij, [3,4,2,1]);
        mij = permutedims(mij, [3,4,2,1]);
        pij = permutedims(pij, [3,4,2,1]);
        visij = permutedims(visij, [3,4,2,1]);
        #subgridij = permutedims(subgridij, [3,4,2,1]);

        numTsteps = size(mij)[3];

        mij = reshape(mij, (3,3,numTsteps*numTrajs));
        roij = reshape(roij, (3,3,numTsteps*numTrajs));
        pij = reshape(pij, (3,3, numTsteps*numTrajs));
        visij = reshape(visij, (3,3,numTsteps*numTrajs));
        #subgridij = reshape(subgridij, [3,3,numTsteps*numTrajs]);

        #generate learning inputs
        symTB = fill(zeros(3,3), (numSymBasis, numTsteps*numTrajs));
        skewSymTB = fill(zeros(3,3), (numSkewSymBasis, numTsteps*numTrajs));
        invariants = zeros(numInvariants, numTsteps*numTrajs);
        calcYifengSymmetricTensorBasis!(symTB, mij, roij);
        calcSkewSymmetricTensorBasis!(skewSymTB, mij, roij);
        # TODO this is garbage, skewSymTB[1,:] != W calcYifengInvariants!(invariants, symTB[1,:], skewSymTB[1,:], symTB[2,:]); #S,W,g

        x = invariants;
        numTrajs = size(x)[end];

        #learning "labels" are pressure hessian at each timestep/traj
        y = pij;

        normalizations = hcat(mean(x[1:numInvariants,:],dims=2),
                              var(x[1:numInvariants,:],dims=2));
    end

    # tsteps = Array{Float64, 1}(range(0.0, step=dt*numTstepsToSkip, length=numTsteps));
    # dt *= numTstepsToSkip;

    # miniBatchSize = 200;
    # maxiters = 20;
    #optimizer = Flux.Optimise.Optimiser(ExpDecay(0.001, 0.5, 1000, 1e-7), ADAM(0.001));
    optimizer = ADAM(0.0001);
    # see https://fluxml.ai/Flux.jl/stable/training/optimisers/

hiddenDim = 30;
θ, dnn_pressure = Flux.destructure(Chain(Dense(numInvariants,numSymBasis)));
#    θ, dnn_pressure = Flux.destructure(Chain(Dense(numInvariants,hiddenDim,tanh),
#                                             Dense(hiddenDim, hiddenDim, tanh),
#                                             Dense(hiddenDim,numSymBasis)));

    #u will be vgt (and g=ρ^Tρ eventually) in order to calculate tensor basis
    # we then transform to invariants, and pass that in to NN, recombining
    # just before we predict
    function predictPressureHessian(u,p,sampleIndxs)
        numSamples = size(u)[end];

        predCoefs = dnn_pressure(p)(u);

        intermediate = [[predCoefs[i,j]*symTB[i,sampleIndxs[j]] for i in 1:numSymBasis] for j in 1:numSamples];
        predictedPressureHessian = [sum(intermediate[j]) for j in 1:numSamples];

        return predictedPressureHessian;
    end

    function loss(features, labels, ps, sampleIndxs)
        numSamples = size(features)[end]
        pred = predictPressureHessian(features,
                                      ps,
                                      sampleIndxs);
        loss = [norm(y[:,:,i] - pred[i])^2 for i in 1:numSamples];
        return sum(loss)/numSamples;
    end

    function ψ(x)
       return log(1+x+x^2/2.0);
    end

    function dψdx(x)
       return (2(x+1))/(x^2 + 2x + 2);
    end

    function lossAndGrads(batchX, batchY, p, sampleIndxs)
        numSamples = size(batchX)[end];
        predCoefs = dnn_pressure(p)(batchX);
        intermediate = [[predCoefs[i,j]*symTB[i,sampleIndxs[j]] for i in 1:numSymBasis] for j in 1:numSamples];
        predictedPressureHessian = [sum(intermediate[j]) for j in 1:numSamples];
        δ = [(batchY[:,:,i] - predictedPressureHessian[i]) for i in 1:numSamples];
        contraction = hcat([[sum(δ[i].*symTB[j,i]) for i in 1:numSamples] for j in 1:numSymBasis]...);
        dnn_dθ(x) = pullback(ps->dnn_pressure(ps)(x), p)[2](1);
        grads = (2/numSamples)*Flux.pullback(ps->dnn_pressure(ps)(batchX),p)[2](contraction')[1]/numSamples;
        loss = sum([norm(δ[i])^2 for i in 1:numSamples])/numSamples;
        return loss,grads;
    end

    function truncatedLossAndGrads(batchX, batchY, p, sampleIndxs)
        numSamples = size(batchX)[end];
        predCoefs = dnn_pressure(p)(batchX);
        intermediate = [[predCoefs[i,j]*symTB[i,sampleIndxs[j]] for i in 1:numSymBasis] for j in 1:numSamples];
        predictedPressureHessian = [sum(intermediate[j]) for j in 1:numSamples];
        δ = [(batchY[:,:,i] - predictedPressureHessian[i]) for i in 1:numSamples];
        contraction = hcat([[-dψdx(sum(δ[i].*symTB[j,i])) for i in 1:numSamples] for j in 1:numSymBasis]...);
        grads = (2/numSamples)*pullback(ps->dnn_pressure(ps)(batchX),p)[2](contraction')[1]/numSamples;
        loss = sum([ψ(norm(δ[i])^2) for i in 1:numSamples])/numSamples;
        return loss,grads;
    end

    println("Setup complete, beginning learning");

function train(train_x,train_y, maxiters, miniBatchSize; seed = nothing)
    lossArray = zeros(maxiters);
    timeSinceLastUpdate = time();

    function customTrain!(ps, train_x, train_y, opt, sampleIdxs)
        l,gs = lossAndGrads(train_x, train_y, ps, sampleIdxs);
        ps += 0.1*gs;
        #Flux.update!(opt, ps, gs);
        return l;
    end

    if (seed != nothing)
        ps = seed;
    else
        ps = θ;
    end
    debugIters = 1;

    sampleIdxs = rand(1:numTrajs, miniBatchSize);
    for i in 1:maxiters
        lossArray[i] = customTrain!(ps, train_x[:,sampleIdxs], train_y[:,:,sampleIdxs], optimizer, sampleIdxs);
        if (i % debugIters == 0)
            @printf("%e   %d, %e\n", mean(lossArray[i-(debugIters-1):i]), i, (time()-timeSinceLastUpdate)/debugIters);
            timeSinceLastUpdate=time();
        end
    end
    return lossArray;
end

max_stars_repo_name: cmhyett/MLCoarseGrainedVGT

max_issues_repo_name: cmhyett/MLCoarseGrainedVGT

max_forks_repo_name: cmhyett/MLCoarseGrainedVGT


text:
#=
                       ------ MAIN INTERFACE -----
This code is a companion to the paper Surrogate Monte Carlo (SMC) -arxiv 2021. This is the main
code used to generate SMC data. However the default setting here are selected
for the code to run withing minutes (<10 min).
Authors: Andrei Da Silva & A.Christian Silva
e-mail: csilva@idatafactory.com
version on 10/2020
=#

using Random
using CSV
using DataFrames
using StatsBase

include("CDFJulia.jl")  # random numbers generation
include("utilities.jl") # opt functions
##
# Parameters
cgoal = 1e-3  # opt goal
nlag  = 50    # max lag for autocorrelation
n2lag = nlag  # max lag for autocorrelation of abs()
maxN  = 1e6   # max numb of terms in loop
##
# Data: downloaded from yahoo finance
sp = CSV.read("spyData.csv", DataFrame)
##
# Algorithm
# Step 1: Generate random numbers that have the same unconditional probability distribution
Random.seed!(28)
rs = getrnd(sp.R, length(sp.R))
#Step 2: Calculate the components of the optimization function
aR  = autocor(sp.R, 1:nlag)
bR  = autocor(abs.(sp.R), 1:n2lag)
cR  = crosscor(sp.R, abs.(sp.R), -nlag:nlag)
bR2 = autocor((sp.R).^2, 1:n2lag)
# Step 3: Calculate the optimization function. This is the starting point of the optimization
cmin = efun6(aR, bR, cR, bR2, rs, nlag, n2lag)
# Step 4: Perform optimization. Random permutation of rs should create final time series
# which minimizes the objective function and it is consequenty dynamically similar to the original
# time series
# ra: artificial time series
# E:  optimization function value
@time ra,E = annelOpt(aR,bR,cR,bR2,cmin,rs,nlag,n2lag,cgoal,maxN)
##
# How good is the opt ?
# Example:
#(1) Compare the autocorrelation of the abs(ra) with the "real" autocorrelation
B = autocor(abs.(ra), 1:n2lag)
plot(B)
plot!(bR)
#(2) Compare the crosscorrelation of the ra and abs(ra) with the measured crosscorrelation
C  = crosscor(ra, abs.(ra), -nlag:nlag)
plot(C)
plot!(cR)

max_stars_repo_name: silvaac/SMCarticle

max_issues_repo_name: silvaac/SMCarticle

max_forks_repo_name: silvaac/SMCarticle


text:
using Test
using jInvVisPyPlot
using jInv.Mesh
using PyPlot

# tests for regular mesh
domain = [1 2 0 .1 3 4]
n      = [81 70 30]
M2D    = getRegularMesh(domain[1:4],n[1:2])
m2D    = ones(tuple(M2D.n...))
M3D    = getRegularMesh(domain,n)
m3D    = ones(tuple(M3D.n...))

println("=== test plotModel (regular mesh) ===")
figure()
plotModel(m2D,true,M2D,2,[0.8,1.2]);

figure()
plotModel(m3D,true,M3D,2,[0.8,1.2]);

println("==== all passed === ")

max_stars_repo_name: JuliaInv/jInvVisPyPlot

max_issues_repo_name: JuliaInv/jInvVisPyPlot

max_forks_repo_name: JuliaInv/jInvVisPyPlot


text:
"""
    function affine_hit_and_run(
        warmup_points::Matrix{Float64},
        lbs::Vector{Float64},
        ubs::Vector{Float64};
        sample_iters = 100 .* (1:5),
        workers = [myid()],
        chains = length(workers),
        seed = rand(Int),
    )

Run a hit-and-run style sampling that starts from `warmup_points` and uses
their affine combinations for generating the run directions to sample the space
delimited by `lbs` and `ubs`.  The points that represent fluxes in
`warmup_points` should be organized in columns, i.e. `warmup_points[:,1]` is
the first warmup flux.

There are total `chains` of hit-and-run runs, each on a batch of
`size(warmup_points, 2)` points. The runs are scheduled on `workers`, for good
load balancing `chains` should be ideally much greater than `length(workers)`.

Each run continues for `maximum(sample_iters)` iterations; the numbers in
`sample_iters` represent the iterations at which the whole "current" batch of
points is collected for output. For example, `sample_iters=[1,4,5]` causes the
process run for 5 iterations, returning the sample batch that was produced by
1st, 4th and last (5th) iteration.

Returns a matrix of sampled fluxes (in columns), with all collected samples
horizontally concatenated. The total number of samples (columns) will be
`size(warmup_points,2) * chains * length(sample_iters)`.

# Example
```
using COBREXA
using Tulip

model = load_model(StandardModel, model_path)

warmup, lbs, ubs = warmup_from_variability(model, Tulip.Optimizer, 100)
samples = affine_hit_and_run(warmup, lbs, ubs, sample_iters = 1:3)
```
"""
function affine_hit_and_run(
    warmup_points::Matrix{Float64},
    lbs::Vector{Float64},
    ubs::Vector{Float64};
    sample_iters = 100 .* (1:5),
    workers = [myid()],
    chains = length(workers),
    seed = rand(Int),
)

    # distribute starting data to workers
    save_at.(workers, :cobrexa_hit_and_run_data, Ref((warmup_points, lbs, ubs)))

    # sample all chains
    samples = hcat(
        pmap(
            chain -> _affine_hit_and_run_chain(
                (@remote cobrexa_hit_and_run_data)...,
                sample_iters,
                seed + chain,
            ),
            CachingPool(workers),
            1:chains,
        )...,
    )

    # remove warmup points from workers
    map(fetch, remove_from.(workers, :cobrexa_hit_and_run_data))

    return samples
end

"""
    _affine_hit_and_run_chain(warmup, lbs, ubs, iters, seed)

Internal helper function for computing a single affine hit-and-run chain.
"""
function _affine_hit_and_run_chain(warmup, lbs, ubs, iters, seed)

    rng = StableRNG(seed % UInt)
    points = copy(warmup)
    d, n_points = size(points)
    result = Matrix{Float64}(undef, size(points, 1), 0)

    iter = 0

    for iter_target in iters

        while iter < iter_target
            iter += 1

            new_points = copy(points)

            for i = 1:n_points

                mix = rand(rng, n_points) .+ _constants.tolerance
                dir = points * (mix ./ sum(mix)) - points[:, i]

                # iteratively collect the maximum and minimum possible multiple
                # of `dir` added to the current point
                lambda_max = Inf
                lambda_min = -Inf
                for j = 1:d
                    dl = lbs[j] - points[j, i]
                    du = ubs[j] - points[j, i]
                    idir = 1 / dir[j]
                    if dir[j] < -_constants.tolerance
                        lower = du * idir
                        upper = dl * idir
                    elseif dir[j] > _constants.tolerance
                        lower = dl * idir
                        upper = du * idir
                    else
                        lower = -Inf
                        upper = Inf
                    end
                    lambda_min = max(lambda_min, lower)
                    lambda_max = min(lambda_max, upper)
                end

                lambda = lambda_min + rand(rng) * (lambda_max - lambda_min)
                !isfinite(lambda) && continue # avoid divergence
                new_points[:, i] = points[:, i] .+ lambda .* dir

                # TODO normally, here we would check if sum(S*new_point) is still
                # lower than the tolerance, but we shall trust the computer
                # instead.
            end

            points = new_points
        end

        result = hcat(result, points)
    end

    result
end


text:
using Test, UUIDs

using AnalyzeRegistry
using AnalyzeRegistry: parse_project

@testset "AnalyzeRegistry" begin
    general = general_registry()
    @test isdir(general)
    @test all(isdir, find_packages())
    # Test some properties of the `Measurements` package.  NOTE: they may change
    # in the future!
    measurements = analyze_from_registry(joinpath(general, "M", "Measurements"))
    @test measurements.uuid == UUID("eff96d63-e80a-5855-80a2-b1b0885c5ab7")
    @test measurements.reachable
    @test measurements.docs
    @test measurements.runtests
    @test !measurements.buildkite
    @test !isempty(measurements.lines_of_code)
    # Test results of a couple of packages.  Same caveat as above
    packages = [joinpath(general, p...) for p in (("C", "Cuba"), ("P", "PolynomialRoots"))]
    @test Set(packages) == Set(find_packages("Cuba", "PolynomialRoots")) == Set(find_packages(["Cuba", "PolynomialRoots"]))
    @test packages ⊆ find_packages()
    results = analyze_from_registry(packages)
    cuba, polyroots = results
    @test length(filter(p -> p.reachable, results)) == 2
    @test length(filter(p -> p.runtests, results)) == 2
    @test cuba.drone
    @test !polyroots.docs # Documentation is in the README!
    # We can also use broadcasting!
    @test Set(results) == Set(analyze_from_registry.(packages))

    # check `analyze_from_registry!` directly
    mktempdir() do root
        measurements2 = analyze_from_registry!(root, joinpath(general, "M", "Measurements"))
        @test measurements == measurements2
        @test isdir(joinpath(root, "eff96d63-e80a-5855-80a2-b1b0885c5ab7")) # not cleaned up yet
    end
end

@testset "`analyze`" begin
    pkg = analyze(pkgdir(AnalyzeRegistry))
    @test pkg.repo == "" # can't find repo from source code
    @test pkg.uuid == UUID("e713c705-17e4-4cec-abe0-95bf5bf3e10c")
    @test pkg.reachable == true # default
    @test pkg.docs == false
    @test pkg.runtests == true # here we are!
    @test pkg.github_actions == true
    @test length(pkg.license_files) == 1
    @test pkg.license_files[1].licenses_found == ["MIT"]
    @test pkg.license_files[1].license_filename == "LICENSE"
    @test pkg.license_files[1].license_file_percent_covered > 90
    @test pkg.license_files isa Vector{<:NamedTuple}
    @test keys(pkg.license_files[1]) == (:license_filename, :licenses_found, :license_file_percent_covered)
    @test isempty(pkg.licenses_in_project)
    @test !isempty(pkg.lines_of_code)
    @test pkg.lines_of_code isa Vector{<:NamedTuple}
    @test keys(pkg.lines_of_code[1]) == (:directory, :language, :sublanguage, :files, :code, :comments, :blanks)
    idx = findfirst(row -> row.directory=="src" && row.language==:Julia && row.sublanguage===nothing, pkg.lines_of_code)
    @test idx !== nothing
    @test pkg.lines_of_code[idx].code > 200

    # the tests folder isn't a package!
    # But this helps catch issues in error paths for when things go wrong
    bad_pkg = analyze(".")
    @test bad_pkg.repo == ""
    @test bad_pkg.uuid == UUID(UInt128(0))
    @test !bad_pkg.cirrus
    @test isempty(bad_pkg.license_files)
    @test isempty(bad_pkg.licenses_in_project)
end

@testset "`subdir` support" begin
    snoop_path = only(find_packages("SnoopCompileCore"))
    snoop = analyze_from_registry(snoop_path)
    @test !isempty(snoop.subdir)
    @test snoop.name == "SnoopCompileCore" # this test would fail if we were parsing the wrong Project.toml (that of SnoopCompile)
end

@testset "`parse_project`" begin
    bad_project = (; name="Invalid Project.toml", uuid=UUID(UInt128(0)), licenses_in_project=String[])
    # malformatted TOML file
    @test parse_project("missingquote") == bad_project

    # bad UUID
    @test parse_project("baduuid") == bad_project

    # non-existent folder
    @test parse_project("rstratarstra") == bad_project

    # proper Project.toml
    this_project = (; name = "AnalyzeRegistry", uuid = UUID("e713c705-17e4-4cec-abe0-95bf5bf3e10c"), licenses_in_project = String[])
    @test parse_project(joinpath(@__DIR__, "..")) == this_project

    # has `license = "MIT"`
    project_1 = (; name = "AnalyzeRegistry", uuid = UUID("e713c705-17e4-4cec-abe0-95bf5bf3e10c"), licenses_in_project=["MIT"])
    @test parse_project("license_in_project") == project_1

    # has `license = ["MIT", "GPL"]`
    project_2 = (; name = "AnalyzeRegistry", uuid = UUID("e713c705-17e4-4cec-abe0-95bf5bf3e10c"), licenses_in_project=["MIT", "GPL"])
    @test parse_project("licenses_in_project") == project_2
end

@testset "`show`" begin
    # this is mostly to test that `show` doesn't error
    str = sprint(show, analyze(pkgdir(AnalyzeRegistry)))
    @test occursin("* uuid: e713c705-17e4-4cec-abe0-95bf5bf3e10c", str)
    @test occursin("* OSI approved: true", str)
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.DLR_Visualization (Modia3D/renderer/DLR_Visualization/_module.jl)
#

@enum ShapeType box=1 sphere=2 cylinder=3 cone=4 capsule=5 coordsys=6 spring=7 gearwheel=8 pipe=9 grid=10 beam=11

const mvecSize  = MVector{3,Cdouble}(0.0,0.0,0.0)

@inline function visualizeShape(simVis::SimVis_Renderer, frame::Composition.Object3D, shape::ShapeType, id::Ptr{NOTHING}, material::Graphics.Material,
                                Lx::Float64, Ly::Float64, Lz::Float64, extras::MVector{3,Float64}=Basics.ZeroMVector)
  mvecSize[1] = Lx
  mvecSize[2] = Ly
  mvecSize[3] = Lz

  SimVis_setBaseObject(simVis, id, Cint(0), Cint(shape), frame.r_abs, frame.R_abs, mvecSize, material.color, Cint(material.wireframe),
                       Cint(material.reflectslight), material.shininess, extras, material.transparency, Cint(0), Cint(material.shadowMask))
end

visualizeSphere(data::Graphics.Sphere, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, sphere, id, data.material, data.Dx, data.Dx, data.Dx)
visualizeSolidSphere(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, sphere, id, data.material, data.geo.Dx, data.geo.Dx, data.geo.Dx)
#
visualizeEllipsoid(data::Graphics.Ellipsoid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, sphere, id, data.material, data.Dx, data.Dy, data.Dz)
visualizeSolidEllipsoid(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, sphere, id, data.material, data.geo.Dx, data.geo.Dy, data.geo.Dz)
#
visualizeBox(data::Graphics.Box, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, box, id, data.material, data.Lx, data.Ly, data.Lz)
visualizeSolidBox(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, box, id, data.material, data.geo.Lx, data.geo.Ly, data.geo.Lz)
#
visualizeCylinder(data::Graphics.Cylinder, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, cylinder, id, data.material, data.Dx, data.Dy, data.Lz)
visualizeSolidCylinder(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, cylinder, id, data.material, data.geo.Dx, data.geo.Dy, data.geo.Lz)
#
visualizeCone(data::Graphics.Cone, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, cone, id, data.material, data.Dx, data.Dy, data.Lz, data.extras)
visualizeSolidCone(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, cone, id, data.material, data.geo.Dx, data.geo.Dy, data.geo.Lz, data.geo.extras)
#
visualizeCapsule(data::Graphics.Capsule, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, capsule, id, data.material, data.Dx, data.Dy, data.Lz)
visualizeSolidCapsule(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, capsule, id, data.material, data.geo.Dx, data.geo.Dy, data.geo.Lz)
#
visualizeSpring(data::Graphics.Spring, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, spring, id, data.material, data.Dx, data.Dy, data.Lz, data.extras)
#
visualizeGearWheel(data::Graphics.GearWheel, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, gearwheel, id, data.material, data.Dx, data.Dy, data.Lz, data.extras)
#
visualizePipe(data::Graphics.Pipe, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, pipe, id, data.material, data.Dx, data.Dy, data.Lz, data.extras)
visualizeSolidPipe(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, pipe, id, data.material, data.geo.Dx, data.geo.Dy, data.geo.Lz, data.geo.extras)
#
visualizeBeam(data::Graphics.Beam, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, beam, id, data.material, data.Lx,data.Dy, data.Lz)
visualizeSolidBeam(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, beam, id, data.material, data.geo.Lx,data.geo.Dy, data.geo.Lz)
#
const DefaultMaterial = Graphics.Material()
visualizeCoordinateSystem(data::Graphics.CoordinateSystem, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, coordsys, id, DefaultMaterial, data.Lx, data.Lx, data.Lx)

#
visualizeGrid(data::Graphics.Grid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
          visualizeShape(simVis, frame, grid, id, DefaultMaterial, data.Lx, data.Ly, data.Lx, data.extras)
#
const emptyShaderName=" "
visualizeFileMesh(data::Graphics.FileMesh, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
     SimVis_setFileObject(simVis, id, Cint(0), frame.r_abs, frame.R_abs,
          data.scaleFactor, Cint(data.material.reflectslight), data.material.shininess, data.material.transparency, Cint(data.material.wireframe), Cint(0),
          data.filename, Cint(data.smoothNormals), data.useMaterialColor, MVector{3,Cint}(data.material.color),
          Cint(data.material.shadowMask), emptyShaderName)
#
visualizeSolidFileMesh(data::Solids.Solid, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
     SimVis_setFileObject(simVis, id, Cint(0), frame.r_abs, frame.R_abs,
          data.geo.scaleFactor, Cint(data.material.reflectslight), data.material.shininess, data.material.transparency, Cint(data.material.wireframe), Cint(0),
          data.geo.filename, Cint(data.geo.smoothNormals), data.geo.useGraphicsMaterialColor, MVector{3,Cint}(data.material.color),
          Cint(data.material.shadowMask), emptyShaderName)

#### Functions only operational for professional edition of SimVis

#
visualizeTextShape(data::Graphics.TextShape, frame::Composition.Object3D, id::Ptr{NOTHING}, simVis::SimVis_Renderer) =
    SimVis_setTextObject(simVis, id, Cint(data.axisAlignment), data.text, 0.0, Cint(0), frame.r_abs, frame.R_abs,
                         data.font.charSize, data.font.fontFileName, data.font.color, data.font.transparency,
                         data.offset, Cint(data.alignment), Cint(0))


text:
using Catlab.CategoricalAlgebra, Catlab.Present, Catlab.Graphs
using Nauty
import LightGraphs
const lg = LightGraphs
using DataStructures: OrderedSet

"""
The general `to_lg` function is overkill when we are just trying to compute the
canonical label of a (simple) DiGraph. This output can be used with Nauty.jl by
giving an initial partition of all vertices in the same set.
"""
function graph_to_lg(g::Graph)::lg.DiGraph
  res = lg.DiGraph(nv(g))
  for (s, t) in zip(g[:src], g[:tgt])
    lg.add_edge!(res, s, t)
  end
  return lg
end

"""
Produce a simple digraph and initial coloring that encodes an ACSet. I.e. a
faithful functor from ACSet -> SimpleColoredDiGraph (for any schema C)

The fact that different vertices come from different ACSet objects can be
encoded as an initial coloring.

How do we denote different morphisms between the same pair of objects?
- subdivide the graph, putting a vertex in the middle of every edge
- color the 'edge vertex' by which morphism it corresponds to.

For example, the multi-digraph V₁ ⟶ V₂ is translated to the vertex-labeled
digraph: V₁ <- src -- E₁ -- tgt -> V₂, where the V's have the same color, all
src nodes have the same color, etc.
"""
function to_lg(g::StructACSet, p::Presentation
  )::Tuple{lg.DiGraph, Vector{Int32}, Vector{Int32},
           Dict{Symbol, UnitRange}, Dict{Symbol, Vector{Any}}}
  colors, curr, oinds = Any[], 1, Dict{Symbol}{UnitRange}()
  lgr = lg.DiGraph(sum(g.obs))

  # Handle Ob
  for o in [g.args[1] for g in p.generators[:Ob]]
  append!(colors, repeat([o], nparts(g, o)))
  oinds[o] = curr:length(colors)
  curr = length(colors) + 1
  end

  # Handle Hom
  for h in p.generators[:Hom]
    hom_name, d_, cd_  = h.args
    domcodom = [d_.args[1], cd_.args[1]]
    for st in enumerate(g[hom_name])
      s_ind, t_ind = [oinds[name][val] for (name, val) in zip(domcodom, st)]
      lg.add_vertex!(lgr)
      push!(colors, hom_name)
      lg.add_edge!(lgr, s_ind, lg.nv(lgr))
      lg.add_edge!(lgr, lg.nv(lgr), t_ind)
    end
    oinds[hom_name] = curr:length(colors)
    curr = length(colors) + 1
  end

  # Handle Data
  attrdict = Dict{Pair{Symbol, Any},Int}() # store the vertex for each value
  attrvals = Dict{Symbol, Vector{Any}}()   # store all values in canonical order
  for dt in p.generators[:AttrType]
    d_name, vals = dt.args[1], Set()
    # Find out how many unique elements there are of this type, across all attrs
    for h in p.generators[:Attr]
      if h.args[3].args[1] == d_name
        union!(vals, g[h.args[1]])
      end
    end
    # Give each val a unique color and a unique vertex
    attrvals[d_name] = sort(collect(vals))
    for v in attrvals[d_name]
      lg.add_vertex!(lgr)
      col = d_name => v
      push!(colors, col)
      attrdict[col] = lg.nv(lgr)
    end
    oinds[d_name] = curr:length(colors)
    curr = length(colors) + 1
  end

  # Handle Attr
  for h in p.generators[:Attr]
    attr_name, d_, cd_  = h.args
    dom_,codom_ = [d_.args[1], cd_.args[1]]
    for (s, v) in enumerate(g[attr_name])
      s_ind = oinds[dom_][s]
      t_ind = attrdict[codom_ => string(v)]
      lg.add_vertex!(lgr)
      push!(colors, attr_name)
      lg.add_edge!(lgr, s_ind, lg.nv(lgr))
      lg.add_edge!(lgr, lg.nv(lgr), t_ind)
    end
    oinds[attr_name] = curr:length(colors)
    curr = length(colors) + 1
  end

  # Compute labels and partition from color data
  #---------------------------------------------
  # Get an array of arrays of nodes which are all the same color
  colorsarray = Vector{Int32}[findall(==(c), colors)
      for c in OrderedSet(colors)]
  # Assign numeric values corresponding to the colors, starting from '0'
  labels = vcat(((y -> map(x->x-1, y)).(colorsarray))...)
  # Encode that a partition ends with a '0'
  partition = vcat([begin z[end]=0; z end
                    for z in ones.(Cint,size.(colorsarray))]...)
  return (lgr, labels, partition, oinds, attrvals)
end

"""
Get canonical ACSet from Nauty's canong output
"""
function from_canong(g::Vector{UInt64}, p::Presentation, G::StructACSet
                    )::StructACSet
  _, _, _, oinds, attrdict = to_lg(G,p)
  gadj = Nauty.label_to_adj(g)

  # Each hom with the offset of its domain and codomain
  homdata = [let (hn, d, cd)=h.args;
             (hn, [oinds[x.args[1]][1]-1 for x in [d,cd]]) end
             for h in p.generators[:Hom]]
  # The ordered data for each hom
  homvals = [last.(sort([[findfirst(>(0), v) - offset for (offset, v)
          in zip(offsets, [gadj[:,i], gadj[i,:]])]
          for i in oinds[h]])) for (h,offsets) in homdata]
  # Each attr with the offset of its domain and codomain
  attrdata = [let (hn, d, cd)=h.args;
             (hn, cd.args[1], oinds[x.args[1]][1]-1 for x in [d,cd]) end
             for h in p.generators[:Attr]]
  # The ordered data for each attr
  attrvals = [last.(sort([[findfirst(>(0), v) - offset for (offset, v)
              in zip(offsets, [gadj[:,i], gadj[i,:]])]
              for i in oinds[h]])) for (h,_,offsets) in attrdata]
  I = deepcopy(G)
  for (h, vs) in zip(first.(homdata), homvals)
    set_subpart!(I, h, vs)
  end
  for ((h,cd), vs) in zip(attrdata, attrvals)
    set_subpart!(I, h, [attrdict[cd][v] for v in vs])
  end
  return I
end

"""Convert to Nauty.jl input, then interpret result back into Catlab language"""
function canonical_iso_nauty(g::StructACSet, p::Presentation)::StructACSet
  lgrph, lab, prt, _, _ = to_lg(g, p)
  cf = coloured_digraph_canonical_form(lgrph, lab, prt)
  return from_canong(cf, p, g)
end

"""Digraph isomorphism testing in Nauty.jl"""
function coloured_digraph_canonical_form(g, labels, partition)
  options = Nauty.doptionblk_mutable()
  options.defaultptn = false
  options.getcanon = true
  return Nauty.densenauty(Nauty.lg_to_nauty(g), options, labels, partition).canong
end


text:
using Test, LinearAlgebra, MPStates, Random

include("./tensor_operations_tests.jl")
include("./mps_tests.jl")
include("./mps_operations_tests.jl")
include("./mpo_tests.jl")
include("./mpo_operations_tests.jl")
include("./dmrg_tests.jl")


text:
@testset "Make sure that EV is filled" begin
    EV = evs.EV
    dEV = evs.dEV
    dEVσ = evs.dEVσ
    roy = 0.25
    geoid = 4.7
    itype = (geoid, roy,)

    # test that full VFI does something
    fill!(EV, 0.0)
    ShaleDrillingModel.solve_vf_terminal!(EV, dEV, dEVσ, wp)
    solve_vf_infill!(evs, tmpv, prim, θt, σv, false, itype)
    @test !all(EV .== 0.0)

    # test that full VFI does something and gets gradient
    fill!(EV, 0.0)
    fill!(dEV, 0.0)
    solve_vf_infill!(evs, tmpv, prim, θt, σv, true, itype)
    @test !all(EV .== 0.0)
    @test !all(dEV .== 0.0)
end

# ---------------- Regime 1 VFI and PFI - test gradient ------------------

@testset "Gradient of EV for infill" begin
    θ1 = similar(θt)
    θ2 = similar(θt)
    fdEV = similar(evs.dEV)
    roy = 0.25
    geoid = 4.7
    itype = (geoid, roy,)

    for k = 1:length(θt)
        h = peturb(θt[k])
        θ1 .= θt
        θ2 .= θt
        θ1[k] -= h
        θ2[k] += h
        hh = θ2[k] - θ1[k]

        solve_vf_terminal!(evs, prim)
        solve_vf_infill!(evs, tmpv, prim, θ1, σv, false, itype)
        fdEV[:,:,k,:] .= -evs.EV

        solve_vf_terminal!(evs, prim)
        solve_vf_infill!(evs, tmpv, prim, θ2, σv, false, itype)
        fdEV[:,:,k,:] .+= evs.EV
        fdEV[:,:,k,:] ./= hh
    end
    solve_vf_terminal!(evs, prim)
    solve_vf_infill!(evs, tmpv, prim, θt, σv, true, itype)

    @views maxv, idx = findmax(abs.(fdEV[:,:,2:end,:].-evs.dEV[:,:,2:end,:]))
    @views sub = CartesianIndices(fdEV[:,:,2:end,:])[idx]
    println("worst value is $maxv at $sub for dθ[2:end]")

    maxv, idx = findmax(abs.(fdEV.-evs.dEV))
    sub = CartesianIndices(fdEV)[idx]
    println("worst value is $maxv at $sub")
    # @test fdEV ≈ evs.dEV
    @test 0.0 < maxv < 1.0
    @test all(isfinite.(evs.dEV))
    @test all(isfinite.(fdEV))
end


text:
using HMMGradients
using BenchmarkTools
using LinearAlgebra
using SparseArrays
using Random

const SUITE = BenchmarkGroup()

sizes =  Dict()
sizes["small"] = ( 100, 10) #Nt,Ns,Ny
sizes["big"]   = (1000,100)

k = "forward!"
SUITE[k] = BenchmarkGroup(["forward!"])
for T in [Float32]
  for s in keys(sizes)
        SUITE[k][T, s] = @benchmarkable forward!(alpha,c,Nt,a,A,y) setup=begin
            Random.seed!(0)
            Nt, Ns = sizes[$s]
            A = rand($T, Ns, Ns)
            a = rand($T, Ns)
            y = randn($T, Nt, Ns)
            alpha = rand($T, Nt, Ns)
            c = randn($T, Nt)
        end
    end
end

k = "get_grad!"
SUITE[k] = BenchmarkGroup(["get_grad!"])
for T in [Float32]
  for s in keys(sizes)
        SUITE[k][T, s] = @benchmarkable HMMGradients.get_grad!(grada,gradA,grady,betat0,betat1,Nt,a,A,y,alpha,c) setup=begin
            Random.seed!(0)
            Nt, Ns = sizes[$s]
            A,gradA = rand($T, Ns, Ns),zeros($T, Ns, Ns)
            a,grada = rand($T, Ns),zeros($T, Ns)
            y,grady = rand($T, Nt, Ns),zeros($T, Nt, Ns)
            betat1 = ones($T,Ns)
            betat0 = zeros($T,Ns)
            alpha,c = forward(Nt,a,A,y)

        end
    end
end

k = "logforward!"
SUITE[k] = BenchmarkGroup(["forward!"])
for T in [Float32]
  for s in keys(sizes)
        SUITE[k][T, s] = @benchmarkable logforward!(alpha,Nt,a,A,y) setup=begin
            Random.seed!(0)
            Nt, Ns = sizes[$s]
            A = rand($T, Ns, Ns)
            a = rand($T, Ns)
            y = rand($T, Nt, Ns)
            alpha = rand($T, Nt, Ns)
        end
    end
end

k = "get_loggrad!"
SUITE[k] = BenchmarkGroup(["get_loggrad!"])
for T in [Float32]
  for s in keys(sizes)
        SUITE[k][T, s] = @benchmarkable HMMGradients.get_loggrad!(grada,gradA,grady,betat0,betat1,Nt,a,A,y,logalpha,logML) setup=begin
            Random.seed!(0)
            Nt, Ns = sizes[$s]
            A,gradA = rand($T, Ns, Ns),-$T(Inf).*ones($T, Ns, Ns)
            a,grada = rand($T, Ns),-$T(Inf).*ones($T, Ns)
            y,grady = rand($T, Nt, Ns),-$T(Inf).*ones($T, Nt, Ns)
            betat1 = zeros($T,Ns)
            betat0 = -$T(Inf).*ones($T,Ns)
            logalpha,logML = logforward(Nt,a,A,y)
        end
    end
end
#run(SUITE)