text:
using Bridge
using StaticArrays
import Bridge: b, σ, B, β, a, constdiff
const ℝ = SVector{N,T} where {N,T}
import Base.display
"""
    LotkaVolterraDiffusion <: ContinuousTimeProcess{ℝ{2}}
Struct defining a stochastic Lotka Volterra model
"""
struct LotkaVolterraDiffusion{T} <: ContinuousTimeProcess{ℝ{2,T}}
    α::T
    β::T
    γ::T
    δ::T
    σ1::T
    σ2::T
    function LotkaVolterraDiffusion(α::T, β, γ, δ, σ1, σ2) where T
        new{T}(α, β, γ, δ, σ1, σ2)
    end
end
function b(t, x, P::LotkaVolterraDiffusion{T}) where T
    ℝ{2}(P.α*x[1] - P.β*x[1]*x[2], P.δ*x[1]*x[2] - P.γ*x[2])
end
function σ(t, x, P::LotkaVolterraDiffusion{T}) where T
    SDiagonal(P.σ1, P.σ2)
end
constdiff(::LotkaVolterraDiffusion) = true
clone(P::LotkaVolterraDiffusion, θ) = LotkaVolterraDiffusion(θ...)
params(P::LotkaVolterraDiffusion) = [P.α, P.β, P.γ, P.δ, P.σ1, P.σ2]
domain(P::LotkaVolterraDiffusion) = LowerBoundedDomain((0.0, 0.0), (1,2))

# <---------------------------------------------
# this is optional, needed for conjugate updates

nonhypo(P::LotkaVolterraDiffusion, x) = x
@inline hypo_a_inv(P::LotkaVolterraDiffusion, t, x) = SDiagonal(1.0/P.σ1^2, 1.0/P.σ2^2)
num_non_hypo(P::Type{<:LotkaVolterraDiffusion}) = 2

phi(::Val{0}, t, x, P::LotkaVolterraDiffusion) = (zero(x[1]), zero(x[2]))
phi(::Val{1}, t, x, P::LotkaVolterraDiffusion) = (x[1], zero(x[2]))
phi(::Val{2}, t, x, P::LotkaVolterraDiffusion) = (-x[1]*x[2], zero(x[2]))
phi(::Val{3}, t, x, P::LotkaVolterraDiffusion) = (zero(x[1]), x[1]*x[2])
phi(::Val{4}, t, x, P::LotkaVolterraDiffusion) = (zero(x[1]), -x[2])
phi(::Val{5}, t, x, P::LotkaVolterraDiffusion) = (zero(x[1]), zero(x[2]))
phi(::Val{6}, t, x, P::LotkaVolterraDiffusion) = (zero(x[1]), zero(x[2]))

#
# <---------------------------------------------

struct LotkaVolterraDiffusionAux{R,S1,S2} <: ContinuousTimeProcess{ℝ{2,R}}
    α::R
    β::R
    γ::R
    δ::R
    σ1::R
    σ2::R
    t::Float64
    u::S1
    T::Float64
    v::S2

    function LotkaVolterraDiffusionAux(α::R, β, γ, δ, σ1, σ2, t::Float64, u::S1,
                                       T::Float64, v::S2) where {R,S1,S2}
        new{R,S1,S2}(α, β, γ, δ, σ1, σ2, t, u, T, v)
    end
end

function B(t, P::LotkaVolterraDiffusionAux)
#    ℝ{2}(P.α*x[1] - P.β*x[1]*x[2], P.δ*x[1]*x[2] - P.γ*x[2])
    @SMatrix [-0.0 -P.β*P.γ/P.δ; P.α*P.δ/P.β -0.0]
end

# mean = ℝ{2}(P.γ/P.δ, P.α/P.β)
function β(t, P::LotkaVolterraDiffusionAux)
    ℝ{2}(P.γ/P.δ*P.α, -P.α/P.β*P.γ)
end

function σ(t, P::LotkaVolterraDiffusionAux)
    SDiagonal(P.σ1, P.σ2)
end

function σ(t, x, P::LotkaVolterraDiffusionAux)
    SDiagonal(P.σ1, P.σ2)
end
depends_on_params(::LotkaVolterraDiffusionAux) = (1, 2, 3, 4, 5, 6)

constdiff(::LotkaVolterraDiffusionAux) = true
b(t, x, P::LotkaVolterraDiffusionAux) = B(t,P) * x + β(t,P)
a(t, P::LotkaVolterraDiffusionAux) = σ(t,P) * σ(t, P)'

clone(P::LotkaVolterraDiffusionAux, θ) = LotkaVolterraDiffusionAux(θ..., P.t,
                                                                   P.u, P.T, P.v)

clone(P::LotkaVolterraDiffusionAux, θ, v) = LotkaVolterraDiffusionAux(θ..., P.t,
                                                                      v, P.T, v)
params(P::LotkaVolterraDiffusionAux) = [P.α, P.β, P.γ, P.δ, P.σ1, P.σ2]


text:
using Measures

const CategoricalAesthetic =
    Union{Nothing, IndirectArray}

const NumericalAesthetic =
    Union{Nothing, AbstractMatrix, AbstractVector}

const NumericalOrCategoricalAesthetic =
    Union{CategoricalAesthetic, NumericalAesthetic}

@varset Aesthetics begin
    x,            Union{NumericalOrCategoricalAesthetic, Distribution}
    y,            Union{NumericalOrCategoricalAesthetic, Distribution}
    z,            Union{Nothing, Function, NumericalAesthetic}
    xend,         NumericalAesthetic
    yend,         NumericalAesthetic

    size,         Union{CategoricalAesthetic,Vector,Nothing}
    shape,        Union{CategoricalAesthetic,Vector,Nothing}
    color,        Union{CategoricalAesthetic,Vector,Nothing}
    alpha,        NumericalOrCategoricalAesthetic
    linestyle,    Union{CategoricalAesthetic,Vector,Nothing}

    label,        CategoricalAesthetic
    group,        CategoricalAesthetic

    xmin,         NumericalAesthetic
    xmax,         NumericalAesthetic
    ymin,         NumericalAesthetic
    ymax,         NumericalAesthetic

    # hexagon sizes used for hexbin
    xsize,        NumericalAesthetic
    ysize,        NumericalAesthetic

    # fixed lines
    xintercept,   NumericalAesthetic
    yintercept,   NumericalAesthetic
    intercept,    NumericalAesthetic
    slope,        NumericalAesthetic

    # boxplots
    middle,       NumericalAesthetic
    lower_hinge,  NumericalAesthetic
    upper_hinge,  NumericalAesthetic
    lower_fence,  NumericalAesthetic
    upper_fence,  NumericalAesthetic
    outliers,     NumericalAesthetic
    width,        NumericalAesthetic

    # subplots
    xgroup,       CategoricalAesthetic
    ygroup,       CategoricalAesthetic

    # guides
    xtick,        NumericalAesthetic
    ytick,        NumericalAesthetic
    xgrid,        NumericalAesthetic
    ygrid,        NumericalAesthetic
    color_key_colors,     Maybe(AbstractDict)
    color_key_title,      Maybe(AbstractString)
    color_key_continuous, Maybe(Bool)
    color_function,       Maybe(Function)
    titles,               Maybe(Dict{Symbol, AbstractString})
    shape_key_title,      Maybe(AbstractString)
    size_key_title,       Maybe(AbstractString)
    size_key_vals,        Maybe(AbstractDict)

    # mark some ticks as initially invisible
    xtickvisible,         Maybe(Vector{Bool})
    ytickvisible,         Maybe(Vector{Bool})

    # scale at which ticks should become visible
    xtickscale,           Maybe(Vector{Float64})
    ytickscale,           Maybe(Vector{Float64})

    # plot viewport extents
    xviewmin,     Any
    xviewmax,     Any
    yviewmin,     Any
    yviewmax,     Any

    # labeling functions
    x_label,      Function, showoff
    y_label,      Function, showoff
    xtick_label,  Function, showoff
    ytick_label,  Function, showoff
    color_label,  Function, showoff
    xgroup_label, Function, showoff
    ygroup_label, Function, showoff
    shape_label, Function, showoff
    size_label,   Function, showoff

    # pseudo-aesthetics
    pad_categorical_x, Union{Missing,Bool}, missing
    pad_categorical_y, Union{Missing,Bool}, missing
end

# Calculating fieldnames at runtime is expensive
const valid_aesthetics = fieldnames(Aesthetics)

function show(io::IO, data::Aesthetics)
    maxlen = 0
    print(io, "Aesthetics(")
    for name in valid_aesthetics
        val = getfield(data, name)
        if !ismissing(val) && val != nothing
            print(io, "\n  ", string(name), "=")
            show(io, getfield(data, name))
        end
    end
    print(io, "\n)\n")
end

# Alternate aesthetic names
const aesthetic_aliases =
    Dict{Symbol, Symbol}(:colour        => :color,
                                 :x_min         => :xmin,
                                 :x_max         => :xmax,
                                 :y_min         => :ymin,
                                 :y_max         => :ymax,
                                 :x_group       => :xgroup,
                                 :y_group       => :ygroup,
                                 :x_viewmin     => :xviewmin,
                                 :x_viewmax     => :xviewmax,
                                 :y_viewmin     => :yviewmin,
                                 :y_viewmax     => :yviewmax,
                                 :x_group_label => :xgroup_label,
                                 :y_group_label => :ygroup_label,
                                 :x_tick        => :xtick,
                                 :y_tick        => :ytick,
                                 :x_grid        => :xgrid,
                                 :y_grid        => :ygrid)

# Index as if this were a data frame
getindex(aes::Aesthetics, i::Integer, j::AbstractString) = getfield(aes, Symbol(j))[i]

# Return the set of variables that are non-nothing.
function defined_aesthetics(aes::Aesthetics)
    vars = Set{Symbol}()
    for name in valid_aesthetics
        getfield(aes, name) === nothing || push!(vars, name)
    end
    vars
end

# Checking aesthetics and giving reasonable error messages.

# Raise an error if any of the given aesthetics are not defined.
#
# Args:
#   who: A string naming the caller which is printed in the error message.
#   aes: An Aesthetics object.
#   vars: Symbol that must be defined in the aesthetics.
#
undefined_aesthetics(aes::Aesthetics, vars::Symbol...) =
        setdiff(Set(vars), defined_aesthetics(aes))

function assert_aesthetics_defined(who::AbstractString, aes::Aesthetics, vars::Symbol...)
    undefined_vars = undefined_aesthetics(aes, vars...)
    isempty(undefined_vars) || error("The following aesthetics are required by ",who,
            " but are not defined: ", join(undefined_vars,", "),"\n")
end

function assert_aesthetics_undefined(who::AbstractString, aes::Aesthetics, vars::Symbol...)
    defined_vars = intersect(Set(vars), defined_aesthetics(aes))
    isempty(defined_vars) || error("The following aesthetics are defined but incompatible with ",
            who,": ",join(defined_vars,", "),"\n")
end

function assert_aesthetics_equal_length(who::AbstractString, aes::Aesthetics, vars::Symbol...)
    defined_vars = filter(var -> !(getfield(aes, var) === nothing), [vars...])

    if !isempty(defined_vars)
        n = length(getfield(aes, first(defined_vars)))
        for var in defined_vars
            length(getfield(aes, var)) != n && error(
                    "The following aesthetics are required by ",who,
                    " to be of equal length: ",join(defined_vars,", "),"\n")
        end
    end
    nothing
end

# Replace values in a with non-nothing values in b.
#
# Args:
#   a: Destination.
#   b: Source.
#
# Returns: nothing
#
# Modifies: a
#
function update!(a::Aesthetics, b::Aesthetics)
    for name in valid_aesthetics
        issomething(getfield(b, name)) && setfield(a, name, getfield(b, name))
    end
    nothing
end

# Serialize aesthetics to JSON.

# Args:
#  a: aesthetics to serialize.
#
# Returns:
#   JSON data as a string.
#
json(a::Aesthetics) = join([string(a, ":", json(getfield(a, var))) for var in aes_vars], ",\n")

# Concatenate aesthetics.
#
# A new Aesthetics instance is produced with data vectors in each of the given
# Aesthetics concatenated, nothing being treated as an empty vector.
#
# Args:
#   aess: One or more aesthetics.
#
# Returns:
#   A new Aesthetics instance with vectors concatenated.
#
function concat(aess::Aesthetics...)
    cataes = Aesthetics()
    for aes in aess
        for var in valid_aesthetics
            if var in [:xviewmin, :yviewmin]
                mu, mv = getfield(cataes, var), getfield(aes, var)
                setfield!(cataes, var,
                          mu === nothing ? mv :
                             mv == nothing ? mu :
                                 min(mu, mv))
            elseif var in [:xviewmax, :yviewmax]
                mu, mv = getfield(cataes, var), getfield(aes, var)
                setfield!(cataes, var,
                          mu === nothing ? mv :
                             mv == nothing ? mu :
                                 max(mu, mv))
            else
                setfield!(cataes, var,
                          cat_aes_var!(getfield(cataes, var), getfield(aes, var)))
            end
        end
    end
    cataes
end

cat_aes_var!(a::(Nothing), b::(Nothing)) = a
cat_aes_var!(a::(Nothing), b::Union{Function,AbstractString}) = b
cat_aes_var!(a::(Nothing), b) = copy(b)
cat_aes_var!(a, b::(Nothing)) = a
cat_aes_var!(a::Function, b::Function) = a === string || a == showoff ? b : a

function cat_aes_var!(a::Dict, b::Dict)
    merge!(a, b)
    a
end

cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{T}) where {T <: Base.Callable} = append!(a, b)
cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{U}) where {T <: Base.Callable, U <: Base.Callable} =
        a = [a...,b...]

# Let arrays of numbers clobber arrays of functions. This is slightly odd
# behavior, comes up with with function statistics applied on a layer-wise
# basis.
cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{U}) where {T <: Base.Callable, U} = b
cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{U}) where {T, U <: Base.Callable} = a
cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{T}) where {T} = append!(a, b)
cat_aes_var!(a, b) = a

function cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{U}) where {T, U}
    V = promote_type(T, U)
    ab = Array{V}(undef, length(a) + length(b))
    i = 1
    for x in a
        ab[i] = x
        i += 1
    end
    for x in b
        ab[i] = x
        i += 1
    end

    return ab
end

function cat_aes_var!(xs::IndirectArray{T,1}, ys::IndirectArray{S,1}) where {T, S}
    TS = promote_type(T, S)
    return append!(IndirectArray(xs.index, convert(Array{TS},xs.values)),
                   IndirectArray(ys.index, convert(Array{TS},ys.values)))
end

cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{U}) where {T<:Measure, U<:Measure} = [a..., b...]

cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{U}) where {T<:Measure, U} =
    isabsolute(T) ? [a..., b...] : b

cat_aes_var!(a::AbstractVector{T}, b::AbstractVector{U}) where {T, U<:Measure} =
    isabsolute(U) ? a : [a..., b...]

# Summarizing aesthetics

# Produce a matrix of Aesthetic or Data objects partitioning the original
# Aesthetics or Data object by the cartesian product of xgroup and ygroup.
#
# This is useful primarily for drawing facets and subplots.
#
# Args:
#   aes: Aesthetics or Data objects to partition.
#
# Returns:
#   A Array{Aesthetics} of size max(1, length(xgroup)) by
#   max(1, length(ygroup))
#
function by_xy_group(aes::T, xgroup, ygroup,
                     num_xgroups, num_ygroups) where T <: Union{Data, Aesthetics}
    @assert xgroup === nothing || ygroup === nothing || length(xgroup) == length(ygroup)

    n = num_ygroups
    m = num_xgroups

    xrefs = xgroup === nothing ? [1] : xgroup
    yrefs = ygroup === nothing ? [1] : ygroup

    aes_grid = Array{T}(undef, n, m)
    staging = Array{AbstractArray}(undef, n, m)
    for i in 1:n, j in 1:m
        aes_grid[i, j] = T()
    end

    xgroup === nothing && ygroup === nothing && return aes_grid

    function make_pooled_array(::Type{IndirectArray{T,N,A,V}}, arr::AbstractArray) where {T,N,A,V}
        uarr = unique(arr)
        return IndirectArray(A(indexin(arr, uarr)), V(uarr))
    end
    make_pooled_array(::Type{IndirectArray{T,R,N,RA}},
            arr::IndirectArray{T,R,N,RA}) where {T,R,N,RA} = arr

    for var in fieldnames(T)
        # Skipped aesthetics. Don't try to partition aesthetics for which it
        # makes no sense to pass on to subplots.
        if var == :xgroup || var == :ygroup||
           var == :xtick || var == :ytick ||
           var == :xgrid || var == :ygrid ||
           var == :x_viewmin || var == :y_viewmin ||
           var == :x_viewmax || var == :y_viewmax ||
           var == :color_key_colors
            continue
        end

        vals = getfield(aes, var)
        if isa(vals, AbstractArray) && length(vals)>1
            if xgroup !== nothing && length(vals) !== length(xgroup) ||
               ygroup !== nothing && length(vals) !== length(ygroup)
                continue
            end

            for i in 1:n, j in 1:m
                staging[i, j] = similar(vals, 0)
            end

            for (i, j, v) in zip(cycle(yrefs), cycle(xrefs), vals)
                push!(staging[i, j], v)
            end

            for i in 1:n, j in 1:m
                if typeof(vals) <: IndirectArray
                    setfield!(aes_grid[i, j], var,
                              make_pooled_array(typeof(vals), staging[i, j]))
                else
                    if !applicable(convert, typeof(vals), staging[i, j])
                        T2 = eltype(vals)
                        if T2 <: Color T2 = Color end
                        da = Array{Union{Missing,T2}}(undef, length(staging[i, j]))
                        copy!(da, staging[i, j])
                        setfield!(aes_grid[i, j], var, da)
                    else
                        setfield!(aes_grid[i, j], var,
                                  convert(typeof(vals), copy(staging[i, j])))
                    end
                end
            end
        else
            for i in 1:n, j in 1:m
                setfield!(aes_grid[i, j], var, vals)
            end
        end
    end

    aes_grid
end

function inherit(a::Aesthetics, b::Aesthetics;
                 clobber=[])
    acopy = copy(a)
    inherit!(acopy, b, clobber=clobber)
    return acopy
end

function inherit!(a::Aesthetics, b::Aesthetics;
                  clobber=[])
    clobber_set = Set{Symbol}(clobber)
    for field in valid_aesthetics
        aval = getfield(a, field)
        bval = getfield(b, field)
        if field in clobber_set
            setfield!(a, field, bval)
        elseif aval === missing || aval === nothing || aval === string || aval == showoff
            setfield!(a, field, bval)
        elseif field == :xviewmin || field == :yviewmin
            bval != nothing && (aval == nothing || aval > bval) && setfield!(a, field, bval)
        elseif field == :xviewmax || field == :yviewmax
            bval != nothing && (aval == nothing || aval < bval) && setfield!(a, field, bval)
        elseif typeof(aval) <: Dict && typeof(bval) <: Dict
            merge!(aval, getfield(b, field))
        end
    end
    nothing
end

ext: jl


text:
module Orth

ACE.@extimports
ACE.@aceimports

# ------------ Products of Spherical Harmonics

import ACE.RPI.Rotations3D: ClebschGordan, cg_conditions
import ACE.SphericalHarmonics: SHBasis

ClebschGordan

struct SHProdCoeffs{T}
   cg::ClebschGordan{T}
end

SHProdCoeffs(T = Float64) = SHProdCoeffs(ClebschGordan(T))

struct SHProdFixlimi{T}
   l1::Int
   m1::Int
   l2::Int
   m2::Int
   M::Int
   cg::ClebschGordan{T}
end

(coeffs::SHProdCoeffs)(l1, m1, l2, m2) =
      SHProdFixlimi(l1, m1, l2, m2, m1+m2, coeffs.cg)

function Base.getindex(P::SHProdFixlimi{T}, L, M) where {T}
   if cg_conditions(P.l1, P.m1, P.l2, P.m2, L, M) &&
         iseven(P.l1 + P.l2 - L)   # condition on <l1 0 l2 0 | L 0>
      return ( sqrt( (2*P.l1+1)*(2*P.l2+1) / (4*pi*(2*L+1)) )
                 * P.cg(P.l1, P.m1, P.l2, P.m2, L, M)
                 * P.cg(P.l1, 0, P.l2, 0, L, 0) )
   else
      return zero(T)
   end
end

function Base.iterate(P::SHProdFixlimi)
   if !( (abs(P.m1) <= P.l1) && (abs(P.m2) <= P.l2) )  # all CG coeffs are zero
      return nothing
   end
   M = P.M
   L = max( abs(P.l1 - P.l2), abs(M) )  # L ≧ |l1-l2|
   return ( (L, M, P[L, M]), L )
end

function Base.iterate(P::SHProdFixlimi, state)
   L = state; M = P.M   # m1+m2=M is built-in
   L += 1
   if isodd(P.l1 + P.l2 - L)  # skip this iterate if <l1 0 l2 0 | L 0 > = 0
      return iterate(P, L)
   end
   if L > P.l1 + P.l2     # require L ≦ l1+l2, otherwise terminate
      return nothing
   end
   return ( (L, M, P[L, M]), L )
end

# ------------ Products of R_n Y_l^m  basis functions

import ACE.OrthPolys: OrthPolyProdCoeffs, EndlessVector

struct OnepProductCoeffs{T}
   shcoeffs::SHProdFixlimi{T}
   rcoeffs::OrthPolyProdCoeffs{T}
end

struct OnepProdCoeffMatrix{T}
   Psh::SHProdFixlimi{T}
   Pr::EndlessVector{T}
end

(coeffs::OnepProductCoeffs)(n1, l1, m1, n2, l2, m2) =
      OnepProdCoeffMatrix( coeffs.shcoeffs(l1, m1, l2, m2),
                           coeffs.rcoeffs(n1, n1) )

(P::OnepProdCoeffMatrix)(N, L, M) = P.Pr[N] * P.Psh[L, M]

# function Base.iterate(P::OnepProdCoeffMatrix)
#    L = abs(P.l1 - P.l2)
#    M = - L
#    return ( (L, M, P[L, M]), (L, M) )
# end
#
#
# function Base.iterate(P::SHProdFixlimi, state)
# end

end


text:
module MLJTSVDInterface

import TSVD
import MLJModelInterface
using Random: MersenneTwister, AbstractRNG, GLOBAL_RNG

const PKG = "TSVD"
const MMI = MLJModelInterface

"""
    TSVDTransformer()

Dimensionality reduction using truncated SVD.

This transformer performs linear dimensionality reduction by means of truncated singular value
decomposition (SVD). Contrary to PCA, this estimator does not center the data before computing
the singular value decomposition. This means it can work with sparse matrices efficiently.

"""
MMI.@mlj_model mutable struct TSVDTransformer <: MLJModelInterface.Unsupervised
    nvals::Int = 2
    maxiter::Int = 1000
    rng::Union{Int, AbstractRNG} = GLOBAL_RNG
end

struct TSVDTransformerResult
    singular_values::Vector{Float64}
    components::Matrix{Float64}
    is_table::Bool
end

as_matrix(X) = MMI.matrix(X)
as_matrix(X::AbstractArray) = X

_get_rng(rng::Int) = MersenneTwister(rng)
_get_rng(rng) = rng

function MMI.fit(transformer::TSVDTransformer, verbosity, Xuser)
    X = as_matrix(Xuser)
    rng = _get_rng(transformer.rng)

    U, s, V = TSVD.tsvd(
        X,
        transformer.nvals;
        maxiter=transformer.maxiter,
        initvec = convert(Vector{float(eltype(X))}, randn(rng, size(X,1)))
    )
    is_table = ~isa(Xuser, AbstractArray)
    fitresult = TSVDTransformerResult(s, V, is_table)
    cache = nothing

    return fitresult, cache, NamedTuple()
end

# for returning user-friendly form of the learned parameters:
function MMI.fitted_params(::TSVDTransformer, fitresult)
    singular_values = fitresult.singular_values
    components = fitresult.components
    is_table = fitresult.is_table
    return (singular_values = singular_values, components = components, is_table=is_table)
end

function MMI.transform(::TSVDTransformer, fitresult, Xuser)
    X = as_matrix(Xuser)
    Xtransformed = X * fitresult.components

    if fitresult.is_table
        Xtransformed = MMI.table(Xtransformed)
    end

    return Xtransformed
end

## META DATA

MMI.metadata_pkg(TSVDTransformer,
             name="$PKG",
             uuid="9449cd9e-2762-5aa3-a617-5413e99d722e",
             url="https://github.com/JuliaLinearAlgebra/TSVD.jl",
             is_pure_julia=true,
             license="MIT",
             is_wrapper=false
)

MMI.metadata_model(TSVDTransformer,
               input_scitype = Union{MMI.Table(MMI.Continuous),AbstractMatrix{MMI.Continuous}},
               output_scitype = Union{MMI.Table(MMI.Continuous),AbstractMatrix{MMI.Continuous}},
               docstring = "Truncated SVD dimensionality reduction",         # brief description
               path = "MLJTSVDInterface.TSVDTransformer"
               )

end # module


text:
#df,cf,bf,ins,ons=lps_prover("spec/fact.spec", 3,3)

df,cf,bf,eisf,insf,onsf=lps("spec/fact.jl", "spec/fact.spec", :fact, 3,3)
db,cb,bb,eisb,insb,onsb=lps("spec/fibo.jl", "spec/fibo.spec", :fibo, 4,4)
d1,c1,b1,eis1,ins1,ons1=lps("spec/id.jl", "spec/id.spec", :id, 3,3)
da,ca,ba,eisa,insa,onsa=lps("spec/add.jl", "spec/add.spec", :add, 4,4)

println("all specs ran")

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

# Single experiment, move to ensemble further on
# Some good parameter values are stored as comments right now
# because this is really good practice

using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using Plots
gr()

function corona!(du,u,p,t)
    S,E,I,R,N,D,C = u
    F, β0,α,κ,μ,σ,γ,d,λ = p
    dS = -β0*S*F/N - β(t,β0,D,N,κ,α)*S*I/N -μ*S # susceptible
    dE = β0*S*F/N + β(t,β0,D,N,κ,α)*S*I/N -(σ+μ)*E # exposed
    dI = σ*E - (γ+μ)*I # infected
    dR = γ*I - μ*R # removed (recovered + dead)
    dN = -μ*N # total population
    dD = d*γ*I - λ*D # severe, critical cases, and deaths
    dC = σ*E # +cumulative cases

    du[1] = dS; du[2] = dE; du[3] = dI; du[4] = dR
    du[5] = dN; du[6] = dD; du[7] = dC
end
β(t,β0,D,N,κ,α) = β0*(1-α)*(1-D/N)^κ
S0 = 14e6
u0 = [0.9*S0, 0.0, 0.0, 0.0, S0, 0.0, 0.0]
p_ = [10.0, 0.5944, 0.4239, 1117.3, 0.02, 1/3, 1/5,0.2, 1/11.2]
R0 = p_[2]/p_[7]*p_[6]/(p_[6]+p_[5])
tspan = (0.0, 21.0)
prob = ODEProblem(corona!, u0, tspan, p_)
solution = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 1)

tspan2 = (0.0,60.0)
prob = ODEProblem(corona!, u0, tspan2, p_)
solution_extrapolate = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 1)

# Ideal data
tsdata = Array(solution)
# Add noise to the data
noisy_data = tsdata + Float32(1e-5)*randn(eltype(tsdata), size(tsdata))

plot(abs.(tsdata-noisy_data)')

### Neural ODE

ann_node = FastChain(FastDense(7, 64, tanh),FastDense(64, 64, tanh), FastDense(64, 64, tanh), FastDense(64, 7))
p = Float64.(initial_params(ann_node))

function dudt_node(u,p,t)
    S,E,I,R,N,D,C = u
    F,β0,α,κ,μ,σ,γ,d,λ = p_
    dS,dE,dI,dR,dD = ann_node([S/N,E,I,R,N,D/N,C],p)

    dN = -μ*N # total population
    dC = σ*E # +cumulative cases

    [dS,dE,dI,dR,dN,dD,dC]
end
prob_node = ODEProblem(dudt_node, u0, tspan, p)
s = concrete_solve(prob_node, Tsit5(), u0, p, saveat = solution.t)

function predict(θ)
    Array(concrete_solve(prob_node, Vern7(), u0, θ, saveat = 1,
                         abstol=1e-6, reltol=1e-6,
                         sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

# No regularisation right now
function loss(θ)
    pred = predict(θ)
    sum(abs2, (noisy_data[2:4,:] .- pred[2:4,:])), pred # + 1e-5*sum(sum.(abs, params(ann)))
end

loss(p)

const losses = []
callback(θ,l,pred) = begin
    push!(losses, l)
    if length(losses)%50==0
        println(losses[end])
    end
    false
end

res1_node = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters = 500)
res2_node = DiffEqFlux.sciml_train(loss, res1_node.minimizer, BFGS(initial_stepnorm=0.01), cb=callback, maxiters = 10000)

prob_node2 = ODEProblem(dudt_node, u0, tspan, res2_node.minimizer)
s = solve(prob_node2, Tsit5(), saveat = 1)
scatter(solution, vars=[2,3,4], label=["True Exposed" "True Infected" "True Recovered"])
plot!(s, vars=[2,3,4], label=["Estimated Exposed" "Estimated Infected" "Estimated Recovered"])

# Plot the losses
plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss")

# Extrapolate out
prob_node_extrapolate = ODEProblem(dudt_node,u0, tspan2, res2_node.minimizer)
_sol_node = solve(prob_node_extrapolate, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 1)
p_node = scatter(solution_extrapolate, vars=[2,3,4], legend = :topleft, label=["True Exposed" "True Infected" "True Recovered"], title="Neural ODE Extrapolation")
plot!(p_node,_sol_node, lw=5, vars=[2,3,4], label=["Estimated Exposed" "Estimated Infected" "Estimated Recovered"])
plot!(p_node,[20.99,21.01],[0.0,maximum(hcat(Array(solution_extrapolate[2:4,:]),Array(_sol_node[2:4,:])))],lw=5,color=:black,label="Training Data End")

savefig("neuralode_extrapolation.png")
savefig("neuralode_extrapolation.pdf")

### Universal ODE Part 1

ann = FastChain(FastDense(3, 64, tanh),FastDense(64, 64, tanh), FastDense(64, 1))
p = Float64.(initial_params(ann))

function dudt_(u,p,t)
    S,E,I,R,N,D,C = u
    F, β0,α,κ,μ,σ,γ,d,λ = p_
    z = ann([S/N,I,D/N],p) # Exposure does not depend on exposed, removed, or cumulative!
    dS = -β0*S*F/N - z[1] -μ*S # susceptible
    dE = β0*S*F/N + z[1] -(σ+μ)*E # exposed
    dI = σ*E - (γ+μ)*I # infected
    dR = γ*I - μ*R # removed (recovered + dead)
    dN = -μ*N # total population
    dD = d*γ*I - λ*D # severe, critical cases, and deaths
    dC = σ*E # +cumulative cases

    [dS,dE,dI,dR,dN,dD,dC]
end
prob_nn = ODEProblem(dudt_,u0, tspan, p)
s = concrete_solve(prob_nn, Tsit5(), u0, p, saveat = 1)

plot(solution, vars=[2,3,4])
plot!(s[2:4,:]')

function predict(θ)
    Array(concrete_solve(prob_nn, Vern7(), u0, θ, saveat = solution.t,
                         abstol=1e-6, reltol=1e-6,
                         sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

# No regularisation right now
function loss(θ)
    pred = predict(θ)
    sum(abs2, noisy_data[2:4,:] .- pred[2:4,:]), pred # + 1e-5*sum(sum.(abs, params(ann)))
end

loss(p)

const losses = []
callback(θ,l,pred) = begin
    push!(losses, l)
    if length(losses)%50==0
        println(losses[end])
    end
    false
end

res1_uode = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters = 500)
res2_uode = DiffEqFlux.sciml_train(loss, res1_uode.minimizer, BFGS(initial_stepnorm=0.01), cb=callback, maxiters = 10000)

loss(res2_uode.minimizer)

prob_nn2 = ODEProblem(dudt_,u0, tspan, res2_uode.minimizer)
uode_sol = solve(prob_nn2, Tsit5(), saveat = 1)
plot(solution, vars=[2,3,4])
plot!(uode_sol, vars=[2,3,4])

# Plot the losses
plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss")

# Collect the state trajectory and the derivatives
X = noisy_data
# Ideal derivatives
DX = Array(solution(solution.t, Val{1}))

# Extrapolate out
prob_nn2 = ODEProblem(dudt_,u0, tspan2, res2_uode.minimizer)
_sol_uode = solve(prob_nn2, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 1)
p_uode = scatter(solution_extrapolate, vars=[2,3,4], legend = :topleft, label=["True Exposed" "True Infected" "True Recovered"], title="Universal ODE Extrapolation")
plot!(p_uode,_sol_uode, lw = 5, vars=[2,3,4], label=["Estimated Exposed" "Estimated Infected" "Estimated Recovered"])
plot!(p_uode,[20.99,21.01],[0.0,maximum(hcat(Array(solution_extrapolate[2:4,:]),Array(_sol_uode[2:4,:])))],lw=5,color=:black,label="Training Data End")

savefig("universalode_extrapolation.png")
savefig("universalode_extrapolation.pdf")

### Universal ODE Part 2: SInDy to Equations

# Create a Basis
@variables u[1:3]
# Lots of polynomials
polys = Operation[]
for i ∈ 0:2, j ∈ 0:2, k ∈ 0:2
    push!(polys, u[1]^i * u[2]^j * u[3]^k)
end

# And some other stuff
h = [cos.(u)...; sin.(u)...; unique(polys)...]
basis = Basis(h, u)

X = noisy_data
# Ideal derivatives
DX = Array(solution(solution.t, Val{1}))
S,E,I,R,N,D,C = eachrow(X)
F,β0,α,κ,μ,_,γ,d,λ = p_
L = β.(0:tspan[end],β0,D,N,κ,α).*S.*I./N
L̂ = vec(ann([S./N I D./N]',res2_uode.minimizer))
X̂ = [S./N I D./N]'

scatter(L,title="Estimated vs Expected Exposure Term",label="True Exposure")
plot!(L̂,label="Estimated Exposure")
savefig("estimated_exposure.png")
savefig("estimated_exposure.pdf")

# Create an optimizer for the SINDY problem
opt = SR3()
# Create the thresholds which should be used in the search process
thresholds = exp10.(-6:0.1:1)

# Test on original data and without further knowledge
Ψ_direct = SInDy(X[2:4, :], DX[2:4, :], basis, thresholds, opt = opt, maxiter = 50000) # Fail
println(Ψ_direct.basis)
# Test on ideal derivative data ( not available )
Ψ_ideal = SInDy(X[2:4, 5:end], L[5:end], basis, thresholds, opt = opt, maxiter = 50000) # Succeed
println(Ψ_ideal.basis)
# Test on uode derivative data
Ψ = SInDy(X̂[:, 2:end], L̂[2:end], basis, thresholds,  opt = opt, maxiter = 10000, normalize = true, denoise = true) # Succeed
println(Ψ.basis)

# Build a ODE for the estimated system
function approx(u,p,t)
    S,E,I,R,N,D,C = u
    F, β0,α,κ,μ,σ,γ,d,λ = p_
    z = Ψ([S/N,I,D/N]) # Exposure does not depend on exposed, removed, or cumulative!
    dS = -β0*S*F/N - z[1] -μ*S # susceptible
    dE = β0*S*F/N + z[1] -(σ+μ)*E # exposed
    dI = σ*E - (γ+μ)*I # infected
    dR = γ*I - μ*R # removed (recovered + dead)
    dN = -μ*N # total population
    dD = d*γ*I - λ*D # severe, critical cases, and deaths
    dC = σ*E # +cumulative cases

    [dS,dE,dI,dR,dN,dD,dC]
end

# Create the approximated problem and solution
a_prob = ODEProblem{false}(approx, u0, tspan2, p_)
a_solution = solve(a_prob, Tsit5())

p_uodesindy = scatter(solution_extrapolate, vars=[2,3,4], legend = :topleft, label=["True Exposed" "True Infected" "True Recovered"])
plot!(p_uodesindy,a_solution, lw = 5, vars=[2,3,4], label=["Estimated Exposed" "Estimated Infected" "Estimated Recovered"])
plot!(p_uodesindy,[20.99,21.01],[0.0,maximum(hcat(Array(solution_extrapolate[2:4,:]),Array(_sol_uode[2:4,:])))],lw=5,color=:black,label="Training Data End")

savefig("universalodesindy_extrapolation.png")
savefig("universalodesindy_extrapolation.pdf")


text:
# using CPLEX
# using GLPK
# using Cbc
# using Clp
# using JuMP

include("../src/utils/model.jl")
# include("../src/solvers/mip.jl")
include("../src/solvers/lp_timing.jl")
include("../src/solvers/earliest_timing.jl")

# EXEMPLE DE LIGNE DE COMMANDE TESTÉE :
# ./bin/seqata.jl val -i data/ampl/01.ampl -p 3,4,5,6,7,8,9,1,10,2 -x clp
#

# ===========
# Préparation des arguments
args = Args.parse_commandline()
const TEST_LEVEL = 1
Args.set(:level, TEST_LEVEL) # sera passé temporairement à 0
Args.set(:outdir, "/tmp") # equivalent à l'option -d/tmp

# ===========
# Construction de l'instance mini
inst = instance_build_mini10()
@test inst.nb_planes == 10

# ==========
# Ordre des avions à tester
planes_str = "p3,p4,p5,p6,p7,p8,p9,p1,p10,p2"
# planes = matchall(r"[\w]+", planes_str)
# planes = collect( (m.match for m = eachmatch(r"[\w]+", planes_str) ))
planes = [ m.match for m = eachmatch(r"[\w]+", planes_str) ]

test_name = "Création et résolution du timing d'une solution pour algo :earliest"
lg1(test_name, "... ")
@testset "$test_name" begin

      sol = Solution(inst, update=false, algo=:earliest)

      @test isa(sol.solver, EarliestTimingSolver)
      # Random.shuffle!(sol.planes); solve!(sol)
      shuffle!(sol, do_update=true)
      @test sol.cost > 1000.0  # cost >> 700 si mélange presque surement

      # On impose l'ordre des avions comme souhaité
      set_from_names!(sol, planes)

      # On doit retrouver cet ordre dans la solution
      @test join(get_names(sol), ",") == planes_str

      # On résout le sous-problème de timing de cette solution
      solve!(sol)
      # @test sol.cost == 700.0 # avec :lp
      @test sol.cost == 2830.0 # avec :earliest

end
ln1(" fait.")

# ====================================================================
test_name = "Création et résolution du timing d'une solution pour algo :lp avec :cplex"
print(test_name, "... ")
@testset "$test_name" begin

      Args.set("timing_algo_solver", :lp)
      # Args.set("external_mip_solver", :clp)
      # Args.set("planes", planes_str) # idem que l'option -p

      # sol = Solution(inst, update=false, algo=:lp)
      sol = Solution(inst, update=false)   # alpo=:lp par défaut
      @test isa(sol.solver, LpTimingSolver)

      # @test Symbol(lowercase(solver_name(sol.solver.model))) == :clp

      # On mélange les avions puis on met à jour la solution
      # print("Mélange et évaluation de la solution (cos>1000 ?)...")
      shuffle!(sol, do_update=true)
      @test sol.cost > 1000.0  # cost >> 700 car mélange presque surement

      # On impose l'ordre des avions comme souhaité
      set_from_names!(sol, planes)

      # On doit retrouver cet ordre dans la solution
      @test join(get_names(sol), ",") == planes_str

      # On résout le sous-problème de timing de cette solution
      solve!(sol)
      @test sol.cost == 700.0
end
ln1(" fait.")

# ====================================================================
# test_name = "Création et résolution du timing d'une solution pour algo :lp avec :clp"
lg1(test_name, "... ")
@testset "$test_name" begin
      # Args.set("external_mip_solver", :clp)
      sol = Solution(inst, update=false, algo=:lp)

      @test isa(sol.solver, LpTimingSolver)

      shuffle!(sol, do_update=true)
      @test sol.cost > 1000.0  # cost >> 700 si mélange presque surement

      # On impose l'ordre des avions comme souhaité
      set_from_names!(sol, planes)

      # On doit retrouver cet ordre dans la solution
      @test join(get_names(sol), ",") == planes_str

      # On résout le sous-problème de timing de cette solution
      solve!(sol)
      @test sol.cost == 700.0

end
ln1(" fait.\n")


text:
@testset "Equivalenz TwoLevelChain - DissipativeLambaChain" begin
	# ----------------------------------------
	# At zero γd both models must coincide
	Γ = rand()

  problemA = WaveguideProblem(TwoLevelChain(1, Γ = Γ),
		WavePacket(SoftBoxMode(τ = 1.0), Coherent(1.0)), 2.
	)
  solA     = solve(problemA, reltol = 1e-10)[2]

  problemB = WaveguideProblem(DissipativeLambdaChain(1, Γ = Γ, γd = 0.0),
		WavePacket(SoftBoxMode(τ = 1.0), Coherent(1.0)), 2.
	)
  solB     = solve(problemB, reltol = 1e-10)[2]

  # Project on the groundstate
  σ_GG_A = transition(NLevelBasis(2), 1, 1)
  σ_GG_B = transition(NLevelBasis(3), 1, 1)
  @test sum(abs, expect(σ_GG_A, solA) .- expect(σ_GG_B, solB) ) < 1e-5

  # Project on the excited state
  σ_WW_A = transition(NLevelBasis(2), 2, 2)
  σ_WW_B = transition(NLevelBasis(3), 2, 2)
  @test sum(abs, expect(σ_WW_A, solA) .- expect(σ_WW_B, solB) ) < 1e-5
end

@testset "TwoLevelAtoms - Chirality" begin
	# ----------------------------------------
	# The dynamics of an atom is not impacted by latter atoms
  Γ = rand()
  σ_WW_i = transition(NLevelBasis(2), 2, 2)

  problem1 = WaveguideProblem(
    TwoLevelChain(1, Γ = Γ), ContinuousWave(1.0), 3.0
  )
  sol1 = solve(problem1, reltol = 1e-10)[2]

  problem2 = WaveguideProblem(
    TwoLevelChain(2, Γ = Γ), ContinuousWave(1.0), 3.0
  )
  sol2 = solve(problem2, reltol = 1e-10)[2]
	sol21, sol22 = ptrace.(sol2, 2), ptrace.(sol2, 1)

  problem3 = WaveguideProblem(
    TwoLevelChain(3, Γ = Γ), ContinuousWave(1.0), 3.0
  )
  sol3 = solve(problem3, reltol = 1e-10)[2]
	sol31, sol32 = [ptrace(sol, [2, 3]) for sol ∈ sol3], [ptrace(sol, [1, 3]) for sol ∈ sol3]

	@test sum(abs, tracedistance.( sol1,  sol21 )) < 1e-5
	@test sum(abs, tracedistance.( sol1,  sol31 )) < 1e-5
	@test sum(abs, tracedistance.( sol22, sol32 )) < 1e-5
end

@testset "DissipativeLambdaChain - Chirality" begin
	# ----------------------------------------
	# The dynamics of an atom is not impacted by latter atoms
  Γ, γd = rand(), rand()
  σ_WW_i = transition(NLevelBasis(3), 2, 2)

  problem1 = WaveguideProblem(
    DissipativeLambdaChain(1, γd = γd, Γ = Γ), ContinuousWave(1.0), 3.0
  )
  sol1 = solve(problem1, reltol = 1e-10)[2]

  problem2 = WaveguideProblem(
    DissipativeLambdaChain(2, γd = γd, Γ = Γ), ContinuousWave(1.0), 3.0
  )
  sol2 = solve(problem2, reltol = 1e-10)[2]
	sol21, sol22 = ptrace.(sol2, 2), ptrace.(sol2, 1)

  problem3 = WaveguideProblem(
    DissipativeLambdaChain(3, γd = γd, Γ = Γ), ContinuousWave(1.0), 3.0
  )
  sol3 = solve(problem3, reltol = 1e-10)[2]
	sol31, sol32 = [ptrace(sol, [2, 3]) for sol ∈ sol3], [ptrace(sol, [1, 3]) for sol ∈ sol3]

  @test sum(abs, tracedistance.( sol1,  sol21 )) < 1e-5
	@test sum(abs, tracedistance.( sol1,  sol31 )) < 1e-5
	@test sum(abs, tracedistance.( sol22, sol32 )) < 1e-5
end


text:
using DelaySSAToolkit, Catalyst
using Test
rn = @reaction_network begin
    1/(1+Y^2), 0 --> X
    1/(1+Y),   Y --> 0
end
jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)
# states(rn)
u0 = [0,0]
de_chan0 = [[]]
tf = 400.
tspan = (0,tf)
τ = 20.
dprob = DiscreteProblem(jumpsys, u0, tspan)
# jumpsys.dep_graph

delay_trigger_affect! = function (integrator, rng)
  append!(integrator.de_chan[1], τ)
end
delay_trigger = Dict(1=>delay_trigger_affect!)
delay_complete = Dict(1=>[2=>1, 1=>-1], 2=>[2=>1])
delay_interrupt = Dict()
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)

meanlist=[]
medianlist=[]
minlist=[]
maxlist=[]
stdlist=[]
algo_list =[DelayDirect(),DelayRejection(),DelayMNRM(),DelayDirectCR()]
for algo in algo_list
    djprob = DelayJumpProblem(jumpsys, dprob,algo,delayjumpset, de_chan0, save_positions=(false,false))
    a=@benchmark solve(djprob, SSAStepper())
    push!(meanlist,mean(a).time/1e9)
    push!(medianlist,median(a).time/1e9)
    push!(minlist,minimum(a).time/1e9)
    push!(maxlist,maximum(a).time/1e9)
    push!(stdlist,std(a).time/1e9)
    print(djprob.aggregator," OK","\n")
end

df=DataFrame(mean=meanlist,median=medianlist,min=minlist,max=maxlist,std=stdlist)
CSV.write("test/zxytest/results/oscillation_save_positions=f.csv",df)

df=CSV.read("test/zxytest/results/oscillation_save_positions=f.csv",DataFrame)
medianlist=df.median

medianvalue=[string(round(mt,digits=6),"s") for mt in medianlist]

algo_name = ["DelayDirect","DelayRejection","DelayMNRM","DelayDirectCR"]
using Plots
p1=bar(algo_name,medianlist,legend=:false,title="oscillation",ylabel="mediantime")
scatter!(algo_name, 0.0003 .+ medianlist , markeralpha=0, series_annotations=medianvalue)
savefig(p1,"test/zxytest/results/oscillation_median.png")


text:
# S
#imple policy that always allocates to the local pool.
struct LocalTracker
    threshold::Float64
end
LocalTracker() = LocalTracker(0.8)

# Ignore all hints.
Base.push!(::LocalTracker, block, x...) = nothing
Base.delete!(::LocalTracker, block) = nothing
update(::LocalTracker, _) = nothing
softevict!(::LocalTracker, manager, block) = nothing
prefetch!(::Block, ::LocalTracker, manager; kw...) = nothing
evict!(A, ::LocalTracker, manager) = nothing

# Only allocate from local pool.
function policy_new_alloc(t::LocalTracker, manager, bytes, id, _::AllocationPriority)
    # See if it's time to GC.
    allocated, total = getstate(getheap(manager, LocalPool()))
    if allocated / total > t.threshold
        println("Manually invoking the garbage collector!")
        @show manager
        GC.gc(true)
    end
    @return_if_exists unsafe_alloc_direct(LocalPool(), manager, bytes, id)
    return nothing
end

function defrag!(manager, policy::LocalTracker)
    cb = () -> unsafe_cleanup!(manager)
    defrag!(getheap(manager, LocalPool()); queued_callback = cb) do _, newblock, oldblock
        primary = getprimary(manager, oldblock)
        @assert primary === oldblock
        _ = unsafe_setprimary!(manager, oldblock, newblock; unsafe = true)
        return nothing
    end
    return nothing
end


text:
module Redpitaya

using Sockets

import Base: reset

export RedPitaya, receive, query, start, stop, send,
       receiveBinaryFloat,receiveBinaryInt16,receiveASCIIArray

mutable struct RedPitaya
  delim::String
  socket::TCPSocket

  function RedPitaya(host, port=5000)
    socket = connect(host,port)
    return new("\r\n", socket)
  end
end

"""
Send a command to the RedPitaya
"""
function send(rp::RedPitaya,cmd::String)
  write(rp.socket,cmd*rp.delim)
end

"""
Receive a String from the RedPitaya
"""
function receive(rp::RedPitaya)
  readline(rp.socket)#[1:end-2]
end

function readBinaryHeader(rp::RedPitaya)
  n = read!(rp.socket,Array{UInt8}(undef,1))
  n = read!(rp.socket,Array{UInt8}(undef,1))
  m = parse(Int64,String(copy(n)))
  numBytesStr = read!(rp.socket,Array{UInt8}(undef,m))
  numBytes = parse(Int64,String(copy(numBytesStr)))
  return numBytes
end

function receiveBinaryFloat(rp::RedPitaya)
  numBytes = readBinaryHeader(rp)
  u = read!(rp.socket, Array{Float32}(undef, div(numBytes,4)))
  uFl = [bswap(a) for a in u]
end

function receiveBinaryInt16(rp::RedPitaya)
  numBytes = readBinaryHeader(rp)
  u = read!(rp.socket, Array{Int16}(undef, div(numBytes,2)))
  uFl = [Float32(bswap(a)) for a in u]
end

function receiveASCIIArray(rp::RedPitaya)
  u = receive(rp)
  uFl = [parse(Float64,o) for o in split(u[2:end-1],",")]
end

"""
Perform a query with the RedPitaya. Return String
"""
function query(rp::RedPitaya,cmd::String)
  send(rp,cmd)
  return receive(rp)
end

"""
Perform a query with the RedPitaya. Parse result as type T
"""
function query(rp::RedPitaya,cmd::String,T::Type)
  a = query(rp,cmd)
  return parse(T,a)
end

include("IO.jl")
include("Generator.jl")
include("Acquire.jl")
include("Utils.jl")
include("Measurements.jl")
include("Custom.jl")

end #module


text:
include("split.jl")
include("apply.jl")
include("combine.jl")
include("timearray.jl")
include("readwrite.jl")
include("timeseriesrc.jl")


text:
"""
    function default()

The original cowfile

# Example

```jldoctest
julia> cowsay("Moo", cow=Cowsay.default)
 _____
< Moo >
 -----
        \\   ^__^
         \\  (oo)\\_______
            (__)\\       )\\/\\
                ||----w |
                ||     ||

```
"""
function default(;eyes="oo", tongue="  ", thoughts="\\")
thecow =
"""
        $thoughts   ^__^
         $thoughts  ($eyes)\\_______
            (__)\\       )\\/\\
             $tongue ||----w |
                ||     ||
"""
return thecow
end


text:
using SciMLTutorials
tutorials_dir = joinpath(dirname(@__DIR__), "tutorials")
SciMLTutorials.weave_file(joinpath(tutorials_dir, "Testing"), "test.jmd")


text:
using DrWatson
@quickactivate
using BSON
using DataFrames
using ValueHistories
using Statistics
using EvalMetrics
using CSV

function compute_score(info, score="AUC")
	(scores, labels) = (info[:val_scores], info[:val_labels])
	if score == "AUC"
		roc = EvalMetrics.roccurve(labels, scores)
		auc = EvalMetrics.auc_trapezoidal(roc...)
		return auc
	elseif score == "AUPRC"
		prc = EvalMetrics.prcurve(labels, scores)
		auprc = EvalMetrics.auc_trapezoidal(prc...)
		return auprc
	elseif score == "TPR@5"
		t5 = EvalMetrics.threshold_at_fpr(labels, scores, 0.05)
		cm5 = ConfusionMatrix(labels, scores, t5)
		tpr5 = EvalMetrics.true_positive_rate(cm5)
		return tpr5
	elseif score == "F1@5"
		t5 = EvalMetrics.threshold_at_fpr(labels, scores, 0.05)
		cm5 = ConfusionMatrix(labels, scores, t5)
		f5 = EvalMetrics.f1_score(cm5)
		return f5
	elseif score == "LOSS"
		try
			return minimum(get(info[:history][:validation_likelihood])[2])
		catch e
			try
				return minimum(get(info[:history][:validation_loss])[2])
			catch ee
				@error "no validation likelihood or validation loss"
			end
		end
	else
		@error "unknown score"
	end
end

function fix_info_name(name, score="reconstruction")
	spl = split(name, "lr=0_")
	name = (length(spl)==2) ? spl[1]*"lr=0.0001_"*spl[2] : name
	params = DrWatson.parse_savename("_"*name)[2]
	params = merge(params, Dict("score"=>score))
	name = DrWatson.savename(params, digits=5)*".bson"
	return name
end

function models_and_params(path_to_model)
	directories = []
	models = []

	for (root, dirs, files) in walkdir(path_to_model)
		for file in files
			push!(directories, root)
		end
	end

	for dir in unique(directories)
		fs = filter(x->(startswith(x, "model")), readdir(dir))
		info = map(x->x[7:end], fs)
		info = fix_info_name.(info) # model_name -> info name
		#par = map(x -> DrWatson.parse_savename("_"*x)[2], info)
		!isempty(fs) ? push!(models, (dir, fs, info)) : nothing
	end
	return models
end

function create_df(models; score::String="LOSS", images::Bool=true)

	if images
		df = DataFrame(
			path = String[], #path to model / encodings
			params = String[],
			info_path = String[],
			dataset = String[],
			ac = Int64[],
			seed = Int64[],
			criterion = Float32[]
		)
	else
		df = DataFrame(
			path = String[], #path to model / encodings
			params = String[],
			info_path = String[],
			dataset = String[],
			seed = Int64[],
			criterion = Float32[]
		)
	end

	i = 1

	for model in models
		(root, mod_path, infos) = model
		for (mod,info) in zip(mod_path, infos)
			try
				path = joinpath(root, info)
				info = BSON.load(path)
				if images
					update = [
						joinpath(root,mod),
						string(info[:parameters]),
						path,
						info[:dataset], # category <- mvtec
						info[:anomaly_class], # 1 <- mvtec
						info[:seed],
						compute_score(info, score)
					]
				else
					update = [
						joinpath(root,mod),
						string(info[:parameters]),
						path,
						info[:dataset],
						info[:seed],
						compute_score(info, score)
					]
				end
				push!(df, update)
			catch e
				println("error $(e)")
				println("info not found #$(i)")
				i += 1
			end
		end
	end
	return df
end

function return_best_n(df, rev=false, n=10)
	# filtering
	df = df[map(x->(sum(df.params .== x)==5), df.params),:]
	#df = df[df.dataset .== dataset, :]
	dff = sort(by(df, [:params, :dataset, :ac], :criterion => mean), :criterion_mean, rev=rev)
	df_top = []
	for dataset in unique(dff.dataset)
		for ac in 1:max(dff.ac...)
			top = first(dff[(dff.dataset .== dataset) .& (dff.ac .== ac), :], n)
			for i=1:n
				tmp = df[df.params .== top.params[i],:]
				ind = i.*ones(size(tmp,1))
				push!(df_top, hcat(tmp, DataFrame(ind=ind)))
			end
		end
	end
	return vcat(df_top...)
end

function return_best_n1(df, rev=false, n=10)
	#df = df[df.dataset .== dataset, :]
	dff = sort(by(df, [:params, :dataset], :criterion => mean), :criterion_mean, rev=rev)
	df_top = []
	for dataset in unique(dff.dataset)
		top = first(dff[dff.dataset .== dataset, :], n)
		for i=1:n
			tmp = df[df.params .== top.params[i],:]
			ind = i.*ones(size(tmp,1))
			push!(df_top, hcat(tmp, DataFrame(ind=ind)))
		end
	end
	return vcat(df_top...)
end

#path = "/home/skvarvit/generativead/GenerativeAD.jl/data/experiments/images_mnistc/vae/"
#models = models_and_params(path)
#df = create_df(models, images=true)
#print(df)
#CSV.write("vae_tab.csv", df)

cheetsheet = Dict("images_mnistc"=> true, "images_mvtec"=> true)

encoders =["vae"]  # ["vae","wae", "wae_vamp"]

for type in ["images_mvtec"] #["images", "tabular"]
	for encoder in encoders
		for score in ["LOSS"] #, "AUPRC", "TPR@5", "F1@5"] # LOSS->validation_likelihood
			@info "working on $(type)-$(encoder)-$(score)"
			models = models_and_params("/home/skvarvit/generativead/GenerativeAD.jl/data/experiments/$(type)/$(encoder)") #shortcut
			#models = models_and_params(datadir("experiments/$(type)/$(encoder)"))
			df = create_df(models, score=score, images=cheetsheet[type])
			@info "dataframe created -> size $(size(df))"
			# change name just because i would be easier to separater model name later
			#encoder = (encoder == "wae_vamp") ? "wae-vamp" : encoder
			CSV.write(datadir("tables/$((encoder == "wae_vamp") ? "wae-vamp" : encoder)_$(score)_$(type)_tab.csv"), df)
			CSV.write(datadir("tables/$((encoder == "wae_vamp") ? "wae-vamp" : encoder)_$(score)_$(type)_best_tab.csv"), return_best_n(df, score!="LOSS", 10))
			println("$(encoder)-$(score)-$(type) ... done")
		end
	end
end


text:
using StatsBase
function dstats(x, rnames="")
    k = size(x,2)
    if rnames==""
        rnames = 1:k
        rnames = rnames'
    end
    m = mean(x,1)
    s = std(x,1)
    sk = m-m
    k = m-m
    for i = 1:size(x,2)
        sk[:,i] = skewness(x[:,i])
        k[:,i] = kurtosis(x[:,i])
    end
    mn = minimum(x,1)
    mx = maximum(x,1)
    q05 = zeros(k)
    q95 = zeros(k)
    for i = 1:size(x,2) q05[i] = quantile(x[:,i], 0.05) end
    for i = 1:size(x,2) q95[i] = quantile(x[:,i], 0.95) end
    cnames = ["  mean", "  std", "skew", "kurt", "min", "max", "q05", "q95"]
    stats = [m' s' sk' k' mn' mx' q05' q95']
    prettyprint(stats, cnames, rnames);
    return stats
end


text:
#
# vect-tests.jl -
#
# Tests for vectorized operations.
#
module TestingLazyAlgebraVectorized

using LazyAlgebra
using Test

lpad_with_parentheses(n::Integer, x) = lpad(string("(", x, ")"), n)
lpad_with_parentheses(n::Integer, x, y) = lpad(string("(", x, ",", y, ")"), n)

distance(a::Real, b::Real) = abs(a - b)
distance(a::NTuple{2,Real}, b::NTuple{2,Real}) =
    hypot(a[1] - b[1], a[2] - b[2])
distance(A::AbstractArray{Ta,N}, B::AbstractArray{Tb,N}) where {Ta,Tb,N} =
    maximum(abs.(A - B))

makeselection(n::Integer) = begin
    sel = Array{Int}(undef, 0)
    j = [2,3,5]
    k = 1
    while k ≤ n
        push!(sel, k)
        k += rand(j)
    end
    return sel
end

alphas = (0, 1, -1,  2.71, π)
betas = (0, 1, -1, -1.33, Base.MathConstants.φ)
types = (Float32, Float64)
dims = (3,4,5)

@testset "vnorm $(lpad_with_parentheses(34,T))" for T in types
    S = (T == Float32 ? Float64 : Float32)
    v = randn(T, dims)
    @test vnorminf(v) == maximum(abs.(v))
    @test vnorminf(S, v) == S(maximum(abs.(v)))
    @test vnorm1(v) ≈ sum(abs.(v))
    @test vnorm2(v) ≈ sqrt(sum(v.*v))
    z = complex.(randn(T, dims), randn(T, dims))
    @test vnorminf(z) ≈ maximum(abs.(z))
    @test vnorm1(z) ≈ sum(abs.(real.(z)) + abs.(imag.(z)))
    @test vnorm2(z) ≈ sqrt(sum(abs2.(z)))
end # testset
@testset "vcopy, vswap $(lpad_with_parentheses(27,T))" for T in types
    u = randn(T, dims)
    uc = vcopy(u)
    @test distance(u, uc) == 0
    v = randn(T, dims)
    vc = vcopy!(vcreate(v), v)
    @test distance(v, vc) == 0
    vswap!(u, v)
    @test distance(u, vc) == distance(v, uc) == 0
    @test_throws DimensionMismatch  vcopy!(Array{T}(undef, dims .+ 1), u)
end # testset
@testset "vfill $(lpad_with_parentheses(34,T))" for T in types
    a = randn(T, dims)
    @test distance(vfill!(a,0), zeros(T,dims)) == 0
    a = randn(T, dims)
    @test distance(vfill!(a,0), vzero!(a)) == 0
    a = randn(T, dims)
    @test distance(vfill!(a,1), ones(T,dims)) == 0
    a = randn(T, dims)
    @test distance(vfill!(a,π), fill!(similar(a), π)) == 0
    ac = vcopy(a)
    @test distance(vzeros(a), zeros(T,dims)) == 0
    @test distance(a, ac) == 0
    @test distance(vones(a), ones(T,dims)) == 0
    @test distance(a, ac) == 0
end # testset
@testset "vscale                                  " begin
    for T in types
        a = randn(T, dims)
        b = vcreate(a)
        for α in (0, -1, 1, π, 2.71)
            d = T(α)*a
            @test distance(vscale(α,a), d) == 0
            @test distance(vscale(a,α), d) == 0
            @test distance(vscale!(b,α,a), d) == 0
            @test distance(vscale!(b,a,α), d) == 0
            c = vcopy(a)
            @test distance(vscale!(c,α), d) == 0
            vcopy!(c, a)
            @test distance(vscale!(α,c), d) == 0
        end
    end
    v = ones(dims)
    @test_throws ErrorException vscale!(1,2)
    @test_throws ErrorException vscale!(v,2,3)
    @test_throws ErrorException vscale!(1,v,3)
    @test_throws ErrorException vscale!(1,2,v)
    @test_throws ErrorException vscale!(1,2,3)
    for Ta in types, Tb in types
        a = randn(Ta, dims)
        ac = vcopy(a)
        b = Array{Tb}(undef, dims)
        e = max(eps(Ta), eps(Tb))
        for α in alphas
            d = α*a
            @test distance(vscale!(b,α,a), d) ≤ 8e
            @test distance(vscale(α,a), d) ≤ 8e
            @test distance(a, ac) == 0
        end
    end
end # testset
@testset "vupdate $(lpad_with_parentheses(32,Ta,Tb))" for Ta in types, Tb in types
    Tmin = sizeof(Ta) ≤ sizeof(Tb) ? Ta : Tb
    Tmax = sizeof(Ta) ≥ sizeof(Tb) ? Ta : Tb
    a = randn(Ta, dims)
    b = randn(Tb, dims)
    sel = makeselection(length(a))
    atol, rtol = zero(Tmin), sqrt(eps(Tmin))
    for α in alphas
        @test vupdate!(vcopy(a),α,b) ≈
            a + Tmax(α)*b atol=atol rtol=rtol norm=vnorm2
        c = vcopy(a)
        c[sel] .+= Tmax(α)*b[sel]
        @test vupdate!(vcopy(a),sel,α,b) ≈ c atol=atol rtol=rtol norm=vnorm2
    end
end # testset
@testset "vproduct $(lpad_with_parentheses(31,Ta,Tb))" for Ta in types, Tb in types
    a = randn(Ta, dims)
    b = randn(Tb, dims)
    sel = makeselection(length(a))
    c = vcreate(a)
    e = max(eps(Ta), eps(Tb))
    @test distance(vproduct!(c,a,b), (a .* b)) ≤ 2e
    c = randn(Tb, dims)
    d = vcopy(c)
    d[sel] = a[sel] .* b[sel]
    @test distance(vproduct!(c,sel,a,b), d) ≤ 2e
end # testset
@testset "vcombine $(lpad_with_parentheses(31,T))" for T in types
    a = randn(T, dims)
    b = randn(T, dims)
    d = vcreate(a)
    for α in alphas,
        β in betas
        @test distance(vcombine!(d,α,a,β,b), (T(α)*a + T(β)*b)) == 0
    end
end # testset
@testset "vdot $(lpad_with_parentheses(35,Ta,Tb))" for Ta in types, Tb in types
    Tmin = sizeof(Ta) ≤ sizeof(Tb) ? Ta : Tb
    Tmax = sizeof(Ta) ≥ sizeof(Tb) ? Ta : Tb
    a = randn(Ta, dims)
    b = randn(Tb, dims)
    sel = makeselection(length(a))
    w = zeros(Tmax, dims)
    w[sel] = rand(length(sel))
    # check implementation for real-valued vectors
    @test vdot(a,b) ≈ sum(a.*b)
    @test vdot(w,a,b) ≈ sum(w.*a.*b)
    @test vdot(sel,a,b) ≈ sum(a[sel].*b[sel])
    # check ⟨a,b⟩ = ⟨b,a⟩ for real-valued vectors
    @test vdot(a,b) == vdot(b,a)
    @test vdot(w,a,b) ≈ vdot(w,b,a)
    @test vdot(sel,a,b) == vdot(sel,b,a)
    @test sqrt(vdot(a,a)) ≈ vnorm2(a)
end # testset
@testset "vdot (Complex{$Ta},Complex{$Tb})" for Ta in types, Tb in types
    Tmin = sizeof(Ta) ≤ sizeof(Tb) ? Ta : Tb
    Tmax = sizeof(Ta) ≥ sizeof(Tb) ? Ta : Tb
    a = complex.(randn(Ta, dims), randn(Ta, dims))
    b = complex.(randn(Tb, dims), randn(Tb, dims))
    sel = makeselection(length(a))
    w = zeros(Tmax, dims)
    w[sel] = rand(length(sel))
    # check implementation for complex-valued vectors
    @test vdot(a,b) ≈ sum(conj.(a).*b)
    @test vdot(Tmax,a,b) ≈ real(sum(conj.(a).*b))
    @test vdot(w,a,b) ≈ sum(w.*conj.(a).*b)
    @test vdot(Tmax,w,a,b) ≈ real(sum(w.*conj.(a).*b))
    @test vdot(sel,a,b) ≈ sum(conj.(a[sel]).*b[sel])
    @test vdot(Tmax,sel,a,b) ≈ real(sum(conj.(a[sel]).*b[sel]))
    # check ⟨a,b⟩ = conj(⟨b,a⟩) for complex-valued vectors
    @test vdot(a,b) == conj(vdot(b,a))
    @test vdot(w,a,b) ≈ conj(vdot(w,b,a))
    @test vdot(sel,a,b) == conj(vdot(sel,b,a))
end # testset
nothing

end # module


text:
module JLLWrappers

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
    @eval Base.Experimental.@compiler_options compile=min optimize=0 infer=false
end

if VERSION >= v"1.6.0-DEV"
    using Preferences
end

# We need to glue expressions together a lot
function excat(exs::Union{Expr,Nothing}...)
    ex = Expr(:block)
    for exn in exs
        exn === nothing && continue
        if Meta.isexpr(exn, :block)
            append!(ex.args, exn.args)
        else
            push!(ex.args, exn)
        end
    end
    return esc(ex)
end

include("toplevel_generators.jl")
include("wrapper_generators.jl")
include("runtime.jl")

end # module

max_stars_repo_name: giordano/JLLWrappers

max_issues_repo_name: giordano/JLLWrappers

max_forks_repo_name: giordano/JLLWrappers


text:
function covmatrix(model::TimeSeriesModel, n, Δ)
    d = ndims(model)
    store = allocate_memory_EI_F(typeof(model), n , Δ)
    acv!(store, model)
    Areal = real.(extract_acv(store))
    C = Matrix{eltype(Areal)}(undef, n * d, n * d)
    for ii = 1:d # increase row
        for jj = 1:d # increase col
            # compute cross correlation
            C[(1:n).+(ii-1)*n, (1:n).+(jj-1)*n] =
                ii > jj ? # if in bottom triangle
                Toeplitz(
                    fftshift(Areal[indexLT(ii,jj,d), :])[n+1:end], # positive lags down
                    fftshift(Areal[indexLT(ii,jj,d), :])[n+1:-1:2], # negative lags across
                ) :
                ii < jj ? # if in top triangle
                Toeplitz(
                    fftshift(Areal[indexLT(ii,jj,d), :])[n+1:-1:2], # negative lags down
                    fftshift(Areal[indexLT(ii,jj,d), :])[n+1:end], # positive lags across
                ) : # if on diagonal
                Toeplitz(Areal[indexLT(ii,jj,d), 1:n], Areal[indexLT(ii,jj,d), 1:n]) # negative and positive lags equal
        end
    end
    if isposdef(C)
        return C
    else
        @warn "Covariance matrix not positive definite, adding a small number to diagonals to try and fix."
        for i ∈ 1:size(C,1)
            C[i,i] += 1e-10
        end
        if isposdef(C)
            return C
        else
            error("C not positive semi definite.")
        end
    end
end

function covmatrix(model::TimeSeriesModel{1}, n, Δ)
    store = allocate_memory_EI_F(typeof(model), n , Δ)
    acv!(store, model)
    Areal = real.(extract_acv(store))
    return Matrix(Toeplitz(Areal[1:n], Areal[1:n]))
end
struct GaussianProcess{T}
    X::T
    Δ::Float64
    function GaussianProcess(model::TimeSeriesModel, n::Int, Δ::Real)
        Δ > 0 || throw(ArgumentError("Δ should be positive."))
        n > 0 || throw(ArgumentError("n should be positive."))
        C = covmatrix(model, n, Δ)
        X = normaldist(model,C)
        new{typeof(X)}(X,Δ)
    end
end
normaldist(::TimeSeriesModel{D}, C) where {D} = reshape(MvNormal(C), size(C,1)÷D, D)
normaldist(::TimeSeriesModel{1}, C) = MvNormal(C)

struct TimeSeries{T,N}
    ts::Array{T,N}
    Δ::Float64
end
Base.ndims( ts::TimeSeries) = size(ts.ts, 2)
Base.length(ts::TimeSeries) = size(ts.ts, 1)

Distributions.rand(gp::GaussianProcess) = TimeSeries(rand(gp.X),gp.Δ)

@recipe function f(ts::TimeSeries)
    0:ts.Δ:ts.Δ*(length(ts)-1), ts.ts
end

function simulate_gp(model::TimeSeriesModel, n::Int, Δ::Number, nreps::Int)
    nreps > 0 || throw(ArgumentError("nreps should be positive."))
    GP = GaussianProcess(model, n , Δ)
    return [rand(GP) for _ in 1:nreps]
end


text:
using PrettyTables

function get_fermi()
    rex = r"@@@ Average Fock Time:\s([0-9]*\.?[0-9]*)"
    fpath = joinpath(@__DIR__, "fermi/output.dat")
    timings = zeros(22)
    i = 1
    for l = eachline(fpath)
        m = match(rex, l)
        if m !== nothing
            timings[i] = m.captures[1] |> String |> x->parse(Float64,x)
            i += 1
        end
    end

    @assert i == 23 # meaning that 22 entries were found

    return timings
end

function get_psi4()
    Frex = r"HF:\sForm\sF\s+:.+?([0-9]*\.?[0-9]*)w\s+([0-9]+)\scalls"
    Grex = r"HF:\sForm\sG\s+:.+?([0-9]*\.?[0-9]*)w\s+([0-9]+)\scalls"
    cd(joinpath(@__DIR__, "psi4"))

    timings = zeros(22)

    # For each molecule
    for i = eachindex(timings)
        Ncalls = []
        Ft = []
        Gt = []
        for l = eachline("S$i/timer.dat")
            Fm = match(Frex, l)
            Gm = match(Grex, l)

            if Fm !== nothing
                t = Fm.captures[1] |> String |> x->parse(Float64,x)
                N = Fm.captures[2] |> String |> x->parse(Int,x)
                push!(Ncalls, N)
                push!(Ft, t)
            end

            if Gm !== nothing
                t = Gm.captures[1] |> String |> x->parse(Float64,x)
                N = Gm.captures[2] |> String |> x->parse(Int,x)
                push!(Ncalls, N)
                push!(Gt, t)
            end
        end
        # Make sure all data is consistent
        @assert all(Ncalls .== Ncalls[1])
        @assert all(Ft .== Ft[1])
        @assert all(Gt .== Gt[1])

        # Compute timing
        timings[i] = (Ft[1] + Gt[1]) / Ncalls[1]
    end

    cd("..")
    return timings
end

if abspath(PROGRAM_FILE) == @__FILE__
    fermi = get_fermi()
    psi4 = get_psi4()

    @pt :header = ["Fermi", "Psi4"] [fermi psi4]
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using Catlab
using Catlab.WiringDiagrams
using Catlab.Doctrines
import Catlab.Doctrines: ⊗, id
import Base: ∘
include("SIR-simple__functions.jl")
⊗(a::WiringDiagram, b::WiringDiagram) = otimes(a, b)
∘(a::WiringDiagram, b::WiringDiagram) = compose(b, a)
⊚(a,b) = b ∘ a

SIR_simple__global__sir__0__beta__neg1, SIR_simple__global__sir__0__S__neg1, SIR_simple__global__sir__0__I__neg1, SIR_simple__global__sir__0__R__neg1, SIR_simple__global__sir__0__dt__neg1, SIR_simple__global__sir__0__gamma__neg1 = Ob(FreeSymmetricMonoidalCategory, :SIR_simple__global__sir__0__beta__neg1, :SIR_simple__global__sir__0__S__neg1, :SIR_simple__global__sir__0__I__neg1, :SIR_simple__global__sir__0__R__neg1, :SIR_simple__global__sir__0__dt__neg1, :SIR_simple__global__sir__0__gamma__neg1)
SIR_simple__global__sir__0__infected__0, SIR_simple__global__sir__0__recovered__0 = Ob(FreeSymmetricMonoidalCategory, :SIR_simple__global__sir__0__infected__0, :SIR_simple__global__sir__0__recovered__0)
SIR_simple__global__sir__0__I__0, SIR_simple__global__sir__0__S__0, SIR_simple__global__sir__0__R__0 = Ob(FreeSymmetricMonoidalCategory, :SIR_simple__global__sir__0__I__0, :SIR_simple__global__sir__0__S__0, :SIR_simple__global__sir__0__R__0)

id_SIR_simple__global__sir__0__R__neg1 = id(Ports([SIR_simple__global__sir__0__R__neg1]))
id_SIR_simple__global__sir__0__S__neg1 = id(Ports([SIR_simple__global__sir__0__S__neg1]))
id_SIR_simple__global__sir__0__I__neg1 = id(Ports([SIR_simple__global__sir__0__I__neg1]))

IN_0 = WiringDiagram(Hom(:L0_REWIRE, SIR_simple__global__sir__0__beta__neg1 ⊗ SIR_simple__global__sir__0__S__neg1 ⊗ SIR_simple__global__sir__0__I__neg1 ⊗ SIR_simple__global__sir__0__R__neg1 ⊗ SIR_simple__global__sir__0__dt__neg1 ⊗ SIR_simple__global__sir__0__gamma__neg1, SIR_simple__global__sir__0__beta__neg1 ⊗ SIR_simple__global__sir__0__S__neg1 ⊗ SIR_simple__global__sir__0__I__neg1 ⊗ SIR_simple__global__sir__0__R__neg1 ⊗ SIR_simple__global__sir__0__dt__neg1 ⊗ SIR_simple__global__sir__0__gamma__neg1 ⊗ SIR_simple__global__sir__0__I__neg1 ⊗ SIR_simple__global__sir__0__dt__neg1 ⊗ SIR_simple__global__sir__0__R__neg1 ⊗ SIR_simple__global__sir__0__S__neg1 ⊗ SIR_simple__global__sir__0__I__neg1))
WD_SIR_simple__global__sir__0__recovered__0 = WiringDiagram(Hom(SIR_simple__sir__assign__recovered__0, SIR_simple__global__sir__0__gamma__neg1 ⊗ SIR_simple__global__sir__0__I__neg1 ⊗ SIR_simple__global__sir__0__dt__neg1, SIR_simple__global__sir__0__recovered__0))
WD_SIR_simple__global__sir__0__infected__0 = WiringDiagram(Hom(SIR_simple__sir__assign__infected__0, SIR_simple__global__sir__0__beta__neg1 ⊗ SIR_simple__global__sir__0__S__neg1 ⊗ SIR_simple__global__sir__0__I__neg1 ⊗ SIR_simple__global__sir__0__R__neg1 ⊗ SIR_simple__global__sir__0__dt__neg1, SIR_simple__global__sir__0__infected__0))
OUT_1 = IN_0 ⊚ (WD_SIR_simple__global__sir__0__infected__0 ⊗ WD_SIR_simple__global__sir__0__recovered__0 ⊗ id_SIR_simple__global__sir__0__R__neg1 ⊗ id_SIR_simple__global__sir__0__S__neg1 ⊗ id_SIR_simple__global__sir__0__I__neg1)

IN_1 = WiringDiagram(Hom(:L1_REWIRE, SIR_simple__global__sir__0__infected__0 ⊗ SIR_simple__global__sir__0__recovered__0 ⊗ SIR_simple__global__sir__0__R__neg1 ⊗ SIR_simple__global__sir__0__S__neg1 ⊗ SIR_simple__global__sir__0__I__neg1, SIR_simple__global__sir__0__I__neg1 ⊗ SIR_simple__global__sir__0__infected__0 ⊗ SIR_simple__global__sir__0__recovered__0 ⊗ SIR_simple__global__sir__0__S__neg1 ⊗ SIR_simple__global__sir__0__infected__0 ⊗ SIR_simple__global__sir__0__R__neg1 ⊗ SIR_simple__global__sir__0__recovered__0))
WD_SIR_simple__global__sir__0__R__0 = WiringDiagram(Hom(SIR_simple__sir__assign__R__0, SIR_simple__global__sir__0__R__neg1 ⊗ SIR_simple__global__sir__0__recovered__0, SIR_simple__global__sir__0__R__0))
WD_SIR_simple__global__sir__0__S__0 = WiringDiagram(Hom(SIR_simple__sir__assign__S__0, SIR_simple__global__sir__0__S__neg1 ⊗ SIR_simple__global__sir__0__infected__0, SIR_simple__global__sir__0__S__0))
WD_SIR_simple__global__sir__0__I__0 = WiringDiagram(Hom(SIR_simple__sir__assign__I__0, SIR_simple__global__sir__0__I__neg1 ⊗ SIR_simple__global__sir__0__infected__0 ⊗ SIR_simple__global__sir__0__recovered__0, SIR_simple__global__sir__0__I__0))
OUT_2 = OUT_1 ⊚ IN_1 ⊚ (WD_SIR_simple__global__sir__0__I__0 ⊗ WD_SIR_simple__global__sir__0__S__0 ⊗ WD_SIR_simple__global__sir__0__R__0)


text:
# Define standard forward and reverse propagation to switch of expression definitions
# for expressions.
function binary_switch(ids; is_forward = true)
    if length(ids) <= 3
        if is_forward
            out = Expr(:if, Expr(:call, :(==), :id, ids[1]), :(return fprop!(t, Val($(ids[1])), g, c, k)))
        else
            out = Expr(:if, Expr(:call, :(==), :id, ids[1]), :(return rprop!(t, Val$((ids[1])), g, c, k)))
        end
        (length(ids) > 1) && push!(out.args, binary_switch(ids[2:end]))
        return out
    else
        mid = length(ids) >>> 1
        return Expr(:if, Expr(:call, :(<=), :id, ids[mid]), binary_switch(ids[1:mid]), binary_switch(ids[mid+1:end]))
    end
end

function Node(aux_info, d::JuMP._Derivatives.NodeData, child_vec, c::UnitRange{Int}, op)
    nt = d.nodetype
    i = d.index
    (nt == JuMP._Derivatives.CALL)          && return _create_call_node(i, child_vec, c, op)
    (nt == JuMP._Derivatives.CALLUNIVAR)    && return _create_call_node_uni(i, child_vec, c, op)
    (nt == JuMP._Derivatives.VALUE)         && return Node(Constant(), i)
    (nt == JuMP._Derivatives.PARAMETER)     && return Node(Parameter(), i)
    (nt == JuMP._Derivatives.SUBEXPRESSION) && return Node(Subexpression(), i)
    (nt == JuMP._Derivatives.VARIABLE)      && return !is_auxilliary_variable(aux_info, i) ? Node(Variable(), i) : Node(Select(), i)
    (nt == JuMP._Derivatives.LOGIC)         && error("Unable to load JuMP expression. Logical operators not currently supported.")
    (nt == JuMP._Derivatives.COMPARISON)    && error("Unable to load JuMP expression. Comparisons not currently supported.")
    error("Node type = $nt not expected from JuMP.")
end

function _convert_node_list(aux_info, x::Vector{JuMP._Derivatives.NodeData}, op)
    y = Vector{Node}(undef, length(x))
    adj = JuMP._Derivatives.adjmat(x)
    child_vec = rowvals(adj)
    for i in eachindex(x)
        y[i] = Node(aux_info, x[i], child_vec, nzrange(adj, i), op)
    end
    return y
end

# Access gradient sparsity of JuMP storage.
sparsity(d::JuMP._FunctionStorage) = d.grad_sparsity
sparsity(d::JuMP._SubexpressionStorage) = d.sparsity

# Compute gradient sparsity from JuMP storage.
function _compute_sparsity(d::JuMP._FunctionStorage, sparse_dict::Dict{Int,Vector{Int}}, is_sub, subexpr_indx)
    dep_subexpression = Int[]
    variable_dict = Dict{Int,Bool}()
    for n in d.nd
        if n.nodetype == JuMP._Derivatives.VARIABLE
            if !haskey(variable_dict, n.index)
                variable_dict[n.index] = true
            end
        end
        if n.nodetype == JuMP._Derivatives.SUBEXPRESSION
            push!(dep_subexpression, n.index)
        end
    end
    sparsity = collect(keys(variable_dict))
    unique!(dep_subexpression)
    sort!(dep_subexpression)
    for s in dep_subexpression
        append!(sparsity, sparse_dict[s])
    end
    unique!(sparsity)
    sort!(sparsity)
    if is_sub
        sparse_dict[subexpr_indx] = sparsity
    end
    sparsity, dep_subexpression
end
function _compute_sparsity(d::JuMP._SubexpressionStorage, sparse_dict::Dict{Int,Vector{Int}}, is_sub, subexpr_indx)
    dep_subexpression = Int[]
    variable_dict = Dict{Int,Bool}()
    for n in d.nd
        if n.nodetype == JuMP._Derivatives.VARIABLE
            if !haskey(variable_dict, n.index)
                variable_dict[n.index] = true
            end
        end
        if n.nodetype == JuMP._Derivatives.SUBEXPRESSION
            push!(dep_subexpression, n.index)
        end
    end
    sparsity = collect(keys(variable_dict))
    unique!(dep_subexpression)
    sort!(dep_subexpression)
    for s in dep_subexpression
        append!(sparsity, sparse_dict[s])
    end
    unique!(sparsity)
    sort!(sparsity)
    if is_sub
        sparse_dict[subexpr_indx] = sparsity
    end
    sparsity, dep_subexpression
end

function linearity(d::JuMP._Derivatives.Linearity)
    (d == JuMP._Derivatives.LINEAR)            && return LIN_LINEAR
    (d == JuMP._Derivatives.PIECEWISE_LINEAR)  && return LIN_PIECEWISE_LINEAR
    (d == JuMP._Derivatives.NONLINEAR)         && return LIN_NONLINEAR
    LIN_CONSTANT                               # assumes d is then JuMP._Derivatives.CONSTANT
end
function linearity(nd::Vector{JuMP._Derivatives.NodeData}, adj::SparseMatrixCSC{Bool,Int}, d::Vector{JuMP._Derivatives.Linearity})
    x = JuMP._Derivatives.classify_linearity(nd, adj, d)
    linearity.(x)
end

function OperatorRegistry(d::JuMP._Derivatives.UserOperatorRegistry)
    mv_id = collect(keys(d.multivariate_operator_to_id))
    mv_operator_to_id = d.multivariate_operator_to_id
    mv_operator_evaluator = d.multivariate_operator_evaluator
    u_operator_id = collect(keys(d.univariate_operator_to_id))
    u_to_id = d.univariate_operator_to_id
    u_to_f = d.univariate_operator_f
    u_fprime = d.univariate_operator_fprime
    u_fprimeprime = d.univariate_operator_fprimeprime
    OperatorRegistry(mv_id, mv_operator_to_id, mv_operator_evaluator, u_operator_id, u_to_id, u_to_f, u_fprime, u_fprimeprime)
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
@testset "prior.jl" begin
    @testset "types" begin
        @test DefaultPrior() isa T.Prior

        my_prior = CustomPrior(Normal(), Normal(), nothing)
        @test my_prior isa T.Prior

        my_prior = CustomPrior(Normal(), Normal(), Exponential(1))
        @test my_prior isa T.Prior
    end

    @testset "constructors" begin
        @test_throws MethodError DefaultPrior(Normal(0, 1))
        @test_throws ArgumentError CustomPrior()
    end

    @testset "CustomPrior" begin
        my_prior = CustomPrior(Normal(), Normal(), nothing)
        @test my_prior.predictors == Normal(0, 1)
        @test my_prior.intercept == Normal(0, 1)
        @test isnothing(my_prior.auxiliary)

        my_prior = CustomPrior(Normal(), Normal(), Exponential(1))
        @test my_prior.predictors == Normal(0, 1)
        @test my_prior.intercept == Normal(0, 1)
        @test my_prior.auxiliary == Exponential(1)
    end
end


text:
using TensorNetworkAD
using Test, Random
using Zygote
using LinearAlgebra: norm

@testset "autodiff" begin
    a = randn(10,10)
    @test Zygote.gradient(norm, a)[1] ≈ num_grad(norm, a)

    foo = x -> sum(Float64[x x; x x])
    @test Zygote.gradient(foo, 1)[1] ≈ num_grad(foo, 1)
end


text:
module MMS
"""
Module for working with the method of manufactured solutions.

"""

using DataStructures

"""
Extract slices from an Array of ordered lengths n1, n2, n3. The sum of lengths
cannot exceed the array length. The dictionary `vars` is used to label each slice.

Example:
        ex = extract_vars([1,2,3,4,5], OrderedDict("a" => 2, "b" => 3))
        #ex["a"] = [1,2]
        #ex["b"] = [3,4,5]

"""
function extract_vars(solution::Array, vars::OrderedDict)
        extracted_vars = OrderedDict()

        offset = 0
        for (key, value) in vars
                extracted_vars[key] = solution[1+offset:offset+value]
                offset += value
        end

        return extracted_vars
end

"""
Return the l2-error for each variable stored in a dict.
The l2-error is computed as \$\\|e\\| = \\sqrt{e^T H e}\$,
where e is the error vector, and H is a norm matrix.
"""
function l2_errors(solution::Dict, mms::Dict, norms::Dict)
        err = Dict()
        for (var, value) in solution
                e = (solution[var] .- mms[var])
                err[var] = sqrt(e' * norms[var] * e)
        end
        return err
end

function abs_errors(solution::Dict, mms::Dict)
        err = Dict()
        for (var, value) in solution
                err[var] = abs.(solution[var] - mms[var])
        end
        return err
end

"""
Return the functional error for each variable stored in a dict.  Consider the
linear functional \$ F(u) =  \\int u dx \$. The (absolute) error is defined as \$e = |F(u) -
\\sum_j (H * v)_j| \$, where H is a quadrature rule and v is the numerical
solution.
"""
function functional_errors(solution::Dict, mms::Dict, norms::Dict)
        err = Dict()
        for (var, value) in solution
                err[var] = abs(mms[var] - sum(norms[var] * solution[var]))
        end
        return err
end

"""
Compute the convergence rate for a list of errors that are computed with
subsequent grid refinement of a factor of two.

"""
function log2_convergence_rates(err::Array)

    q = -log2.(err[2:end]./err[1:end-1])
    prepend!(q, Inf)
    return q
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
module TestData

# using Revise
using Test
using DataFrames
import TypedTables
# using StatsBase
# using JuliaDB
# using SparseArrays
using CategoricalArrays
import Tables
using ScientificTypes

using Random
import Random.seed!
seed!(1234)

import MLJBase: decoder, int, classes, partition, unpack, selectcols, matrix,
    CategoricalElement, selectrows, select, table, nrows

@testset "partition" begin
    train, test = partition(1:100, 0.9)
    @test collect(train) == collect(1:90)
    @test collect(test) == collect(91:100)
    train, test = partition(1:100, 0.9, shuffle=true)
    @test length(train) == 90

    train, test = partition(1:100, 0.9, shuffle=true, rng=1)
    @test length(train) == 90

    train, test = partition(1:100, 0.9, shuffle=true,
                            rng=Random.MersenneTwister(3))
    @test length(train) == 90
end

@testset "unpack" begin
    channing = TypedTables.Table(
        Sex = categorical(["Female", "Male", "Female"]),
        Entry = Int32[965, 988, 850],
        Exit = Int32[1088, 1045, 940],
        Time = Int32[123, 57, 90],
        Cens = Int32[0, 0, 1],
        weights = [1,2,5])

    w, y, X =  unpack(channing,
                  ==(:weights),
                  ==(:Exit),
                  x -> x != :Time;
                  :Exit=>Continuous,
                  :Entry=>Continuous,
                  :Cens=>Multiclass)

    @test w == selectcols(channing, :weights)
    @test y == selectcols(channing, :Exit)
    @test X == selectcols(channing, [:Sex, :Entry, :Cens])
    @test scitype_union(y) <: Continuous
    @test scitype_union(selectcols(X, :Cens)) <: Multiclass

    w, y, X =  unpack(channing,
                      ==(:weights),
                      ==(:Exit),
                      x -> x != :Time;
                      wrap_singles=true,
                      :Exit=>Continuous,
                      :Entry=>Continuous,
                      :Cens=>Multiclass)

    @test selectcols(w, 1)  == selectcols(channing, :weights)
    @test selectcols(y, 1)  == selectcols(channing, :Exit)
    @test X == selectcols(channing, [:Sex, :Entry, :Cens])

    @test_throws(Exception, unpack(channing,
                                   ==(:weights),
                                   ==(:Exit),
                                   ==(:weights),
                                   x -> x != :Time;
                                   :Exit=>Continuous,
                                   :Entry=>Continuous,
                                   :Cens=>Multiclass))

end

@testset "categorical element decoder, classes " begin

    N = 10
    mix = shuffle(0:N - 1)

    Xraw = broadcast(x->mod(x,N), rand(Int, 2N, 3N))
    Yraw = string.(Xraw)

    # to turn a categ matrix into a ordinary array with categorical
    # elements. Need because broacasting the identity gives a
    # categorical array in CategoricalArrays >0.5.2
    function matrix_(X)
        ret = Array{Any}(undef, size(X))
        for i in eachindex(X)
            ret[i] = X[i]
        end
        return ret
    end

    X = categorical(Xraw)
    x = X[1]
    Y = categorical(Yraw)
    y = Y[1]
    V = matrix_(X)
    W = matrix_(Y)
    Xo = levels!(deepcopy(X), mix)
    xo = Xo[1]
    Yo = levels!(deepcopy(Y), string.(mix))
    yo = Yo[1]
    Vo = matrix_(Xo)
    Wo = matrix_(Yo)

    raw(x::CategoricalElement) = x.pool.index[x.level]
    @test raw.(classes(xo)) == xo.pool.levels
    @test raw.(classes(yo)) == yo.pool.levels

    # getting all possible elements from one:
    @test raw.(X) == Xraw
    @test raw.(Y) == Yraw
    @test raw.(classes(xo)) == levels(Xo)
    @test raw.(classes(yo)) == levels(Yo)

    # broadcasted encoding:
    @test int(X) == int(V)
    @test int(Y) == int(W)

    # encoding is right-inverse to decoding:
    d = decoder(xo)
    @test d(int(Vo)) == Vo # ie have the same elements
    e = decoder(yo)
    @test e(int(Wo)) == Wo

    # classes return value is in correct order:
    int(classes(xo)) == 1:length(classes(xo))

#    levels!(Yo, reverse(levels(Yo)))
#    e = decoder(yo)
#    @test e(int(Wo)) == Wo

    # int is based on ordering not index
    v = categorical(['a', 'b', 'c'], ordered=true)
    @test int(v) == 1:3
    levels!(v, ['c', 'a', 'b'])
    @test int(v) == [2, 3, 1]

    # special private int(pool, level) method:
    @test broadcast(int, [d.pool], levels(Xo)) == int.(classes(xo))
    levels!(Xo, reverse(levels(Xo)))
    @test broadcast(int, [d.pool], levels(Xo)) == int.(classes(xo))

end

@testset "table to matrix to table" begin

    B = rand(UInt8, (4, 5))
    @test matrix(DataFrame(B)) == B
    @test matrix(table(B)) == B
    @test matrix(table(B), transpose=true) == B'

    X  = (x1=rand(5), x2=rand(5))

    @test table(X, prototype=DataFrame()) == DataFrame(X)
    T = table((x1=(1,2,3), x2=(:x, :y, :z)))
    selectcols(T, :x1) == [1, 2, 3]

    v = categorical(11:20)
    A = hcat(v, v)
    tab = table(A)
    selectcols(tab, 1) == v

    @test matrix(B) == B
    @test matrix(B, transpose=true) == permutedims(B)

end

## TABLE INDEXING

A = broadcast(x->Char(65+mod(x,5)), rand(Int, 10, 5))
X = CategoricalArrays.categorical(A)

df = DataFrame(A)
tt = TypedTables.Table(df)
# uncomment 4 lines to restore JuliaDB testing:
# db = JuliaDB.table(tt)
# nd = ndsparse((document=[6, 1, 1, 2, 3, 4],
#                word=[:house, :house, :sofa, :sofa, :chair, :house]), (values=["big", "small", 17, 34, 4, "small"],))

df.z  = 1:10
tt = TypedTables.Table(df)
# uncomment 1 line to restore JuliaDB testing:
# db = JuliaDB.table(tt)

@test selectcols(df, 4:6) == df[:,4:6]
@test selectcols(df, [:x1, :z]) == df[:,[:x1, :z]]
@test selectcols(df, :x2) == df.x2
@test selectcols(df, 2) == df.x2
@test selectrows(df, 4:6) == selectrows(df[4:6, :], :)
@test selectrows(df, 1) == selectrows(df[1:1, :], :)
@test select(df, 2, :x2) == df[2,:x2]
s = schema(df)
@test nrows(df) == size(df, 1)

@test selectcols(tt, 4:6) == selectcols(TypedTables.Table(x4=tt.x4, x5=tt.x5, z=tt.z), :)
@test selectcols(tt, [:x1, :z]) == selectcols(TypedTables.Table(x1=tt.x1, z=tt.z), :)
@test selectcols(tt, :x2) == tt.x2
@test selectcols(tt, 2) == tt.x2
@test selectrows(tt, 4:6) == selectrows(tt[4:6], :)
@test nrows(tt) == length(tt.x1)
@test select(tt, 2, :x2) == tt.x2[2]

@testset "vector accessors" begin
    v = rand(Int, 4)
    @test selectrows(v, 2:3) == v[2:3]
    @test selectrows(v, 2) == [v[2]]
    @test nrows(v) == 4

    v = categorical(collect("asdfasdf"))
    @test selectrows(v, 2:3) == v[2:3]
    @test selectrows(v, 2) == [v[2]]
    @test nrows(v) == 8
end

@testset "matrix accessors" begin
    A = rand(5, 10)
    @test selectrows(A, 2:4) == A[2:4,:]
    @test selectrows(A, 2:4) == A[2:4,:]
    @test selectrows(A, 2) == A[2:2,:]

    A = rand(5, 10) |> categorical
    @test selectrows(A, 2:4) == A[2:4,:]
    @test selectrows(A, 2:4) == A[2:4,:]
    @test selectrows(A, 2) == A[2:2,:]

    @test nrows(A) == 5
end

v = categorical(collect("asdfasdf"))
df = DataFrame(v=v, w=v)
@test selectcols(df, :w) == v
tt = TypedTables.Table(df)
# uncomment 1 line to restore JuliaDB testing:
# db = JuliaDB.table(tt)
@test selectcols(tt, :w) == v
# uncomment 1 line to restore JuliaDB testing:
# @test selectcols(db, :w) == v

# uncomment 9 lines to restore JuliaDB testing:
# @test select(nd, 2, :house) isa Missing
# @test select(nd, 1, :house) == "small"
# @test all(select(nd, :, :house) .=== ["small", missing, missing, "small", missing, "big"])
# @test all(select(nd, [2,4], :house) .=== [missing, "small"])
# @test all(selectcols(nd, :house) .=== select(nd, :, :house))
# @test nrows(nd) == 6
# s = schema(nd)
# @test s.names == (:chair, :house, :sofa)
# @test s.types == (Union{Missing,Int64}, Union{Missing,String}, Union{Missing,Int64})

# uncomment 3 lines to restore JuliaDB testing:
# sparsearray = sparse([6, 1, 1, 2, 3, 4], [2, 2, 3, 3, 1, 2],
#                      ["big", "small", 17, 34, 4, "small"])
# @test matrix(nd) == sparsearray

end # module

true


text:
int_rules_8_10 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 8.10*Bessel*functions =#
    @apply_utils Antiderivative(BesselJ(1, ~(a') + ~(b') * ~x), ~x) => -(BesselJ(0, ~a + ~b * ~x)) / ~b <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(BesselJ(~n, ~(a') + ~(b') * ~x), ~x) => (-2 * BesselJ(~n - 1, ~a + ~b * ~x)) / ~b + Antiderivative(BesselJ(~n - 2, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x) && IGtQ((~n - 1) / 2, 0)
    #= Antiderivative(BesselJ((~n),(~a')+(~b')*(~x)),~x) :=  (-1)^n*Antiderivative(BesselJ(-n,a+b*x),x) <-- FreeQ([a,b],x) && ILtQ(n,0) =#
    @apply_utils Antiderivative(BesselJ(~n, ~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) ^ (~n + 1) * HypergeometricPFQ([(~n + 1) / 2], [(~n + 3) / 2, ~n + 1], (-1 / 4) * (~a + ~b * ~x) ^ 2)) / (2 ^ ~n * ~b * Gamma(~n + 2)) <-- FreeQ([~a, ~b, ~n], ~x)
end

max_stars_repo_path: src/rules/8 Special functions/8.10 Bessel functions.jl


text:
using Test
using SomervilleCouncilParser
using SomervilleCouncilParser: validate_date, request_agenda_items, agenda_version
using DataFrames, Dates
using Legolas

const TEST_ASSETS = joinpath(pkgdir(SomervilleCouncilParser), "test", "assets")

include("meeting.jl")
include("agenda.jl")

# If any of the following examples needs adjusting, be sure to update the
# README.md examples! Should probably be handled automatically by,
# e.g., doctests---but until then, this'll have to suffice.
@testset "Project README.md" begin
    meetings = request_meetings("9/1/2021", "9/30/2021")
    @test nrow(meetings) == 13
    @test length(names(meetings)) == 4
    @test length(unique(meetings.name)) == 10

    results = search_agendas_for_content("9/1/2021", "9/30/2021", ["fluff"])
    @test nrow(results.meetings) == 13
    @test length(names(results.meetings)) == 7
    @test nrow(results.items) == 1
    @test isnothing(display_items_by_meeting(results.items))
end


text:
include("helper.jl")
include("minitest.jl")
Minitest.do_tests() do
  #Minitest.include("test_confusion_matrix.jl")
  #Minitest.include("test_entropy.jl")
  Minitest.include("test_compress.jl")
  Minitest.include("test_ncd.jl")
  Minitest.include("test_ncd_matrix.jl")
end


text:
using Documenter, Bobby

makedocs(
    format = :html,
    sitename = "Bobby.jl",
)


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

using DiscreteEvents, Random, Unitful, Test, .Threads, DataStructures
import Unitful: Time, ms, s, minute, hr

x = 2 # set global (Main) variable for mocking fclosure.jl doctest line 90

println(".... testing DiscreteEvents.jl .....")
@testset "clock.jl" begin
    include("test_clock.jl")
end

@testset "threads.jl" begin
    include("test_threads.jl")
end

@testset "channel 1 example" begin
    include("test_channels1.jl")
end

@testset "process.jl" begin
    include("test_process.jl")
end

if (VERSION ≥ v"1.3") && (nthreads() > 1)
    @testset "timer.jl" begin
        include("test_timer.jl")
    end
else
    println("... timer requires Julia ≥ 1.3, nthreads() > 1")
end

@testset "resources.jl" begin
    include("test_resources.jl")
end

# if (VERSION ≥ v"1.3") && (nthreads() > 1)
#     @testset "doctests" begin
#         include("test_docs.jl")
#     end
# else
#     println("... no doctests with Julia $VERSION, nthreads: $(nthreads()) ...")
# end
#
@testset "examples" begin
    include("test_examples.jl")
end

println(".... finished testing DiscreteEvents.jl ....")


text:
include("single_room_undirected.jl")
include("single_room_directed.jl")
include("grid_rooms_undirected.jl")
include("grid_rooms_directed.jl")
include("sequential_rooms_undirected.jl")
include("sequential_rooms_directed.jl")
include("maze_undirected.jl")
include("maze_directed.jl")
include("go_to_target_undirected.jl")
include("go_to_target_directed.jl")
include("door_key_undirected.jl")
include("door_key_directed.jl")
include("collect_gems_undirected.jl")
include("collect_gems_directed.jl")
include("collect_gems_multi_agent_undirected.jl")
include("dynamic_obstacles_undirected.jl")
include("dynamic_obstacles_directed.jl")
include("sokoban/sokoban_undirected.jl")
include("sokoban/sokoban_directed.jl")
include("snake.jl")
include("catcher.jl")
include("transport_undirected.jl")
include("transport_directed.jl")

const ENVS = [
              SingleRoomUndirectedModule.SingleRoomUndirected,
              SingleRoomDirectedModule.SingleRoomDirected,
              GridRoomsUndirectedModule.GridRoomsUndirected,
              GridRoomsDirectedModule.GridRoomsDirected,
              SequentialRoomsUndirectedModule.SequentialRoomsUndirected,
              SequentialRoomsDirectedModule.SequentialRoomsDirected,
              MazeUndirectedModule.MazeUndirected,
              MazeDirectedModule.MazeDirected,
              GoToTargetUndirectedModule.GoToTargetUndirected,
              GoToTargetDirectedModule.GoToTargetDirected,
              DoorKeyUndirectedModule.DoorKeyUndirected,
              DoorKeyDirectedModule.DoorKeyDirected,
              CollectGemsUndirectedModule.CollectGemsUndirected,
              CollectGemsDirectedModule.CollectGemsDirected,
              CollectGemsMultiAgentUndirectedModule.CollectGemsMultiAgentUndirected,
              DynamicObstaclesUndirectedModule.DynamicObstaclesUndirected,
              DynamicObstaclesDirectedModule.DynamicObstaclesDirected,
              SokobanUndirectedModule.SokobanUndirected,
              SokobanDirectedModule.SokobanDirected,
              SnakeModule.Snake,
              CatcherModule.Catcher,
              TransportUndirectedModule.TransportUndirected,
              TransportDirectedModule.TransportDirected,
             ]


text:
# adjoint of diff/cumsum operators for judiVector

export adjoint_diff, adjoint_cumsum

function adjoint_diff(d_obs::judiVector)
    d_out = deepcopy(d_obs)
    for i = 1:d_obs.nsrc
        d_out.data[i][1:end-1,:] = d_obs.data[i][1:end-1,:] - d_obs.data[i][2:end,:]
    end
    return 1f0/d_obs.geometry.dt[1]*d_out
end

function adjoint_cumsum(d_obs::judiVector;dims=1)
    d_out = deepcopy(d_obs)
    for i = 1:d_obs.nsrc
        d_out.data[i] = reverse(cumsum(reverse(d_obs.data[i],dims=dims),dims=dims),dims=dims)
    end
    return d_obs.geometry.dt[1]*d_out
end


text:
using RetroCap
using Test

@testset "recompress function" begin
    RetroCap.with_temp_dir() do tmp_dir
        cd(tmp_dir)
        run(`git clone https://github.com/JuliaRegistries/General.git`)
        RetroCap.Compress.recompress(joinpath(tmp_dir, "General"))
    end
end


text:
# Load required packages
using Random, Statistics, LinearAlgebra, DataFrames, DataFramesMeta, CSV, FixedEffectModels, MixedModels

# Tell Julia where to-be-called functions are
include("corrcov.jl")
include("estimatelearning.jl")
include("initlearning.jl")
include("estimatelearningcore.jl")

# read in data that was created in Stata
df = CSV.read("nlswlearn.csv")
ID = df.idcode
N = 4711
T = 15
S = 1
y = df.ln_wage
x = [ones(N*T) df.exper df.exper.^2 df.collgrad df.race.==1]
Choice = df.Choice

# Estimate learning model
bstart,sig,Delta,DeltaCorr = estimatelearning(ID,N,T,S,y,x,Choice);
#ID is a N*T x 1 vector,
#N is an integer,
#T is an integer,
#S should be set to 1 for now
#y is a N*T x J matrix of outcomes (with 999 where outcomes are unobserved),
#x is a N*T x K x J covariate tensor
#Choice is a N*T x 1 vector of integers (1, ..., J+1) where J+1 refers to "not receive a signal"

# print estimates
println("sqrt Delta = ",sqrt(Delta))
println("sqrt sig = ",sqrt(sig))
println("betahat = ",bstart)

# compare with FE/RE
dfuse = df[df.ln_wage.!=999,:]
# FE
@show reg(dfuse, @formula(ln_wage ~ 1 + exper*exper + collgrad + race1 + fe(idcode)), Vcov.cluster(:idcode))
# RE
categorical!(dfuse, :idcode)
fm1 = fit(MixedModel, @formula(ln_wage ~ 1 + exper*exper + collgrad + race1 + (1|idcode)), dfuse)

# recover beliefs
sig_eps = .092187
df = @transform(df, signal = :ln_wage .- coef(fm1)[1] .- coef(fm1)[2].*:exper .- coef(fm1)[3]*:collgrad .- coef(fm1)[4].*:race1 .- coef(fm1)[5].*:exper.^2,
                    priorEbelief = zeros(length(:ln_wage)),
                    postrEbelief = zeros(length(:ln_wage)),
                    priorVbelief = 0.106297*ones(length(:ln_wage)),
                    postrVbelief = 0.106297*ones(length(:ln_wage)))

for i = 1:N
    for t=1:T
        rowt = (i-1)*T+t
        row1 = (i-1)*T+t+1
        if df.ln_wage[rowt]==999
            df.signal[rowt] = 0
            df.postrEbelief[rowt] = df.priorEbelief[rowt]
            df.postrVbelief[rowt] = df.priorVbelief[rowt]
        else
            df.postrEbelief[rowt] = df.priorEbelief[rowt]*(sig_eps./(sig_eps + df.priorVbelief[rowt])) + df.signal[rowt]*(df.priorVbelief[rowt])./(sig_eps + df.priorVbelief[rowt])
            df.postrVbelief[rowt] = df.priorVbelief[rowt]*(sig_eps./(sig_eps + df.priorVbelief[rowt]))
        end
        if t<T
            df.priorEbelief[row1] = df.postrEbelief[rowt]
            df.priorVbelief[row1] = df.postrVbelief[rowt]
        end
    end
end
#


text:
# various types of boundary conditions
abstract type ConvBoundary end
# no adjustments made
struct Periodic <: ConvBoundary end
# padded with zeros
struct Pad{N} <: ConvBoundary
    padBy::NTuple{N,Int}
end
"""
    Pad(x::Vararg{<:Integer, N}) where N = Pad{N}(x)

N gives the number of dimensions of convolution, while `x` gives the specific amount to pad in each dimension (done on both sides). If the values in `x` are negative, then the support of the filters will be determined automataically
"""
Pad(x::Vararg{<:Integer,N}) where N = Pad{N}(x)
import Base.ndims
ndims(p::Pad{N}) where N = N

function Base.show(io::IO, p::Pad{N}) where N
    print(io, "Pad$(p.padBy)")
end

# as in the DCT2 TODO: implement this in a way that is as space efficient as a DCT2, instead of doubling things
struct Sym <: ConvBoundary end

"""
    originalSize(sz, boundary::ConvBoundary)
Given the size `sz` from the fft, return the size of the input.
"""
function originalSize(sz, boundary::Pad{N}) where N
    return sz .- 2 .* boundary.padBy# ([sz[ii]-2*p for p in boundary.padBy]..., )
end

function originalSize(sz, boundary::Periodic) where N
    return sz# ([sz[1] for p in boundary.padBy]..., )
end

# TODO: when Sym gets implemented using a CFT, this should shrink
function originalSize(sz, boundary::Sym) where N
    return Int.(sz ./ 2) # ([Int(sz[ii]/2) for p in boundary.padBy]..., )
end

"""
    effectiveSize(sz, boundary::ConvBoundary)
Given the input size, figure out the resulting size
"""
function effectiveSize(sz, boundary::Pad{N}) where N
    if boundary.padBy[1] > 0
        return sz .+ 2 .* boundary.padBy, boundary
    elseif boundary.padBy[1] == 0
        return sz, Periodic
    elseif boundary.padBy[1] < 0
        newBoundary = Pad(([sz[i] >> 2 + 1 for i = 1:N]...,))
        return effectiveSize(sz, newBoundary)[1], newBoundary
    end
end

function effectiveSize(sz, boundary::Periodic) where N
    return sz, boundary
end

# TODO: when Sym gets implemented using a CFT, this should shrink
function effectiveSize(sz, boundary::Sym) where N
    return 2 .* sz, boundary
end

# apply the boundary condition appropriately to the input matrix. nd is the
# number of dimensions that we're working with
function applyBC(x, bc::Pad, nd)
    usedInds = ([bc.padBy[ii] .+ (1:size(x, ii)) for ii = 1:ndims(bc)]...,)
    return (pad(x, bc.padBy), usedInds)
end

function applyBC(x, bc::Periodic, nd)
    return (x, axes(x)[1:nd])
end

function applyBC(x, bc::Sym, nd)
    flipThisDim = cat(x, reverse(x, dims=nd), dims=nd)
    if nd == 1
        return flipThisDim, axes(x)[1:nd]
    else
        return applyBC(flipThisDim, bc, nd - 1)[1], axes(x)[1:nd]
    end
end

# padding methods
for (TYPE, CONVERT) in ((AbstractArray, identity),
    (CuArray, cu))
    @eval begin
        # 1D padding
        function pad(x::$TYPE{T,N}, padBy::Union{<:Integer,NTuple{1,<:Integer}}) where {T,N}
            szx = size(x)
            padded = cat($CONVERT(zeros(T, padBy[1], szx[2:end]...)),
                             x,
                             $CONVERT(zeros(T, padBy[1], szx[2:end]...)), dims=(1,))
            return padded
        end

        # 2D padding
        function pad(x::$TYPE{T,N}, padBy::NTuple{2,<:Integer}) where {T,N}
            szx = size(x)
            firstRow = cat($CONVERT(zeros(T, padBy...,  szx[3:end]...)),
                           $CONVERT(zeros(T, padBy[1], szx[2:end]...)),
                           $CONVERT(zeros(T, padBy...,  szx[3:end]...)), dims=2)
            secondRow = cat($CONVERT(zeros(T, szx[1], padBy[2], szx[3:end]...)),
                            x,
                            $CONVERT(zeros(T, szx[1], padBy[2], szx[3:end]...)),dims=(2,))
            thirdRow = cat($CONVERT(zeros(T, padBy...,  szx[3:end]...)),
                           $CONVERT(zeros(T, padBy[1], szx[2:end]...)),
                           $CONVERT(zeros(T, padBy...,  szx[3:end]...)), dims=(2,))
            return cat(firstRow, secondRow, thirdRow, dims=(1,))
        end
    end
end

Zygote.@adjoint function pad(x, padBy::Union{<:Integer,NTuple{1,<:Integer}})
    return pad(x, padBy), function (Δ)
        axΔ = axes(Δ)
        return (Δ[(1 + padBy[1]):(end - padBy[1]), axΔ[2:end]...],
                nothing)
    end
end

Zygote.@adjoint function applyBC(x, bc::Periodic, nd)
    return applyBC(x, bc, nd), function (Δ)
        xbc, usedInds = Δ
        return xbc, nothing, nothing
    end
end

Zygote.@adjoint function applyBC(x, bc::Sym, nd)
    return applyBC(x, bc, nd), function (Δ)
        xbc, usedInds = Δ
        ax = axes(x)
        aΔ = axes(xbc)
        return xbc[ax[1:nd]..., aΔ[nd + 1:end]...], nothing, nothing
    end
end


text:
using ApproximateDPs
using Test

@testset "ApproximateDPs.jl" begin
    # Write your tests here.
end


text:
"""
A point in the problem's search space with the
corresponding fitness value.

`F` is the original problem's fitness type
"""
abstract type FitIndividual{F} end

fitness_type(indi::FitIndividual{F}) where {F} = F

"""
Get the problem parameters (a point in the search space) of the individual.
"""
params(indi::FitIndividual) = indi.params

"""
    fitness(indi::FitIndividual)

Gets the fitness of the individual.
"""
fitness(indi::FitIndividual) = indi.fitness

"""
Candidate solution to the problem.

Candidate can be either a member of the population (`index` > 0) or
a standalone solution (`index` == -1).
Can carry additional information, like the `tag` or the genetic operator applied (`extra`).
"""
mutable struct Candidate{F} <: FitIndividual{F}
    params::Individual
    index::Int          # index of individual in the population, -1 if unassigned # FIXME nothing if unassigned
    fitness::F          # fitness

    extra::Any          # extra information
    tag::Int            # additional information set by the genetic operator

    Candidate{F}(params::Individual, index::Int = -1,
                 fitness::F = NaN,
                 extra::Any = nothing,
                 tag::Int = 0) where {F} =
        new(params, index, fitness, extra, tag)

    Candidate(params::Individual, index::Int = -1,
              fitness::F = NaN,
              extra::Any = nothing,
              tag::Int = 0) where {F} =
        new{F}(params, index, fitness, extra, tag)
end

index(cand::Candidate) = cand.index
tag(cand::Candidate) = cand.tag

Base.copy(c::Candidate) = Candidate(copy(c.params), c.index, c.fitness, c.extra, c.tag)

function Base.copy!(c::Candidate{F}, o::Candidate{F}) where F
    copyto!(c.params, o.params)
    c.index = o.index
    c.fitness = o.fitness # FIXME if vector?
    c.extra = o.extra
    c.tag = o.tag
    return c
end