text:
#######################################################################################################################################################################################################
#
# Changes to the struct
# General
#     2022-Jan-14: add van der Tol model struct
#     2022-Jan-24: fix documentation
#     2022-Feb-07: remove the hierachy from abstract fluorescence model
#     2022-Feb-07: move struct definition as a field of VJPReactionCenter
# Sources
#     van der Tol et al. (2014) Models of fluorescence and photosynthesis for interpreting measurements of solar-induced chlorophyll fluorescence
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Structure that stores van der Tol et al. (2014) fluorescence model parameters.

# Fields

$(TYPEDFIELDS)

"""
struct VanDerTolFluorescenceModel{FT<:AbstractFloat}
    "Fitting parameter K_0"
    K_0::FT
    "Fitting parameter α"
    K_A::FT
    "Fitting parameter β"
    K_B::FT
end

#######################################################################################################################################################################################################
#
# Changes to the constructor
# General
#     2022-Jan-14: migrate from Photosynthesis.jl
#     2022-Feb-24: set default to true to make the response more reasonable
# To do
#     TODO: examine why van der Tol et al has the nonstressed parameter set that are so far away from the stressed one
# Sources
#     van der Tol et al. (2013) Models of fluorescence and photosynthesis for interpreting measurements of solar-induced chlorophyll fluorescence
#
#######################################################################################################################################################################################################
"""

    VanDerTolFluorescenceModel{FT}(drought::Bool = true) where {FT<:AbstractFloat}

Constructor for `VanDerTolFluorescenceModel` fluorescence model, given
- `drought` If true, use parameters trained from drought stressed plant. Default is `true`.

---
# Examples
```julia
vdt = VanDerTolFluorescenceModel{Float64}();
vdt = VanDerTolFluorescenceModel{Float64}(false);
```
"""
VanDerTolFluorescenceModel{FT}(drought::Bool = true) where {FT<:AbstractFloat} = (
    if drought
        return VanDerTolFluorescenceModel{FT}(5.01, 1.93, 10)
    else
        return VanDerTolFluorescenceModel{FT}(2.48, 2.83, 0.114)
    end
);

#######################################################################################################################################################################################################
#
# Changes to this type
# General
#     2022-Jan-24: abstractize the reaction center
#     2022-Jan-25: fix documentation
#     2022-Feb-07: fix documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Structure that stores reaction center information

Hierachy of the `AbstractReactionCenter`
- [`VJPReactionCenter`](@ref)
- [`CytochromeReactionCenter`](@ref)
"""
abstract type AbstractReactionCenter{FT<:AbstractFloat} end

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Jan-15: isolate the reaction center from Leaf in Photosynthesis.jl
#     2022-Jan-25: fix documentation
#     2022-Feb-07: add fluorescence model as a field
#     2022-Mar-01: fix documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Structure that stores reaction center information

# Fields

$(TYPEDFIELDS)

"""
mutable struct VJPReactionCenter{FT<:AbstractFloat} <:AbstractReactionCenter{FT}
    # parameters that do not change with time
    "Fraction of absorbed light used by PSII ETR"
    F_PSII::FT
    "Fluorescence model"
    FLM::VanDerTolFluorescenceModel{FT}
    "Rate constant for thermal dissipation"
    K_D::FT
    "Rate constant for fluorescence"
    K_F::FT
    "Maximal rate constant for photochemistry"
    K_P_MAX::FT
    "max PSII yield (k_npq_r=0, all RC open)"
    Φ_PSII_MAX::FT

    # prognostic variables that change with time
    "Reversible NPQ rate constant (initially zero)"
    k_npq_rev::FT
    "Sustained NPQ rate constant (for seasonal changes, default is zero)"
    k_npq_sus::FT
    "Rate constant for photochemistry"
    k_p::FT
    "Non-Photochemical quenching "
    npq::FT
    "Fluorescence yield"
    ϕ_f::FT
    "Photochemical yield"
    ϕ_p::FT

    # dignostic variables that change with time
    "Dark adapted yield (`Kp=0`)"
    f_m::FT
    "Light adapted yield (`Kp=0`)"
    f_m′::FT
    "Dark-adapted fluorescence yield (`Kp=max`)"
    f_o::FT
    "Light-adapted fluorescence yield in the dark (`Kp=max`)"
    f_o′::FT
    "Energy quenching"
    q_e::FT
    "Photochemical quenching"
    q_p::FT
end

#######################################################################################################################################################################################################
#
# Changes to this constructor
# General
#     2022-Jan-15: isolate the reaction center from Leaf in Photosynthesis.jl
#     2022-Jan-25: fix documentation
#
#######################################################################################################################################################################################################
"""

    VJPReactionCenter{FT}() where {FT<:AbstractFloat}

Constructor of `VJPReactionCenter`

---
# Examples
```julia
rc = VJPReactionCenter{Float64}();
```
"""
VJPReactionCenter{FT}() where {FT<:AbstractFloat} = (
    return VJPReactionCenter{FT}(
                0.5,                                # F_PSII
                VanDerTolFluorescenceModel{FT}(),   # FLM
                0.85,                               # K_D
                0.05,                               # K_F
                4,                                  # K_P_MAX
                4/(0.85+0.05+4),                    # Φ_PSII_MAX
                0,                                  # k_npq_rev
                0,                                  # k_npq_sus
                4,                                  # k_p
                0,                                  # npq
                0,                                  # ϕ_f
                0,                                  # ϕ_p
                0,                                  # f_m
                0,                                  # f_m′
                0,                                  # f_o
                0,                                  # f_o′
                0,                                  # q_e
                0)                                  # q_p
);

#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2022-Jan-18: add the struct of Cytochrome reaction center
#     2022-Jan-25: fix documentation
#     2022-Feb-07: add more fields to work with Photosynthesis v0.3.1
#     2022-Feb-10: add K_X, ϵ_1, and ϵ_2 fields
#     2022-Mar-01: add more fields: η_c and η_l
#     2022-Mar-01: delete Η_C and Η_L, move η_c and η_l to photosynthesis model
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Structure that stores reaction center information

# Fields

$(TYPEDFIELDS)

"""
mutable struct CytochromeReactionCenter{FT<:AbstractFloat} <:AbstractReactionCenter{FT}
    # parameters that do not change with time
    "Fraction of absorbed light used by PSI ETR"
    F_PSI::FT
    "Rate constant of consititutive heat loss from the antennae `[ns⁻¹]`"
    K_D::FT
    "Rate constant of fluorescence `[ns⁻¹]`"
    K_F::FT
    "Rate constant of photochemistry for PS I `[ns⁻¹]`"
    K_PSI::FT
    "Rate constant of photochemistry for PS II `[ns⁻¹]`"
    K_PSII::FT
    "Rate constant of excitation sharing for PS II `[ns⁻¹]`"
    K_U::FT
    "Rate constant of regulated heat loss via oxidized PS I center `[s⁻¹]`"
    K_X::FT
    "Maximal PS I photochemical yield"
    Φ_PSI_MAX::FT

    # prognostic variables that change with time
    "Weight factor that PSI fluorescence reaches sensor (after reabsorption)"
    ϵ_1::FT
    "Weight factor that PSII fluorescence reaches sensor (after reabsorption)"
    ϵ_2::FT
    "Fluorescence yield"
    ϕ_f::FT
    "Photochemical yield"
    ϕ_p::FT

    # dignostic variables that change with time
end

#######################################################################################################################################################################################################
#
# Changes to this constructor
# General
#     2022-Jan-24: add constructor for CytochromeReactionCenter
#     2022-Jan-25: fix documentation
#     2022-Feb-07: make the constructor more readable
#     2022-Feb-10: add K_X, ϵ_1, and ϵ_2 fields
#     2022-Mar-01: add more fields: η_c and η_l
#     2022-Mar-01: delete Η_C and Η_L, move η_c and η_l to photosynthesis model
#
#######################################################################################################################################################################################################
"""
    CytochromeReactionCenter{FT}() where {FT<:AbstractFloat}

Constructor of `CytochromeReactionCenter`

---
# Examples
```julia
rc = CytochromeReactionCenter{Float64}();
```
"""
CytochromeReactionCenter{FT}() where {FT<:AbstractFloat} = (
    return CytochromeReactionCenter{FT}(
                0.41/(0.41+0.44),           # F_PSI
                0.55,                       # K_D
                0.05,                       # K_F
                14.5,                       # K_PSI
                4.5,                        # K_PSII
                0,                          # K_U
                14.5,                       # K_X
                14.5 / (14.5+0.55+0.05),    # Φ_PSI_MAX
                0,                          # ϵ_1
                1,                          # ϵ_2
                0,                          # ϕ_f
                0)                          # ϕ_p
);


text:
import Test: @test, @testset, @test_throws
import Luna
import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes

import LinearAlgebra: norm
@testset "mode-average vs modal" begin
## setup
a = 13e-6
gas = :Ar
pres = 5
τ = 30e-15
λ0 = 800e-9
L = 10e-2
grid = Grid.RealGrid(L, 800e-9, (160e-9, 3000e-9), 0.5e-12)

a0 = a
aL = 3a/4
afun = let a0=a0, aL=aL, L=L
    afun(z) = a0 + (aL-a0)*z/L
end

dens0 = PhysData.density(gas, pres)
densityfun(z) = dens0
energyfun, energyfunω = Fields.energyfuncs(grid)
responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)

## mode-average
Iωavg = let
m = Capillary.MarcatiliMode(afun, gas, pres, loss=false, model=:full);
aeff(z) = Modes.Aeff(m, z=z)
linop, βfun! = LinearOps.make_linop(grid, m, λ0)
inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=600e-9)
Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output, status_period=10)
abs2.(output["Eω"])
end

## modal
Iωmodal = let
modes = (
    Capillary.MarcatiliMode(afun, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
)
linop = LinearOps.make_linop(grid, modes, λ0)
inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=600e-9)
Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                               modes, :y, full=false)
statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output, status_period=10)
abs2.(dropdims(output["Eω"], dims=2))
end

@test norm(Iωavg - Iωmodal)/(sqrt(norm(Iωavg))*sqrt(norm(Iωmodal))) < 0.004
end

@testset "const vs afun" begin
## setup
a = 13e-6
gas = :Ar
pres = 5
τ = 30e-15
λ0 = 800e-9
L = 5e-2
grid = Grid.RealGrid(L, 800e-9, (160e-9, 3000e-9), 0.5e-12)

afun = let a=a
    z -> a
end

dens0 = PhysData.density(gas, pres)
densityfun(z) = dens0
energyfun, energyfunω = Fields.energyfuncs(grid)
responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)

## mode-average
Iωavg = let
m = Capillary.MarcatiliMode(afun, gas, pres, loss=false, model=:full);
aeff(z) = Modes.Aeff(m, z=z)
linop, βfun! = LinearOps.make_linop(grid, m, λ0)
inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output, status_period=10)
abs2.(output["Eω"])
end

## mode-average, constant a
Iωavg_c = let
m = Capillary.MarcatiliMode(a, gas, pres, loss=false, model=:full);
aeff(z) = Modes.Aeff(m, z=z)
inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)
Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
Luna.run(Eω, grid, linop, transform, FT, output, status_period=10)
abs2.(output["Eω"])
end

@test all(Iωavg .≈ Iωavg_c)
end

ext: jl


text:
using MultivariateStats
using Base.Test

## testing data

function pwdists(X)
	S = Base.sumabs2(X, 1)
	D2 = S .+ S' - 2 * (X'X)
	D2[diagind(D2)] = 0.0
	return sqrt(D2)
end

X0 = randn(3, 5)
G0 = X0'X0
D0 = pwdists(X0)

## conversion between dmat and gram

@assert issym(D0)
@assert issym(G0)

D = gram2dmat(G0)
@test issym(D)
@test_approx_eq D D0

G = dmat2gram(D0)
@test issym(G)
@test_approx_eq gram2dmat(G) D0

## classical MDS

X = classical_mds(D0, 3)
@test_approx_eq pwdists(X) D0


text:
using Memoize
@memoize ack3(m::Integer, n::Integer) = m == 0 ? n + 1 : n == 0 ? ack3(m - 1, 1) : ack3(m - 1, ack3(m, n - 1))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using LAR
include("src/inters.jl")

lineArray = randomLines(400,0.3);
W,EW = lines2lar(lineArray);
model = W,EW

printLar("test/csv/test1", model)
V,EV = readLar("test/csv/test1")
larview(V,EV)
viewexploded(V,EV)

EW = biconnectedComponents(V,EV)
viewexploded(V,EW)


text:
"""
Sun, W., Qu, J., Chen, Y., Di, Y. and Gao, F., 2016.
Heuristic sample reduction method for support vector data description.
Turkish Journal of Electrical Engineering & Computer Sciences, 24(1), pp.298-312.

Heuristic Sampling Reduction (HSR)
Parameter:
k
eps: 0.1 * num_dimensions
"""
mutable struct HSR <: Sampler
    k::Int
    eps::Real
    function HSR(k::Int, eps::Real=0.1)
        k > 0 || throw(ArgumentError("Invalid number of clusters $(k)."))
        eps >= 0 || throw(ArgumentError("Negative threshold parameter $(eps)."))
        return new(k, eps)
    end
end

function sample(sampler::HSR, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    sampler.k = min(size(data, 2), sampler.k)
    threshold = sampler.eps * size(data, 1)
    clustering = kmeans(data, sampler.k)
    dist_center_sorted_idx = sortperm(clustering.costs)
    tree = KDTree(data)
    sample_mask = trues(size(data, 2))
    processed = falses(size(data, 2))
    for i in dist_center_sorted_idx
        if !processed[i]
            # inrange includes i itself
            # -> query sorted result and drop first element
            # Special case: another observation is equal to i
            # -> one of the two remains in sample_mask
            to_remove = inrange(tree, data[:, i], threshold, true)[2:end]
            sample_mask[to_remove] .= false
            processed[to_remove] .= true
            processed[i] = true
        end
    end
    @assert all(processed)
    return sample_mask
end


text:
function vcov(X::AbstractMatrix, v::HC)
    N, p = size(X)
    XX = Base.LinAlg.At_mul_B(X, X)
    return scale!(XX, 1/N)
end

typealias GenLinMod GeneralizedLinearModel

adjfactor!(u, l::LinPredModel, k::HC0) = u[:] = one(Float64)
adjfactor!(u, l::LinPredModel, k::HC1) = u[:] = _nobs(l)./_df_residual(l)
adjfactor!(u, l::LinPredModel, k::HC2) = u[:] = 1./(1.-hatmatrix(l))
adjfactor!(u, l::LinPredModel, k::HC3) = u[:] = 1./(1.-hatmatrix(l)).^2

_nobs(l::LinPredModel) = length(l.rr.y)
_df_residual(l::LinPredModel) = _nobs(l) - length(coef(l))

function adjfactor!(u, lp::LinPredModel, k::HC4)
    h = hatmatrix(lp)
    n = _nobs(lp)
    p = npars(lp)
    for j = 1:n
        delta = min(4, n*h[j]/p)
        @inbounds u[j] = 1/(1-h[j])^delta
    end
end

function adjfactor!(u, lp::LinPredModel, k::HC4m)
    h = hatmatrix(lp)
    n = _nobs(lp)
    p = npars(lp)
    for j = 1:n
        delta = min(1.0, n*h[j]/p) + min(1.5, n*h[j]/p)
        @inbounds u[j] = 1/(1-h[j])^delta
    end
end

function adjfactor!(u, lp::LinPredModel, k::HC5)
    h     = hatmatrix(lp)
    n     = _nobs(lp)
    p     = npars(lp)
    mx    = max(n*0.7*maximum(h)/p, 4)
    for j = 1:n
        alpha =  min(n*h[j]/p, mx)
        @inbounds u[j] = 1/(1-h[j])^alpha
    end
end

nclus(k::CRHC) = length(unique(k.cl))
npars(x::LinPredModel) = length(x.pp.beta0)

function bread(lp::LinPredModel)
    A = inv(cholfact(lp.pp))
    scale!(A, nobs(lp))
end

residuals(l::LinPredModel, k::HC) = residuals(l)
residuals(l::LinPredModel, k::HAC) = residuals(l)
#residuals(l::LinPredModel, k::HC)  = wrkresidwts(l.rr)
#residuals(l::LinPredModel, k::HAC) = wrkresidwts(l.rr)

function residuals(r::GLM.ModResp)
    a = wrkwts(r)
    u = copy(wrkresid(r))
    length(a) == 0 ? u : broadcast!(*, u, u, a)
end

function wrkresidwts(r::GLM.ModResp)
    a = wrkwts(r)
    u = copy(wrkresid(r))
    length(a) == 0 ? u : broadcast!(*, u, u, a)
end

wrkresid(r::GLM.ModResp) = r.wrkresid
wrkwts(r::GLM.ModResp) = r.wrkwts
wrkwts(r::GLM.LmResp) = r.wts
wrkresid(r::GLM.LmResp) = r.y-r.mu

function weightedModelMatrix(l::LinPredModel)
    w = wrkwts(l.rr)
    length(w) > 0 ? ModelMatrix(l).*sqrt(w) : copy(ModelMatrix(l))
end

function hatmatrix(l::LinPredModel)
    # w = wrkwts(l.rr)
    # z = ModelMatrix(l).*sqrt(w)
    z = weightedModelMatrix(l)
    cf = cholfact(l.pp)[:UL]
    Base.LinAlg.A_rdiv_B!(z, cf)
    diag(Base.LinAlg.A_mul_Bt(z, z))
end

function vcov(ll::LinPredModel, k::HC)
    B = meat(ll, k)
    A = bread(ll)
    scale!(A*B*A, 1/nobs(ll))
end

function meat(l::LinPredModel,  k::HC)
    u = residuals(l, k)
    X = ModelMatrix(l)
    z = X.*u
    adjfactor!(u, l, k)
    scale!(Base.LinAlg.At_mul_B(z, z.*u), 1/nobs(l))
end

vcov(x::DataFrameRegressionModel, k::RobustVariance) = vcov(x.model, k)
stderr(x::DataFrameRegressionModel, k::RobustVariance) = sqrt(diag(vcov(x, k)))
stderr(x::LinPredModel, k::RobustVariance) = sqrt(diag(vcov(x, k)))

################################################################################
## Cluster
################################################################################

function clusterize!(M, U, bstarts)
    k, k = size(M)
    s = Array(Float64, k)
    V = Array(Float64, k, k)
    for m = 1:length(bstarts)
        for j = 1:k, i = 1:k
                @inbounds V[i, j] = zero(Float64)
        end
        for i = 1:k
            @inbounds s[i] = zero(Float64)
        end
        for j = 1:k, i = bstarts[m]
            @inbounds s[j] += U[i, j]
        end
        for j = 1:k, i = 1:k
            @inbounds V[i, j] += s[i]*s[j]
        end
        for j = 1:k, i = 1:k
            @inbounds M[i, j] += V[i, j]
        end
    end
end

function getqii(v::CRHC3, e, X, A, bstarts)
    for j in 1:length(bstarts)
        rnge = bstarts[j]
        se = view(e, rnge)
        sx = view(X, rnge, :)
        In = eye(length(rnge))
        e[rnge] =  (In - sx*A*sx')\se
    end
    return e
end

if VERSION < v"0.5.0"
    Base.cholfact(A::Symmetric, args...) = cholfact(A.data, Symbol(A.uplo), args...)
end

function getqii(v::CRHC2, e, X, A, bstarts)
    for j in 1:length(bstarts)
        rnge = bstarts[j]
        se = view(e, rnge)
        sx = view(X, rnge,:)
        In = eye(length(rnge))
        BB = Symmetric(In - sx*A*sx')
        e[rnge] =  cholfact(BB)\se
    end
    return e
end

function _adjresid!(v::CRHC, X, e, chol, bstarts)
    getqii(v, e, X, chol, bstarts)
end

_adjresid!(v::CRHC, X, e, ichol, bstarts, c::Float64) = scale!(c, _adjresid!(v::CRHC, X, e, ichol, bstarts))

function scalar_adjustment(X, bstarts)
    n, k = size(X);
    g    = length(bstarts);
    sqrt((n-1)/(n-k) * g/(g-1))
end

adjresid!(v::CRHC0, X, e, ichol, bstarts) = identity(e)
adjresid!(v::CRHC1, X, e, ichol, bstarts) = e[:] = scalar_adjustment(X, bstarts)*e
adjresid!(v::CRHC2, X, e, ichol, bstarts) = _adjresid!(v, X, e, ichol, bstarts, 1.0)
adjresid!(v::CRHC3, X, e, ichol, bstarts) = _adjresid!(v, X, e, ichol, bstarts, scalar_adjustment(X, bstarts))

function meat(x::LinPredModel, v::CRHC)
    idx = sortperm(v.cl)
    cls = v.cl[idx]
    #ichol = inv(x.pp.chol)
    ichol  = inv(cholfact(x.pp))
    X = ModelMatrix(x)[idx,:]
    e = wrkresid(x.rr)[idx]
    w = wrkwts(x.rr)
    if !isempty(w)
        w = w[idx]
        broadcast!(*, X, X, sqrt(w))
        broadcast!(*, e, e, sqrt(w))
    end
    bstarts = [searchsorted(cls, j[2]) for j in enumerate(unique(cls))]
    adjresid!(v, X, e, ichol, bstarts)
    M = zeros(size(X, 2), size(X, 2))
    clusterize!(M, X.*e, bstarts)
    return scale!(M, 1/nobs(x))
end

function vcov(x::LinPredModel, v::CRHC)
    B = bread(x)
    M = meat(x, v)
    scale!(B*M*B, 1/nobs(x))
end


text:
import TypeDBClient.typedb.protocol as Proto
using TypeDBClient: Label, bytes
using UUIDs: uuid4

@testset "RequestBuilder" begin

no_option = Proto.Options()

let b = TypeDBClient.DatabaseManagerRequestBuilder
    @test b.create_req("a").name == "a"
    @test b.contains_req("a").name == "a"
    @test_nowarn b.all_req()
end

let b = TypeDBClient.DatabaseRequestBuilder
    @test b.schema_req("a").name == "a"
    @test b.delete_req("a").name == "a"
end

let b = TypeDBClient.SessionRequestBuilder
    let r = b.open_req("a", Int32(0), no_option)
        @test r.database == "a"
        @test r._type == Int32(0)
        @test r.options == no_option
    end
    @test b.pulse_req([0x01]).session_id == [0x01]
end

let b = TypeDBClient.TransactionRequestBuilder
    let ts = [Proto.Transaction_Req()]
        @test b.client_msg(ts).reqs == ts
    end
    let id = uuid4(), r = b.stream_req(id)
        @test r.req_id == string(id)
        @test hasproperty(r, :stream_req)
    end
    let id = [0x01], ty = Int32(0), options = no_option, ms = 0,
            r = b.open_req(id, ty, options, ms)
        @test r.open_req.session_id == id
        @test r.open_req._type == ty
        @test r.open_req.options == no_option
        @test r.open_req.network_latency_millis == ms
    end
end

let b = TypeDBClient.QueryManagerRequestBuilder
    for f in (
        :define_req, :undefine_req, :match_req, :match_aggregate_req,
        :match_group_req, :match_group_aggregate_req,
        :insert_req, :delete_req, :update_req,
    )
        func = getfield(b, f)
        let q = func("a").query_manager_req
            underlying_req = getproperty(q, f)
            @test underlying_req.query == "a"
        end
    end
end

let b = TypeDBClient.ConceptManagerRequestBuilder
    @test b.put_entity_type_req("a").concept_manager_req.put_entity_type_req.label == "a"
    @test b.put_relation_type_req("a").concept_manager_req.put_relation_type_req.label == "a"
    @test b.put_attribute_type_req("a", Int32(0)).concept_manager_req.put_attribute_type_req.label == "a"
    @test b.put_attribute_type_req("a", Int32(0)).concept_manager_req.put_attribute_type_req.value_type == Int32(0)
    @test b.get_thing_type_req("a").concept_manager_req.get_thing_type_req.label == "a"
    @test b.get_thing_req("a").concept_manager_req.get_thing_req.iid == bytes("a")
end

let b = TypeDBClient.LogicManagerRequestBuilder
    let r = b.put_rule_req("a", "b", "c")
        @test r.logic_manager_req.put_rule_req.label == "a"
        @test r.logic_manager_req.put_rule_req.when == "b"
        @test r.logic_manager_req.put_rule_req.then == "c"
    end
    @test b.get_rule_req("a").logic_manager_req.get_rule_req.label == "a"
    @test_nowarn b.get_rules_req()
end

# standard labels for testing *type request builders
label = Label("a")
scoped_label = Label("s","a")

let b = TypeDBClient.TypeRequestBuilder

    # For simplicity, conditional scope settting is tested here (`is_abstract_req1)
    # and will be excluded for other tests below.

    let without_scope = b.is_abstract_req(label),
        with_scope = b.is_abstract_req(scoped_label)
        @test without_scope.type_req.label == "a"
        @test hasproperty(without_scope.type_req, :type_is_abstract_req)
        @test !hasproperty(without_scope.type_req, :scope)

        @test with_scope.type_req.scope == "s"
        @test with_scope.type_req.label == "a"
        @test hasproperty(with_scope.type_req, :type_is_abstract_req)
    end

    @test b.set_label_req(label, "b").type_req.type_set_label_req.label == "b"
    @test hasproperty(b.get_supertypes_req(label).type_req, :type_get_supertypes_req)
    @test hasproperty(b.get_subtypes_req(label).type_req, :type_get_subtypes_req)
    @test hasproperty(b.get_supertype_req(label).type_req, :type_get_supertype_req)
    @test hasproperty(b.delete_req(label).type_req, :type_delete_req)
end

let b = TypeDBClient.RoleTypeRequestBuilder

    let r = b.proto_role_type(scoped_label, Proto.Type_Encoding.THING_TYPE)
        @test r.label == "a"
        @test r.encoding == Proto.Type_Encoding.THING_TYPE
    end

    # test assertion for scope requirement
    @test_throws AssertionError b.proto_role_type(label, Proto.Type_Encoding.THING_TYPE)

    @test hasproperty(b.get_relation_types_req(label).type_req, :role_type_get_relation_types_req)
    @test hasproperty(b.get_players_req(label).type_req, :role_type_get_players_req)
end

# TODO Need to continue developing unit tests for these
# module ThingTypeRequestBuilder
# module EntityTypeRequestBuilder
# module RelationTypeRequestBuilder
# module ThingRequestBuilder
# module RelationRequestBuilder
# module AttributeRequestBuilder
# module RuleRequestBuilder

end # RequestBuilder


text:
include("../intcode.jl")

tape = set_memory!(Tape(parse_data(readline(open("input.txt")))), 100000)

function ascii_interact!(tape)
    while true
        out = take!(tape.outch)
        if isa(out, AwaitingInput)
            input = readline()
            for ch in input
                put!(tape.inch, Int64(ch))
                out = take!(tape.outch)
                isa(out, AwaitingInput) || error("unexpected " * string(out))
            end
            put!(tape.inch, 10)
        elseif isa(out, CodeFinished)
            println("FINISHED")
            return
        else
            print(Char(out))
        end
    end
end

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
module REST

using URIs
using Downloads
using MbedTLS
using Dates
using Base64
using XMLDict

const API_VER = "2016-05-31"

const DEFAULT_KWARGS = Dict(:timeout=>300, :throw=>false)

mutable struct StandardHeaders
    content_encoding::Union{String,Nothing}
    content_language::Union{String,Nothing}
    content_length::Union{String,Nothing}
    content_md5::Union{String,Nothing}
    content_type::Union{String,Nothing}
    date::Union{String,Nothing}
    if_modified_since::Union{String,Nothing}
    if_match::Union{String,Nothing}
    if_none_match::Union{String,Nothing}
    if_unmodified_since::Union{String,Nothing}
    range::Union{String,Nothing}
end

function StandardHeaders(; kwargs...)
    obj = StandardHeaders(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)
    for (n,v) in kwargs
        setfield!(obj, n, v)
    end
    obj
end

struct ServiceResponse
    resp::Union{Downloads.Response,Downloads.RequestError}
    data::IOBuffer
end

struct ServiceRequestError <: Exception
    code::Symbol
    resp::ServiceResponse
end

struct ServiceRequest
    account::String
    verb::String
    resource::String
    standard_headers::StandardHeaders
    headers::Dict{String,String}
    data::Union{IO,Nothing}

    function ServiceRequest(account::String, verb::String, resource::String, data::Vector{UInt8}; headers...)
        ServiceRequest(account, verb, resource, IOBuffer(data); headers...)
    end

    function ServiceRequest(account::String, verb::String, resource::String, data::AbstractString; headers...)
        ServiceRequest(account, verb, resource, convert(Vector{UInt8}, codeunits(data)); headers...)
    end

    function ServiceRequest(account::String, verb::String, resource::String, data::Union{IO,Nothing}=nothing; headers...)
        hdrdate = Dates.format(now(Dates.UTC), Dates.RFC1123Format)
        endswith(hdrdate, "GMT") || (hdrdate *= " GMT")

        reqhdrs = Dict{String,String}("x-ms-date"=>hdrdate, "x-ms-version"=>API_VER)
        (verb == "PUT") && (data === nothing) && (reqhdrs["Content-Length"] = "0")

        stdhdrparams = Dict{Symbol,String}()
        for (n,v) in headers
            if (n in fieldnames(StandardHeaders))
                stdhdrparams[n] = v
                reqhdrs[to_http_header_name(n)] = v
            else
                reqhdrs[string(n)] = v
            end
        end
        stdhdr = StandardHeaders(; stdhdrparams...)

        new(account, verb, resource, stdhdr, reqhdrs, data)
    end
end

struct AccountKey
    account::String
    key::String
end

to_http_header_name(n::Symbol) = join(map(titlecase, split(string(n), '_')), '-')

sign_strip(str::String) = replace(strip(str), "\n"=>" ")
sign_hdr(nv::Pair{String,String}) = lowercase(sign_strip(nv[1])) * ":" * sign_strip(nv[2])
sign_hdr(nv::Pair{T1,T2}) where {T1<:AbstractString, T2<:AbstractString} = sign_hdr(String(nv[1])=>String(nv[2]))

function sign_sharedkey(req::ServiceRequest, key::String)
    iob = IOBuffer()

    println(iob, req.verb)

    for n in fieldnames(StandardHeaders)
        v = getfield(req.standard_headers, n)
        # special case rule for Content-Length header in api version 2015-02-21 and later
        (n === :content_length) && (v !== nothing) && !isempty(v) && (parse(Int, v) == 0) && (v = nothing)
        println(iob, v === nothing ? "" : v)
    end

    x_ms_headers = sort(map(sign_hdr, collect(filter((p)->startswith(lowercase(strip(p[1])), "x-ms-"), req.headers))))
    for v in x_ms_headers
        println(iob, v)
    end

    print(iob, canonicalize_resource(req.account, req.resource))

    signingstr = String(take!(iob))
    hmacsign = base64encode(digest(MD_SHA256, signingstr, base64decode(key)))
    req.headers["Authorization"] = "SharedKey $(req.account):$(hmacsign)"
    nothing
end

function canonicalize_resource(resource_account::String, uri::String)
    url = URIs.URI(uri)
    container = url.path
    query = url.query

    iob = IOBuffer()
    print(iob, "/", resource_account, container)
    if !isempty(query)
        qparts = sort(map(sign_hdr, collect(URIs.queryparams(url))))
        for q in qparts
            print(iob, "\n", q)
        end
    end
    String(take!(iob))
end

function execute(req::ServiceRequest, key::String; retry_count::Int=0, retry_interval::Int=60)
    success = false
    svcresp = nothing
    retry = true
    count = 0
    sign_sharedkey(req, key)
    while !success && retry
        count += 1
        (count == 1) || sleep(retry_interval)
        verb = uppercase(req.verb)
        uri = string(URIs.URI(req.resource))
        output = IOBuffer()
        if (verb in ("POST", "PUT")) && (req.data !== nothing)
            #@debug("HTTP request", verb, uri, req.headers, req.data, DEFAULT_KWARGS)
            resp = Downloads.request(uri; method=verb, headers=req.headers, input=req.data, output=output, DEFAULT_KWARGS...)
        else
            #@debug("HTTP request", verb, uri, req.headers, DEFAULT_KWARGS)
            resp = Downloads.request(uri; method=verb, headers=req.headers, output=output, DEFAULT_KWARGS...)
        end
        svcresp = ServiceResponse(resp, output)
        retry = (count <= retry_count)
        issuccess(svcresp) && return svcresp
        @warn("Storage service request failed", resp, retry)
    end

    return svcresp
end

issuccess(resp::ServiceResponse) = isa(resp.resp, Downloads.Response) && (200 <= resp.resp.status <= 206)

parse_service_response(resp::ServiceResponse) = parse_xml(String(service_response(resp)))

function service_response(resp::ServiceResponse)
    issuccess(resp) || throw(ServiceRequestError(:http_error, resp))
    take!(resp.data)
end

end # module REST


text:
blp_date_string(date::Date) :: String = Dates.format(date, dateformat"yyyymmdd")
blpstring(str::AbstractString) = str
blpstring(date::Date) = blp_date_string(date)

function SchemaElementDefinition(element::Element)
    return SchemaElementDefinition(blpapi_Element_definition(element.handle))
end

function is_array(::Type{Element{A,D}}) :: Bool where {A,D}
    A
end

function is_array(::Element{A,D}) :: Bool where {A,D}
    A
end

@generated function get_child_element(element::AbstractElement{A,D}, index) where {A,D}
    if is_complex_datatype(D)

        @assert !A "get_child_element not allowed for array element."

        if index <: Integer
            # index into element using an Integer as key
            return quote
                result_element_handle_ref = Ref{Ptr{Cvoid}}(C_NULL)
                err = blpapi_Element_getElementAt(element.handle, result_element_handle_ref, index - 1)
                error_check(err, "Failed to get element at index $index from $(element.name)")
                return Element(result_element_handle_ref[], element)
            end

        elseif index <: AbstractString
            # index into element using a String as key
            return quote
                result_element_handle_ref = Ref{Ptr{Cvoid}}(C_NULL)
                err = blpapi_Element_getElement(element.handle, result_element_handle_ref, pointer(index), C_NULL)
                error_check(err, "Failed to get element $index from $(element.name)")
                return Element(result_element_handle_ref[], element)
            end
        else
            error("index data type $index not supported")
        end
    end

    error("get_child_element not allowed for data type $D.")
end

@generated function Base.getindex(element::AbstractElement{false,D}, index::AbstractString) where {D}

    # index <: AbstractString
    # used to index into child elements of a complex element that is not an array
    if is_complex_datatype(D)
        return quote
            return get_child_element(element, index)
        end
    end

    error("Base.getindex not implemented for element with datatype $D with is_array = $A.")
end

has_name(element::Element, name) = has_name(element.name, name)
get_name(element::Element) = element.name

function get_choice(element::AbstractElement{false,BLPAPI_DATATYPE_CHOICE})
    result_element_handle_ref = Ref{Ptr{Cvoid}}(C_NULL)
    err = blpapi_Element_getChoice(element.handle, result_element_handle_ref)
    error_check(err, "Failed to get choice from $(get_name(element))")
    return Element(result_element_handle_ref[], element)
end

Base.getindex(element::AbstractElement{false,BLPAPI_DATATYPE_CHOICE}) = get_choice(element)

@generated function has_child_element(element::AbstractElement{false,D}, name::AbstractString) :: Bool where {D}
    if is_complex_datatype(D)
        return quote
            return blpapi_Element_hasElement(element.handle, pointer(name), C_NULL) != 0
        end
    end

    error("has_child_element not implemented for element with datatype $D.")
end

function Base.haskey(element::AbstractElement{false,D}, name::AbstractString) :: Bool where {D}
    has_child_element(element, name)
end

struct ChildElementsIterator{T<:AbstractElement{false,BLPAPI_DATATYPE_SEQUENCE}}
    element::T
    num_elements::Csize_t
end

function each_child_element(element::T) where {T<:AbstractElement{false, BLPAPI_DATATYPE_SEQUENCE}}
    ChildElementsIterator(element, blpapi_Element_numElements(element.handle))
end

Base.isempty(itr::ChildElementsIterator) = itr.num_elements == 0
Base.length(itr::ChildElementsIterator) = itr.num_elements

function Base.iterate(itr::ChildElementsIterator)
    if isempty(itr)
        return nothing
    else
        # (current_element, next_index)
        return (get_child_element(itr.element, UInt64(1)), UInt64(2))
    end
end

function Base.iterate(itr::ChildElementsIterator, state::UInt64)
    if state > itr.num_elements
        return nothing
    else
        return (get_child_element(itr.element, state), state + UInt64(1))
    end
end

"""
    get_element_value(element, [index])

Returns the Julia value for `element`.

`element` cannot be a `BLPAPI_DATATYPE_SEQUENCE` or `BLPAPI_DATATYPE_CHOICE`.
"""
@generated function get_element_value(element::AbstractElement{A,D}) where {A,D}
    if A
        # retrieves a vector
        return quote
            num_values = blpapi_Element_numValues(element.handle)

            result = Vector()
            for i in 1:num_values
                push!(result, get_element_value(element, i))
            end

            return result
        end
    else
        if is_complex_datatype(D)
            error("get_element_value: Can't get element value from a complex element that is not an array.")
        end

        # retrieves a scalar
        return quote
            return get_element_value(element, 1)
        end
    end
end

@generated function get_element_value(element::AbstractElement{A,D}, index::Integer) where {A,D}
    if D == BLPAPI_DATATYPE_BOOL
        parse_value_block = quote
            buffer_ref = Ref{Cint}()
            err = blpapi_Element_getValueAsBool(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Bool element value")
            return buffer_ref[] != 0
        end

    elseif D == BLPAPI_DATATYPE_CHAR
        parse_value_block = quote
            buffer_ref = Ref{Cchar}()
            err = blpapi_Element_getValueAsChar(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Char element value")
            return Char(buffer_ref[])
        end

    elseif D == BLPAPI_DATATYPE_INT32
        parse_value_block = quote
            buffer_ref = Ref{Int32}()
            err = blpapi_Element_getValueAsInt32(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Int32 element value")
            return buffer_ref[]
        end

    elseif D == BLPAPI_DATATYPE_INT64
        parse_value_block = quote
            buffer_ref = Ref{Int64}()
            err = blpapi_Element_getValueAsInt64(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Int32 element value")
            return buffer_ref[]
        end

    elseif D == BLPAPI_DATATYPE_FLOAT32
        parse_value_block = quote
            buffer_ref = Ref{Float32}()
            err = blpapi_Element_getValueAsFloat32(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Float32 element value")
            return buffer_ref[]
        end

    elseif D == BLPAPI_DATATYPE_FLOAT64
        parse_value_block = quote
            buffer_ref = Ref{Float64}()
            err = blpapi_Element_getValueAsFloat64(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Float64 element value")
            return buffer_ref[]
        end

    elseif D == BLPAPI_DATATYPE_STRING
        parse_value_block = quote
            buffer_ref = Ref{Cstring}()
            err = blpapi_Element_getValueAsString(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get String element value")
            return unsafe_string(buffer_ref[])
        end

    elseif D == BLPAPI_DATATYPE_DATE
        parse_value_block = quote
            buffer_ref = Ref{BLPDateTime}()
            err = blpapi_Element_getValueAsDatetime(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Date element value")
            return Dates.Date(buffer_ref[])
        end

    elseif D == BLPAPI_DATATYPE_TIME
        parse_value_block = quote
            buffer_ref = Ref{BLPDateTime}()
            err = blpapi_Element_getValueAsDatetime(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Date element value")
            return Dates.Time(buffer_ref[])
        end

    elseif D == BLPAPI_DATATYPE_DATETIME
        parse_value_block = quote
            buffer_ref = Ref{BLPDateTime}()
            err = blpapi_Element_getValueAsDatetime(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Date element value")
            #return DateTime(buffer_ref[]) # TODO: use TimeZones.
            return buffer_ref[]
        end

    elseif D == BLPAPI_DATATYPE_SEQUENCE
        # get_element_value can be used on sequence only if it is an array
        @assert A "get_element_value: Can't get element value from datatype $D that is not an array (is_array = $A)."

        parse_value_block = quote
            buffer_ref = Ref{Ptr{Cvoid}}(C_NULL)
            err = blpapi_Element_getValueAsElement(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get Date element value")
            return Element(buffer_ref[], element)
        end

    elseif D == BLPAPI_DATATYPE_ENUMERATION
        parse_value_block = quote
            buffer_ref = Ref{Ptr{Cvoid}}(C_NULL)
            err = blpapi_Element_getValueAsName(element.handle, buffer_ref, index-1)
            error_check(err, "Failed to get BLPName element value from Enumeration")
            return Symbol(BLPName(buffer_ref[]))
        end

    else
        error("get_element_value: support for datatype $D not implemented.")
    end

    return quote
        if is_null_element(element)
            return missing
        else
            $parse_value_block
        end
    end
end

@generated function set_element_value!(element::AbstractElement{A,D}, val::V, index::Integer) where {A,D,V}
    if D == BLPAPI_DATATYPE_STRING
        return quote
            err = blpapi_Element_setValueString(element.handle, blpstring(val), index)
            error_check(err, "Failed to push value $val to element $(get_name(element))")
        end

    elseif D == BLPAPI_DATATYPE_BOOL
        @assert val == Bool "Can't set value of type $val to element of type $D."
        return quote
            err = blpapi_Element_setValueBool(element.handle, val, index)
            error_check(err, "Failed to push value $val to element $(get_name(element))")
        end

    elseif D == BLPAPI_DATATYPE_INT32
        @assert val <: Integer
        return quote
            err = blpapi_Element_setValueInt32(element.handle, Int32(val), index)
            error_check(err, "Failed to push value $val to element $(get_name(element))")
        end

    elseif D == BLPAPI_DATATYPE_INT64
        @assert val <: Integer
        return quote
            err = blpapi_Element_setValueInt64(element.handle, Int64(val), index)
            error_check(err, "Failed to push value $val to element $(get_name(element))")
        end

    elseif D == BLPAPI_DATATYPE_FLOAT32
        @assert val <: Number
        return quote
            err = blpapi_Element_setValueFloat32(element.handle, Float32(val), index)
            error_check(err, "Failed to push value $val to element $(get_name(element))")
        end

    elseif D == BLPAPI_DATATYPE_FLOAT64
        @assert val <: Number
        return quote
            err = blpapi_Element_setValueFloat64(element.handle, Float64(val), index)
            error_check(err, "Failed to push value $val to element $(get_name(element))")
        end

    elseif D == BLPAPI_DATATYPE_CHAR
        @assert val == Char
        return quote
            err = blpapi_Element_setValueChar(element.handle, val, index)
            error_check(err, "Failed to push value $val to element $(get_name(element))")
        end

    elseif D == BLPAPI_DATATYPE_ENUMERATION
        # for enums, we need to check the value type informed by the element's schema
        # and also if the value is allowed for the field
        return quote
            schema = SchemaElementDefinition(element)
            schema_type = schema.schema_type
            @assert isa(schema_type, EnumerationSchemaTypeDefinition)

            if schema_type.enumeration.datatype == BLPAPI_DATATYPE_STRING
                val_str = blpstring(val)
                @assert is_value_allowed(element, val_str) "Unvalid value for enum `$(get_name(element))`: `$val_str`. Options are: $(list_enum_options(element))"

                err = blpapi_Element_setValueString(element.handle, val_str, index)
                error_check(err, "Failed to push value $val_str to element $(get_name(element))")
            else
                error("enumeration datatype not supported: $(schema_type.enumeration.datatype)")
            end
        end

    elseif D == BLPAPI_DATATYPE_DATETIME || D == BLPAPI_DATATYPE_DATE
        @assert val == BLPDateTime || val == DateTime || val == Date

        return quote
            err = blpapi_Element_setValueDatetime(element.handle, BLPDateTime(val), index)
            error_check(err, "Failed to push value $val to element $(get_name(element))")
        end
    end

    error("set_element_value! not implemented for datatype $D.")
end

# allows `element[] = val`, when element is not an array
function Base.setindex!(element::Element{false,D}, val::V) where {D,V}
    set_element_value!(element, val, 0)
end

function is_value_allowed(element::Element{A, BLPAPI_DATATYPE_ENUMERATION}, val) where {A}
    schema = SchemaElementDefinition(element)
    schema_type = schema.schema_type
    @assert isa(schema_type, EnumerationSchemaTypeDefinition)

    # check wether val is allowed for this enum
    val_is_allowed = false
    for enum_const in schema_type.enumeration.list
        if enum_const.value == val
            val_is_allowed = true
            break
        end
    end

    return val_is_allowed
end

function list_enum_options(element::Element{A, BLPAPI_DATATYPE_ENUMERATION}) where {A}
    schema = SchemaElementDefinition(element)
    schema_type = schema.schema_type
    @assert isa(schema_type, EnumerationSchemaTypeDefinition)
    return [ enum_const.value for enum_const in schema_type.enumeration.list ]
end

function Base.push!(element::Element{true,D}, val::V) where {D,V}
    set_element_value!(element, val, BLPAPI_ELEMENT_INDEX_END)
end

function Base.push!(element::Element, vals...)
    for val in vals
        push!(element, val)
    end
end

Base.append!(element::Element{true}, vals::Vector) = push!(element, vals...)

# allows `element["child_element_name"] = val`, when element is complex and not an array
@generated function Base.setindex!(element::Element{false,D}, val::V, name::AbstractString) where {D,V}
    if is_complex_datatype(D)
        return quote
            child_element = element[name]
            child_element[] = val
            nothing
        end
    end

    error("setindex! not implemented for datatype $D.")
end

function is_null_value(element::AbstractElement, index::Integer)
    rc = blpapi_Element_isNullValue(element.handle, index-1)

    if rc != 0 && rc != 1
        error("Unexpected result of blpapi_Element_isNullValue on element of type $(typeof(element)): $rc")
    end

    return rc != 0
end

function is_null_element(element::AbstractElement)
    rc = blpapi_Element_isNull(element.handle)
    return rc != 0
end

function as_dict!(dict::Dict, root::Element{A,D}) where {A,D}

    if is_choice_datatype(D)
        as_dict!(dict, get_choice(root))

    elseif D == BLPAPI_DATATYPE_SEQUENCE

        # vector of sequences
        if A
            elements_vec = Vector()
            root_name_symbol = Symbol(get_name(root))
            dict[root_name_symbol] = elements_vec
            for el in get_element_value(root)
                el_dict = Dict()
                as_dict!(el_dict, el)
                if length(el_dict) == 1 && haskey(el_dict, root_name_symbol)
                    push!(elements_vec, el_dict[root_name_symbol])
                else
                    push!(elements_vec, el_dict)
                end
            end
        else
            # a single sequence
            child_dict = Dict()
            dict[Symbol(get_name(root))] = child_dict
            for child_element in each_child_element(root)
                as_dict!(child_dict, child_element)
            end
        end

    else
        @assert is_simple_datatype(D)
        dict[Symbol(get_name(root))] = get_element_value(root)
    end
end

function Base.Dict(element::AbstractElement)
    dict = Dict()
    as_dict!(dict, element)
    return dict
end

function push_named_tuples!(result::T, element_vec::Element{true, BLPAPI_DATATYPE_SEQUENCE}) where {T<:Vector}
    tuple_keys = nothing

    for element in get_element_value(element_vec)

        if tuple_keys == nothing
            tuple_keys = Tuple([ Symbol(get_name(child_element)) for child_element in each_child_element(element) ])
        else
            @assert tuple_keys == Tuple([ Symbol(get_name(child_element)) for child_element in each_child_element(element) ])
        end

        tuple_values = Tuple([ get_element_value(child_element) for child_element in each_child_element(element) ])
        push!(result, (; zip(tuple_keys, tuple_values)...)) # trick based on the docstring for NamedTuple
    end
end

to_named_tuple(element::Element{false, BLPAPI_DATATYPE_SEQUENCE}) = (; [ (Symbol(get_name(child_element)), get_element_value(child_element)) for child_element in each_child_element(element) ]...)

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
using TestFilter
using Test

@testset "TestFilter.jl" begin
    # Write your tests here.
    @test 1 + 1 == 2
    @test "qwe" == "qwe"

    x = add_one(10)

    @test x == 11

end


text:
using GpSelection

my_tests = ["test_GpSelection.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
@diffusion_process LorenzAux{K,R} begin
    :dimensions
    process --> 3
    wiener --> 3

    :parameters
    _ --> (3, K)
    σ --> K

    :auxiliary_info
    t0 --> Float64
    T --> Float64
    vT --> R

    :additional
    constdiff --> true
    linear --> true
end

DiffusionDefinition.B(t, P::LorenzAux) = @SMatrix [
    -P.p1         P.p1        0.0;
     P.p2-P.vT[3]  -1.0    -P.vT[1];
     P.vT[2]      P.vT[1]    -P.p3
]

DiffusionDefinition.β(t, P::LorenzAux) = @SVector[
    0.0,
    P.vT[1]*P.vT[3],
    -P.vT[1]*P.vT[3]
]

DiffusionDefinition.σ(t, P::LorenzAux) = SDiagonal(P.σ, P.σ, P.σ)


text:
# Sum square difference

LIMIT = 100

function solve()
    return sum(1:LIMIT)^2 - sum(n^2 for n in 1:LIMIT)
end

print(Int(solve()))

max_stars_repo_name: wephy/project-euler

max_issues_repo_name: wephy/project-euler

max_forks_repo_name: wephy/project-euler


text:
export stable_subgroups

add_verbose_scope(:StabSub)
add_assert_scope(:StabSub)

###############################################################################
#
#  Tools for ZpnGModules
#
###############################################################################

function base_ring(M::ZpnGModule)
  return M.R
end

function show(io::IO, M::ZpnGModule)
  print(io, "Module over Z/", M.R.n, "Z with structure ", M.V)
end

#
#  Lifts a matrix from F_p to Z/p^nZ
#

function lift(M::fq_nmod_mat, R::Nemo.NmodRing)
  @hassert :StabSub 1 isprime_power(modulus(R))
  N=zero_matrix(R,nrows(M),ncols(M))
  for i=1:nrows(M)
    for j=1:ncols(M)
      N[i,j]=FlintZZ(coeff(M[i,j],0))
    end
  end
  return N
end

function lift(M::gfp_mat, R::Nemo.NmodRing)
  @hassert :StabSub 1 isprime_power(modulus(R))
  N=zero_matrix(R, nrows(M), ncols(M))
  for i=1:nrows(M)
    for j=1:ncols(M)
      N[i,j] = R(lift(M[i,j]))
    end
  end
  return N
end

#  Action of a matrix on an element of the group
function *(x::GrpAbFinGenElem, M::nmod_mat)
  G = parent(x)
  @assert ngens(G) == nrows(M)
  R = base_ring(M)
  coeff = map_entries(R, x.coeff)
  y = coeff*M
  l = lift(y)
  return GrpAbFinGenElem(G, l)
end

#  Smith Normal Form for a ZpnGModule
function Nemo.snf(M::ZpnGModule)
  A = M.V
  G = M.G
  if issnf(A)
    return M, id_hom(A)
  end
  R = base_ring(M)
  S, mS = snf(A)
  image_mS = map_entries(R, mS.map)
  preimage_mS = map_entries(R, mS.imap)
  H = Array{nmod_mat,1}(undef, length(G))
  for i=1:length(G)
    H[i] = image_mS*G[i]*preimage_mS
  end
  return ZpnGModule(S, H), mS
end

function isstable(act::Array{T, 1}, mS::GrpAbFinGenMap) where T <: Map{GrpAbFinGen, GrpAbFinGen}

  S=mS.header.domain
  for s in gens(S)
    x=mS(s)
    for g in act
      if !haspreimage(mS,g(x))[1]
        return false
      end
    end
  end
  return true

end

function issubmodule(M::ZpnGModule, S::nmod_mat)

  #@assert issnf(M.V)
  s, ms = submodule_to_subgroup(M, S)
  for x in gens(s)
    el=ms(x)
    for g in M.G
      if !haspreimage(ms,el*g)[1]
        return false
      end
    end
  end
  return true

end

#
#  Given a group $G$ and a group of automorphisms of G, this function returns the corresponding ZpnGModule
#

function action(V::GrpAbFinGen, act::Array{T,1}) where T<: Map{GrpAbFinGen, GrpAbFinGen}

  expon = Int(exponent(V))
  @hassert :StabSub 1 length(factor(order(V)).fac)==1
  RR = ResidueRing(FlintZZ, expon, cached=false)
  act_mat = Array{nmod_mat,1}(undef, length(act))
  for z = 1:length(act)
    A = zero_matrix(RR, ngens(V), ngens(V))
    for i = 1:ngens(V)
      y = act[z](V[i])
      for j=1:ngens(V)
        A[i,j]=y.coeff[1,j]
      end
    end
    act_mat[z]=A
  end
  return ZpnGModule(V,act_mat)

end

#################################################################################
#
#  Duality
#
#################################################################################

#=
  To get the transpose map of a homomorphism, and so the action on the dual module,
  you need to transpose the matrix and multiply or the elements for a power of p, assuming that the group is in snf.
  Precisely, let p^e be the exponent of the group and let p^t_i be the powers of p appearing in the snf. Then
  ( p^(e-t_1)     )                  ( p^(e-t_1)     )
  (          ...  )  x ^t A  = A  x  (          ...  )
  (              1)                  (              1)
=#

function dual_module(M::ZpnGModule)

  @assert issnf(M.V)
  G1=deepcopy(M.G)
  for i=1:length(G1)
    for j=1:ngens(M.V)-1
      for k=j+1:ngens(M.V)
        x=divexact(M.V.snf[k], M.V.snf[j])
        G1[i][j,k]=divexact(G1[i][j,k].data, x)
        G1[i][k,j]*=x
      end
    end
    transpose!(G1[i])
  end
  return ZpnGModule(M.V, G1)

end

function _dualize(M::nmod_mat, V::GrpAbFinGen, v::Array{fmpz,1})
  #  First, compute the kernel of the corresponding homomorphisms
  K = abelian_group(fmpz[V.snf[end] for j=1:nrows(M)])
  A = lift(transpose(M))
  for j=1:nrows(A)
    for k=1:ncols(A)
      A[j, k] *= v[j]
    end
  end
  mH = Hecke.GrpAbFinGenMap(V, K, A)
  newel = kernel_as_submodule(mH)
  return change_base_ring(base_ring(M), newel)
end

function _dualize_1(M::nmod_mat, snf_struct::Array{fmpz,1})

  A=nullspace(transpose(M))
  B=vcat(transpose(A),zero_matrix(M[1,1].parent, ncols(A),ncols(A)))
  for j=1:ncols(A)
    B[nrows(A)+j,j]=snf_struct[j]
  end
  S=nullspace(B)
  C=vcat(transpose(A),zero_matrix(M[1,1].parent, ncols(A),ncols(A)))
  return S*C

end

function kernel_as_submodule(h::GrpAbFinGenMap)
  G = domain(h)
  H = codomain(h)
  hn, t = hnf_with_transform(vcat(h.map, rels(H)))
  for i=1:nrows(hn)
    if iszero_row(hn, i)
      return view(t, i:nrows(t), 1:ngens(G))
    end
  end
  error("JH")
end

########################################################################################################
#
#  Quotients and subgroups of ZpnGModules
#
########################################################################################################

function quo(M::ZpnGModule, S::nmod_mat)
  Q, mQ=quo(M.V, lift(S), false)
  return ZpnGModule(Q, M.G), mQ
end

function sub(M::ZpnGModule, S::nmod_mat)

  sg, msg=sub(M.V, lift(S), false)
  G=Array{nmod_mat,1}(undef, length(M.G))
  for k=1:length(M.G)
    A=zero_matrix(M.R, ngens(sg), ngens(sg))
    for i=1:ngens(sg)
      x=msg(sg[i])*M.G[k]
      x=haspreimage(msg, x)[2].coeff
      for j=1:ngens(sg)
        A[i,j]=x[1,j]
      end
    end
    G[k]=A
  end
  return ZpnGModule(sg,G), msg

end

function sub(M::ZpnGModule, n::Int)

  if issnf(M.V)
    return _sub_snf(M, n)
  end
  sg, msg = sub(M.V, n, false)
  R = base_ring(M)
  G = Array{nmod_mat,1}(undef, length(M.G))
  big_m = vcat(msg.map, rels(M.V))
  for k = 1:length(M.G)
    A = map_entries(R, msg.map)*M.G[k]
    fl, res = can_solve_with_solution(big_m, lift(A), side = :left)
    @assert fl
    G[k] = map_entries(R, view(res, 1:ngens(sg), 1:ngens(sg)))
  end
  return ZpnGModule(sg,G), msg

end

function _sub_snf(M::ZpnGModule, n::Int)
  V = M.V
  nfmpz = fmpz(n)
  ind = 1
  while V.snf[ind] <= nfmpz
    ind += 1
  end
  invariants = Vector{fmpz}(undef, length(V.snf)-ind+1)
  for i = ind:length(V.snf)
    invariants[i-ind+1] = divexact(V.snf[i], n)
  end
  Gnew = abelian_group(invariants)
  action = nmod_mat[sub(x, ind:ngens(V), ind:ngens(V)) for x in M.G]
  mat_map = zero_matrix(FlintZZ, length(invariants), ngens(V))
  for i = 1:ngens(Gnew)
    mat_map[i, ind+i-1] = n
  end
  mGnew = hom(Gnew, M.V, mat_map)
  return ZpnGModule(Gnew, action), mGnew
end

function _exponent_p_sub(M::ZpnGModule; F::GaloisField = GF(M.p, cached = false))

  @assert issnf(M.V)
  G = M.G
  V = M.V
  p = M.p
  v = fmpz[divexact(V.snf[i], p) for i=1:ngens(V)]
  G1 = Array{gfp_mat,1}(undef, length(G))
  for s=1:length(G1)
    G1[s] = zero_matrix(F, ngens(V), ngens(V))
    for i=1:ngens(V)
      for j=1:ngens(V)
        if G[s][i,j] !=0 && v[i] <= v[j]
          G1[s][i,j] = divexact((G[s][i,j].data)*v[i], v[j])
        end
      end
    end
  end
  return ModAlgAss(G1)

end

function submodule_to_subgroup(M::ZpnGModule, S::nmod_mat)
  return sub(M.V, lift(S), false)
end

##########################################################################
#
#  Minimal Submodules function
#
##########################################################################

function minimal_submodules(M::ZpnGModule, ord::Int=-1)

  R = M.R
  S, mS = snf(M)
  N = _exponent_p_sub(S)
  if ord == -1
    list_sub = minimal_submodules(N)
  else
    list_sub = minimal_submodules(N, ngens(S.V)-ord)
  end
  list = Array{nmod_mat,1}(undef, length(list_sub))
  v = Int[M.p^(valuation(S.V.snf[i], M.p)-1) for i=1:ngens(S.V)]
  W = MatrixSpace(R, 1, ngens(M.V), false)
  for z = 1:length(list)
    list[z] = vcat(nmod_mat[W((mS(S.V(fmpz[lift(list_sub[z][k,i])*v[i] for i=1:ngens(S.V)]))).coeff) for k=1:nrows(list_sub[z])])
  end
  return list

end

###########################################################################################
#
#  Maximal Submodules function
#
###########################################################################################

function maximal_submodules(M::ZpnGModule, ind::Int=-1)

  R = M.R
  S, mS = snf(M)
  N = dual_module(S)
  if ind == -1
    minlist = minimal_submodules(N)
  else
    minlist = minimal_submodules(N,ind)
  end
  list=Array{nmod_mat,1}(undef, length(minlist))
  v=[divexact(fmpz(R.n),S.V.snf[j]) for j=1:ngens(S.V) ]
  for x in minlist
    K = abelian_group([fmpz(R.n) for j=1:nrows(x)])
    A = lift(transpose(x))
    for j=1:nrows(A)
      for k=1:ncols(A)
        A[j,k]*=v[j]
      end
    end
    mH = Hecke.GrpAbFinGenMap(S.V,K,A)
    sg, msg = kernel(mH)
    push!(list, vcat([ (mS(msg(y))).coeff for y in gens(sg)]))
  end
  return list

end

##########################################################################################
#
#  Composition factors and series
#
##########################################################################################
#
#  Given a list of square matrices G, it returns a list of matrices given by the minors
#  (n-s) x (n-s) of the matrices G[i] mod p
#
function _change_ring(G::Array{nmod_mat,1}, F::GaloisField, s::Int)

  G1 = Array{gfp_mat, 1}(undef, length(G))
  n = nrows(G[1])
  for i = 1:length(G)
    M = zero_matrix(F,n-s+1,n-s+1)
    for j = s:n
      for k = s:n
        M[j-s+1,k-s+1] = F(deepcopy((G[i][j,k]).data))
      end
    end
    G1[i] = M
  end
  return G1

end

#
#  Cut the module in submodules with exponent p, returning the quotients p^i M /p^(i+1) M
#
function _mult_by_p(M::ZpnGModule)
  G = M.G
  V = M.V
  p = M.p
  @assert issnf(V)
  #  First, the quotient M/pM
  F = GF(p, cached = false)
  n = ngens(V)
  Gq = _change_ring(G, F, 1)
  spaces = [ModAlgAss(Gq)]
  #
  #  Now, the others
  #
  s = valuation(fmpz(M.R.n),p)
  j = 1
  for i = 2:s
    while !divisible(V.snf[j], p^i)
      j += 1
    end
    GNew = _change_ring(G, F, j)
    push!(spaces, ModAlgAss(GNew))
  end
  return spaces
end

function composition_factors(M::ZpnGModule; dimension::Int = -1)
  N, _ = snf(M)
  list_sub = _mult_by_p(N)
  list = Vector{Tuple{ModAlgAss{GaloisField, gfp_mat, gfp_elem}, Int}}()
  for x in list_sub
    append!(list, composition_factors(x, dimension = dimension))
  end
  final_list = Vector{Tuple{ModAlgAss{GaloisField, gfp_mat, gfp_elem}, Int}}()
  done = falses(length(list))
  for i = 1:length(list)
    if done[i]
      continue
    end
    done[i] = true
    mult = list[i][2]
    for j = i+1:length(list)
      if !done[j] && isisomorphic(list[i][1], list[j][1])
        mult += list[j][2]
        done[j] = true
      end
    end
    push!(final_list, (list[i][1], mult))
  end
  return final_list
end

###############################################################################
#
#  Submodules interface
#
###############################################################################

@doc Markdown.doc"""
    submodules(M::ZpnGModule; typequo, typesub, order)

Given a ZpnGModule $M$, the function returns all the submodules of $M$.

"""
function submodules(M::ZpnGModule; typequo=Int[-1], typesub=Int[-1], ord=-1)

  if typequo!=[-1]
    return submodules_with_quo_struct(M,typequo)
  elseif typesub!=[-1]
    return submodules_with_struct(M,typesub)
  elseif ord!=-1
    return submodules_order(M,ord)
  else
    return submodules_all(M)
  end

end

###############################################################################
#
#  Function to find all the submodules
#
###############################################################################

function submodules_all(M::ZpnGModule)

  R = M.R
  if isone(order(M.V))
    return (nmod_mat[])
  end
  S,mS = snf(M)
  minlist = minimal_submodules(S)
  list = nmod_mat[identity_matrix(R, length(S.V.snf)), zero_matrix(R,1,length(S.V.snf))]

  #
  #  Find minimal submodules, factor out and recursion on the quotients
  #

  for x in minlist
    N, _ = quo(S, x)
    newlist = submodules_all(N)
    for y in newlist
      push!(list,vcat(y,x))
    end
  end

  append!(list,minlist)
  #
  #  Eliminate redundance via Howell form
  #
  w=fmpz[divexact(fmpz(R.n), S.V.snf[j]) for j=1:ngens(S.V)]
  list=_no_redundancy(list, w)

  #
  #  Writing the submodules in terms of the given generators and returning an iterator
  #
  MatSnf=change_base_ring(R, mS.map)
  return (x*MatSnf for x in list)

end

###############################################################################
#
#  Submodules with given structure as a subgroup
#
###############################################################################

function main_submodules_cyclic(M::ZpnGModule, ord::Int)

  lf = composition_factors(M, dimension = 1)
  if isempty(lf)
    return Vector{nmod_mat}()
  end
  list = _submodules_with_struct_cyclic(M, ord, lf)
  for step = 1:ord-1
    c = fmpz(M.p^step)
    list1 = nmod_mat[]
    for x in list
      L, _ = quo(M, x)
      newlist = _submodules_with_struct_cyclic(L, ord-step, lf)
      for i = 1:length(newlist)
        mat_test = lift(newlist[i])
        mul!(mat_test, mat_test, c)
        el = GrpAbFinGenElem(M.V, mat_test)
        if !iszero(el)
          push!(list1, newlist[i])
        end
      end
    end
    list = list1
  end
  return list

end

function _submodules_with_struct_cyclic(M::ZpnGModule, ord::Int, lf)
  R = M.R
  #  We search for an element in p^(ord-1)*G
  s, ms = sub(M, M.p^(ord-1))
  S, mS = snf(s)
  N = _exponent_p_sub(S, F = base_ring(lf[1][1]))
  @vtime :StabSub 1 submod = minimal_submodules(N, 1, lf)
  list1 = Array{nmod_mat,1}(undef, length(submod))
  v = nmod[R(divexact(S.V.snf[i], M.p)) for i = 1:ngens(S.V)]
  for i = 1:length(submod)
    list1[i] = lift(submod[i], R)
    @assert nrows(list1[i]) == 1
    for k = 1:ncols(list1[i])
      list1[i][1, k] *= v[k]
    end
  end
  MatSnf = change_base_ring(R, mS.map*ms.map)
  for j=1:length(list1)
    list1[j] = list1[j]*MatSnf
  end
  return list1

end

function _update_typesub(typesub::Vector{Int})
  i = 1
  while typesub[i] != typesub[end]
    i += 1
  end
  new_typesub = Array{Int, 1}(undef, length(typesub))
  for j = 1:i-1
    new_typesub[j] = typesub[j]
  end
  for j = i:length(new_typesub)
    new_typesub[j] = typesub[j] - 1
  end
  return new_typesub
end

function _submodules_with_struct_main(M::ZpnGModule, typesub::Array{Int,1})
  @assert issnf(M.V)
  p = M.p
  R = base_ring(M)
  #First iteration out of the loop.
  list = _submodules_with_struct(M, typesub)
  @vprint :StabSub 1 "Subgroups at first step: $(length(list)) \n"
  #Some data needed in the loop
  w = Vector{fmpz}(undef, ngens(M.V))
  for i = 1:length(w)
    w[i] = divexact(fmpz(R.n), M.V.snf[i])
  end
  new_typesub = _update_typesub(typesub)
  #Now, the iterative process.
  for i = 1:typesub[end]-1
    list1 = nmod_mat[]
    new_typesub1 = _update_typesub(new_typesub)
    diag = typesub - new_typesub1
    Gtest = snf(abelian_group(Int[p^x for x in diag]))[1]
    order_test = order(Gtest)
    for x in list
      L, _ = quo(M, x)
      newlist = _submodules_with_struct(L, new_typesub)
      @vprint :StabSub 1 "Candidates to be added at this step: $(length(newlist)) \n"
      #First, I sieve for the subgroups with the correct structure.
      for s = 1:length(newlist)
        ord = fmpz(1)
	      for t = 1:nrows(newlist[s])
          ord *= order(GrpAbFinGenElem(M.V, lift(view(newlist[s], t:t, 1:ncols(newlist[s])))))
	      end
        if ord >= order_test
          t1, mt1 = submodule_to_subgroup(M, newlist[s])
          if _special_is_isomorphic(t1, Gtest)
	          push!(list1, newlist[s])
	        end
        end
      end
      @vprint :StabSub 1 "New till this point at this step: $(length(list1)) \n"
    end
    new_typesub = new_typesub1
    list = list1
  end
  return list
end

function _special_is_isomorphic(G::GrpAbFinGen, Gtest::GrpAbFinGen)
  if issnf(G)
    return G.snf == Gtest.snf
  end
  mat = hnf_modular_eldiv(G.rels, exponent(Gtest))
  mm = snf(mat)
  inv = fmpz[mm[i, i] for i = 1:ncols(mm) if !isone(mm[i, i])]
  return inv == Gtest.snf
end

function _submodules_with_struct(M::ZpnGModule, typesub::Array{Int, 1})

  R = base_ring(M)
  s, ms = sub(M, M.p^(typesub[end]-1))
  S, mS = snf(s)
  N = _exponent_p_sub(S)

  i = 1
  while typesub[i] != typesub[end]
    i += 1
  end
  a = length(typesub) - i + 1

  @vtime :StabSub 1 submod = submodules(N, (N.dimension)-a)
  #
  #  Write the modules as elements of S
  #

  list1 = Array{nmod_mat, 1}(undef, length(submod))
  v = fmpz[divexact(S.V.snf[i], M.p) for i = 1:ngens(S.V)]
  for i = 1:length(submod)
    @inbounds list1[i] = lift(submod[i], R)
    for j = 1:nrows(list1[i])
      for k = 1:ncols(list1[i])
        @inbounds list1[i][j,k] *= v[k]
      end
    end
  end
  #and now as elements of M
  auxmat = mS.map*ms.map
  auxmat2 = change_base_ring(R, auxmat)
  for j = 1:length(list1)
    @inbounds list1[j] = list1[j]*auxmat2
  end
  return list1

end

function submodules_with_struct(M::ZpnGModule, typesub::Array{Int,1})
  # If the group is cyclic, it is easier
  if length(typesub) == 1
    return main_submodules_cyclic(M, typesub[1])
  end
  sort!(typesub)
  return _submodules_with_struct_main(M, typesub)
end

function _no_redundancy(list::Array{nmod_mat,1}, w::Array{fmpz,1})

  R = base_ring(list[1])
  n = ncols(list[1])
  #
  #  Howell form of every candidate, embedding them in a free module
  #
  for i = 1:length(list)
    if n > nrows(list[i])
      @inbounds list[i] = vcat(list[i], zero_matrix(R, n-nrows(list[i]), ncols(list[i])))
    end
    for j=1:n
      for k=1:nrows(list[i])
        @inbounds list[i][k,j] *= w[j]
      end
    end
    howell_form!(list[i])
    list[i] = view(list[i], 1:n, 1:n)
  end

  #
  #  Now, check if they are equal
  #
  i=1
  list1 = nmod_mat[list[1]]
  for i = 2:length(list)
    found = false
    for j = 1:length(list1)
      if list[i] == list1[j]
        found = true
        break
      end
    end
    if !found
      push!(list1, list[i])
    end
  end

  #
  #  Write them again not embedded
  #
  for i=1:length(list1)
    for j=1:n
      for k=1:j
        list1[i][k,j] = divexact(list1[i][k,j].data, w[j])
      end
    end
  end
  return list1

end

###############################################################################
#
#  Submodules of given order
#
###############################################################################

function submodules_order(M::ZpnGModule, ord::Int)

  R=M.R
  S,mS=snf(M)
  @assert exponent(S.V)==fmpz(R.n)
  N=Hecke._exponent_p_sub(S)
  lf=composition_factors(N)
  list=nmod_mat[]
  v=fmpz[divexact(S.V.snf[i], M.p) for i=1:ngens(S.V)]
  for i=1:ord-1
    minlist=minimal_submodules(N,i,lf)
    for x in minlist
      A=lift(x,R)
      for s=1:nrows(A)
        for t=1:ncols(A)
          A[s,t]*=v[t]
        end
      end
      L, _=quo(S, A, false)
      newlist=Hecke.submodules_order(L,ord-i)
      for z=1:length(newlist)
        push!(list,vcat(newlist[z],A))
      end
    end
  end

  #
  #  Check for redundancy
  #
  w=fmpz[divexact(fmpz(R.n), S.V.snf[j]) for j=1:ngens(S.V)]
  list=_no_redundancy(list,w)

  #
  #  Write the submodules in terms of the set of given generators
  #

  MatSnf=map(mS.map, R)
  for j=1:length(list)
    list[j]=list[j]*MatSnf #vcat([W(( mS( S.V([list[j][k,i].data for i=1:ngens(S.V)]))).coeff)  for k=1:nrows(list[j])])
  end

  #
  #  Add the minimal submodules
  #

  minlist=minimal_submodules(N,ord, lf)
  for x in minlist
    push!(list, vcat([W((mS( S.V(fmpz[FlintZZ(coeff(x[k,i],0))*((M.p)^(v[i]-1)) for i=1:ngens(S.V)]))).coeff) for k=1:nrows(x) ]))
  end
  return (x for x in list)

end

###############################################################################
#
#  Submodules with given structure of the quotient
#
###############################################################################

function submodules_with_quo_struct(M::ZpnGModule, typequo::Array{Int,1})

  R = base_ring(M)
  p = M.p
  S, mS = snf(M)
  sort!(typequo)
  wish = abelian_group(Int[p^i for i in typequo])

  if isisomorphic(wish, S.V)
    return (nmod_mat[zero_matrix(R, 1, ngens(M.V))])
  end
  if length(typequo) > length(S.V.snf)
    return (nmod_mat[])
  end
  for i = 1:length(typequo)
    if !divisible(S.V.snf[ngens(S.V)+1-i],fmpz((M.p)^typequo[length(typequo)+1-i]))
      return (nmod_mat[])
    end
  end

  #  Dual Module and candidate submodules
  M_dual = dual_module(S)
  candidates = submodules_with_struct(M_dual, typequo)

  #  Dualize the modules
  v = fmpz[divexact(S.V.snf[end],S.V.snf[j]) for j=1:ngens(S.V) ]
  list = (_dualize(x, S.V, v) for x in candidates)

  #  Write the submodules in terms of the given generators
  MatSnf = change_base_ring(R, mS.map)
  return (final_check_and_ans(x, MatSnf, M) for x in list)

end

function final_check_and_ans(x::nmod_mat, MatSnf::nmod_mat, M::ZpnGModule)

  y=x*MatSnf
  @hassert :RayFacElem 1 issubmodule(M, y)
  return y

end

##################################################################################
#
#  Stable Subgroups function
#
##################################################################################

@doc Markdown.doc"""
    stable_subgroups(R::GrpAbFinGen, quotype::Array{Int,1}, act::Array{T, 1}; op=sub)

Given a group $R$, an array of endomorphisms of the group and the type of the quotient, it returns all the stable
subgroups of $R$ such that the corresponding quotient has the required type.
"""
function stable_subgroups(R::GrpAbFinGen, act::Array{T, 1}; op = sub, quotype::Array{Int, 1} = Int[-1], minimal::Bool = false) where T <: Map{GrpAbFinGen, GrpAbFinGen}
  if quotype[1] != -1 && minimal
    error("Cannot compute minimal submodules with prescribed quotient type")
  end
  if quotype[1]!= -1
    #I write quotype as the diagonal entries of the snf form.
    DG = abelian_group(quotype)
    quotype = Int[Int(x) for x in snf(DG)[1].snf]
  end
  c = lcm(quotype[end], exponent(R))
  Q, mQ = quo(R, c, false)
  S, mS = snf(Q)
  #I translate the action to S
  actS = Vector{GrpAbFinGenMap}(undef, length(act))
  for i = 1:length(act)
    imgs = Vector{GrpAbFinGenElem}(undef, ngens(S))
    for j = 1:length(imgs)
      imgs[j] = mS\(act[i](mS(S[j])))
    end
    actS[i] = hom(S, S, imgs, check = false)
  end
  subs_snf = _stable_subgroup_snf(S, actS; quotype = quotype, minimal = minimal)
  #Finally, translate back to R.
  subs = (GrpAbFinGenElem[mQ\mS(x) for x in y] for y in subs_snf)
  return (op(R, x) for x in subs)

end

function _stable_subgroup_snf(R::GrpAbFinGen, act::Array{GrpAbFinGenMap, 1}; quotype::Array{Int,1} = Int[-1], minimal::Bool = false)
  @assert issnf(R)
  c = exponent(R)
  lf = factor(c)
  list = Base.Generator[]
  for p in keys(lf.fac)
    x1 = valuation(c, p)
    G, mG = psylow_subgroup(R, p, false)
    S, mS = snf(G)
    comp = mS*mG

    #We need to distinguish between FqGModule and ZpnGModule (in the first case the algorithm is more efficient)
    if x1 == 1

      F = GF(Int(p), cached = false)
      act_mat = Array{gfp_mat, 1}(undef, length(act))
      for w=1:length(act)
        act_mat[w] = zero_matrix(F, ngens(S), ngens(S))
      end
      for w = 1:ngens(S)
        el = mG(mS(S[w]))
        for z = 1:length(act)
          elz=mS\(haspreimage(mG, act[z](el))[2])
          for l = 1:ngens(S)
            act_mat[z][w,l] = elz[l]
          end
        end
      end
      M = ModAlgAss(act_mat)
      #  Searching for submodules
      if quotype[1] != -1
        ind = 0
        for i = 1:length(quotype)
          if !iscoprime(quotype[i], p)
            ind += 1
          end
        end
        plist = submodules(M, ind)
      else
        if minimal
          plist = minimal_submodules(M)
        else
          plist = submodules(M)
        end
      end
      it = (_lift_and_construct(x, comp) for x in plist)
      push!(list, it)
    else

      RR = ResidueRing(FlintZZ, Int(p)^x1, cached=false)
      act_mat1 = Array{nmod_mat, 1}(undef, length(act))
      for z=1:length(act)
        imgs = GrpAbFinGenElem[]
	      for w = 1:ngens(S)
	        el = act[z](comp(S[w]))
	        fl, el1 = haspreimage(mG, el)
	        @assert fl
          push!(imgs, mS\(el1))
	      end
        act_mat1[z] = map_entries(RR, hom(S, S, imgs).map)
      end

      #  Searching for submodules
      M1 = ZpnGModule(S, act_mat1)
      if quotype[1] != -1
        quotype_p = Int[]
        for i=1:length(quotype)
          v = valuation(quotype[i],p)
          if v>0
            push!(quotype_p, v)
          end
        end
        plist = submodules(M1, typequo = quotype_p)
      else
        if minimal
          plist = minimal_submodules(M1)
        else
          plist = submodules(M1)
        end
      end
      it = (_lift_and_construct(x, comp) for x in plist)
      push!(list, it)
    end
  end
  if isempty(list)
    return Base.Generator([])
  end

  if minimal
    res1 = collect(list[1])
    for j = 2:length(list)
      append!(res1, collect(list[j]))
    end
    return res1
  else
     return ( vcat(cc...) for cc in Iterators.product(list...))
  end
end

function _lift_and_construct(A::Zmodn_mat, mp::GrpAbFinGenMap)
  R = codomain(mp)
  G = domain(mp)
  newsub = GrpAbFinGenElem[]
  for i=1:nrows(A)
    if !iszero_row(A, i)
      y = view(A, i:i, 1:ncols(A))
      el = GrpAbFinGenElem(G, lift(y))
      push!(newsub, mp(el))
    end
  end
  return newsub
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
    to_json(s::SymbolicNode)

Convert a `SymbolicNode` into a JSON string.
"""
function to_json(s::SymbolicNode)
  ref_json = Ref{char_p}(0)
  @mxcall(:MXSymbolSaveToJSON, (MX_handle, Ref{char_p}), s, ref_json)
  return unsafe_string(ref_json[])
end

"""
    from_json(repr :: AbstractString, ::Type{SymbolicNode})

Load a `SymbolicNode` from a JSON string representation.
"""
function from_json(repr :: AbstractString, ::Type{SymbolicNode})
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateFromJSON, (char_p, Ref{MX_handle}), repr, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    load(filename :: AbstractString, ::Type{SymbolicNode})

Load a `SymbolicNode` from a JSON file.
"""
function load(filename :: AbstractString, ::Type{SymbolicNode})
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateFromFile, (char_p, Ref{MX_handle}), filename, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    save(filename :: AbstractString, node :: SymbolicNode)

Save a `SymbolicNode` to a JSON file.
"""
function save(filename :: AbstractString, node :: SymbolicNode)
  @mxcall(:MXSymbolSaveToFile, (MX_handle, char_p), node, filename)
end


text:
using Dates
using Test
using UUIDs

# Avoid colons in the timestamp
const timestamp = Dates.format(now(UTC), dateformat"yyyymmdd-HHMMSS.sss")
const uuid = UUIDs.uuid4()
const folder = "test-$timestamp-$uuid"
println("Using folder \"$folder\" for testing")

include("testsdk.jl")
include("testcli.jl")


text:
function f_real{T}(x::T, N)
    l = 0.1
    a = 0im
    for i=1:N
        a = l*x
    end
end

function f_compl{T}(x::T, N)
    l::Complex128 = 0.1im
    for i::Int=1:N
        a = l::Complex128*x
    end
    #return nothing
end

#N = 1e10

# @time f_real(1.im, 10)
# #@time f_real(1.im, N)
# @time f_real(1., 10)
# #@time f_real(1., N)
# @time f_compl(1., 10)
#@time f_compl(1., N)
@time 1

#const a = Complex128(1,0)
const a = 1.
const N = 10^6

f_compl(a, 1)
@time f_compl(a, N)
#@time f_compl(1.im, N)


text:
module Echidna

include("var_types.jl")
export MOGA_Real, MOGA_Binary, MOGA_Integer

include("core.jl")
export Problem, Algorithm, Solution, Archive, evaluate!,
       insert_solutions!, insert_solution!, compare_pareto_dominance,
       nondominated_cmp, nondominated_sort, crowding_distance,
       nondominated_truncate, reference_point_truncate, copy, copyinto!

include("operators.jl")
export random_candidate, candidate_from_genome, PM, PM!, SBX, SBX!
       tournament_selector

include("algorithms.jl")
export garun, init_pop, iter_generation, NSGAII, NSGAIII

include("problems.jl")
export ZDT1

include("weights.jl")
export generate_regular_reference_points, generate_recursive

include("geom_ops.jl")
export find_extreme_points, find_hyperplane_params, find_axis_intercepts,
       associate_points

include("logging.jl")
export log_population

end # module

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
"onesite tdvp"
@with_kw struct Tdvp <: Algorithm
    tol::Float64 = Defaults.tol
    tolgauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
end

"""
    function timestep(psi, operator, dt, alg,pars = params(psi,operator))

time evolves psi by timestep dt using algorithm alg
"""
function timestep(state::InfiniteMPS, H::Hamiltonian, timestep::Number,alg::Tdvp,pars::Cache=params(state,H))

    newAs = similar(state.AL)

    acjobs = map(enumerate(state.AC)) do (loc,ac)
        @Threads.spawn exponentiate(x->ac_prime(x,loc,state,pars) ,-1im*timestep,ac,Lanczos(tol=alg.tol))
    end
    cjobs = map(enumerate(state.CR)) do (loc,c)
        @Threads.spawn exponentiate(x->c_prime(x,loc,state,pars) ,-1im*timestep,c,Lanczos(tol=alg.tol))
    end

    for loc in 1:length(state)
        (newAcenter,convhist) = fetch(acjobs[loc])
        convhist.converged==0 && @info "time evolving ac($loc) failed $(convhist.normres)"

        (newCenter,convhist) = fetch(cjobs[loc])
        convhist.converged==0 && @info "time evolving c($loc) failed $(convhist.normres)"

        #find Al that best fits these new Acenter and centers
        QAc,_ = leftorth!(newAcenter,alg=TensorKit.QRpos())
        Qc,_ = leftorth!(newCenter,alg=TensorKit.QRpos())
        @tensor Aleft[-1 -2;-3]:=QAc[-1,-2,1]*conj(Qc[-3,1])

        newAs[loc]     = Aleft
    end

    return InfiniteMPS(newAs; tol = alg.tolgauge, maxiter = alg.maxiter,leftgauged = true),pars
end

function timestep(state::Union{FiniteMPS,MPSComoving}, H::Operator, timestep::Number,alg::Tdvp,pars=params(state,H))
    #left to right
    for i in 1:(length(state)-1)
        (state.AC[i],convhist)=let pars = pars,state = state
            exponentiate(x->ac_prime(x,i,state,pars),-1im*timestep/2,state.AC[i],Lanczos(tol=alg.tolgauge))
        end

        (state.CR[i],convhist) = let pars = pars,state = state
            exponentiate(x->c_prime(x,i,state,pars),1im*timestep/2,state.CR[i],Lanczos(tol=alg.tolgauge))
        end
    end

    (state.AC[end],convhist)=let pars = pars,state = state
        exponentiate(x->ac_prime(x,length(state),state,pars),-1im*timestep/2,state.AC[end],Lanczos(tol=alg.tolgauge))
    end

    #right to left
    for i in length(state):-1:2
        (state.AC[i],convhist)= let pars=pars, state = state
            exponentiate(x->ac_prime(x,i,state,pars),-1im*timestep/2,state.AC[i],Lanczos(tol=alg.tolgauge))
        end

        (state.CR[i-1],convhist) = let pars = pars, state = state
            exponentiate(x->c_prime(x,i-1,state,pars),1im*timestep/2,state.CR[i-1],Lanczos(tol=alg.tolgauge))
        end
    end

    (state.AC[1],convhist) = let pars=pars, state = state
        exponentiate(x->ac_prime(x,1,state,pars),-1im*timestep/2,state.AC[1],Lanczos(tol=alg.tolgauge))
    end

    return state,pars
end

"twosite tdvp (works for finite mps's)"
@with_kw struct Tdvp2 <: Algorithm
    tol::Float64 = Defaults.tol
    tolgauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
    trscheme = truncerr(1e-3)
end

#twosite tdvp for finite mps
function timestep(state::Union{FiniteMPS,MPSComoving}, H::Operator, timestep::Number,alg::Tdvp2,pars=params(state,H);rightorthed=false)
    #left to right
    for i in 1:(length(state)-1)
        ac2 = _permute_front(state.AC[i])*_permute_tail(state.AR[i+1])

        (nac2,convhist) = let state=state,pars=pars
            exponentiate(x->ac2_prime(x,i,state,pars),-1im*timestep/2,ac2,Lanczos())
        end

        (nal,nc,nar) = tsvd(nac2,trunc=alg.trscheme)

        state.AC[i] = (nal,complex(nc))
        state.AC[i+1] = (complex(nc),_permute_front(nar))

        if(i!=(length(state)-1))
            (state.AC[i+1],convhist) = let state=state,pars=pars
                exponentiate(x->ac_prime(x,i+1,state,pars),1im*timestep/2,state.AC[i+1],Lanczos())
            end
        end

    end

    #right to left

    for i in length(state):-1:2
        ac2 = _permute_front(state.AC[i-1])*_permute_tail(state.AR[i])

        (nac2,convhist) = let state=state,pars=pars
            exponentiate(x->ac2_prime(x,i-1,state,pars),-1im*timestep/2,ac2,Lanczos())
        end

        (nal,nc,nar) = tsvd(nac2,trunc=alg.trscheme)

        state.AC[i-1] = (nal,complex(nc))
        state.AC[i] = (complex(nc),_permute_front(nar));

        if(i!=2)
            (state.AC[i-1],convhist) = let state=state,pars=pars
                exponentiate(x->ac_prime(x,i-1,state,pars),1im*timestep/2,state.AC[i-1],Lanczos())
            end
        end
    end

    return state,pars
end


text:
using Oceananigans
using OceanTurbulenceParameterEstimation
using Oceananigans.Units
using DataDeps
using GLMakie
using Printf

suite = "six_day_suite"
case_path(case) = @datadep_str("$(suite)_1m/$(case)_instantaneous_statistics.jld2")

cases = [
         "free_convection",
         "weak_wind_strong_cooling",
         "strong_wind_weak_cooling",
         "strong_wind",
         "strong_wind_no_rotation",
        ]

filename = string("catke_simulation_", suite, "_", Δz, "m.jld2")
u_catke_time_series = FieldTimeSeries(filename, "u", boundary_conditions=nothing)
v_catke_time_series = FieldTimeSeries(filename, "v", boundary_conditions=nothing)
T_catke_time_series = FieldTimeSeries(filename, "T", boundary_conditions=nothing)
e_catke_time_series = FieldTimeSeries(filename, "e", boundary_conditions=nothing)

Δz = 4
regrid_size = (1, 1, Int(256/Δz))
observations = [SyntheticObservations(case_path(case); field_names=(:u, :v, :T, :e), regrid_size) for case in cases]

#####
##### Animate!
#####

fig = Figure(resolution=(1800, 800))

T_axs = []
u_axs = []
#e_axs = []

for (i, case) in enumerate(cases)
    ax_l = Label(fig[1, i+1], replace(case, "_" => "\n"), tellwidth=false)

    if i == 1
        ax_T = Axis(fig[2, i+1], ylabel = "z (m)")
        ax_u = Axis(fig[3, i+1], ylabel = "z (m)")
    else
        ax_T = Axis(fig[2, i+1])
        ax_u = Axis(fig[3, i+1])
    end

    #ax_e = Axis(fig[4, i+1])

    #xlims!(ax_e, 0, 0.005)
    xlims!(ax_u, -0.3, 0.5)

    #hidexdecorations!(ax_T)
    #hidexdecorations!(ax_u)
    #hidexdecorations!(ax_e)

    if i > 1
        hideydecorations!(ax_T)
        hideydecorations!(ax_u)
    end

    #hidexdecorations!(ax_e)

    if i > 1
        hidespines!(ax_T, :r, :t, :l)
        hidespines!(ax_u, :r, :t, :l)
    else
        hidespines!(ax_T, :t, :r)
        hidespines!(ax_u, :t, :r)
    end

    #hidespines!(ax_e)

    push!(T_axs, ax_T)
    push!(u_axs, ax_u)
    #push!(e_axs, ax_e)
end

Nt = length(first(observations).times)
slider = Slider(fig[5, :], range=1:Nt, startvalue=1)
n = slider.value

###
### Vertical velocity
###

z = znodes(T_catke_time_series)

labels = ["constant Pr", "variable Pr", "convective adjustment"]

for (j, case) in enumerate(cases)
    obs = observations[j]

    nobs = @lift begin
        model_time = u_catke_time_series.times[$n]
        findfirst(t -> t ≈ model_time, obs.times)
    end

    u_les = @lift interior(obs.field_time_serieses.u[$nobs])[1, 1, :]
    v_les = @lift interior(obs.field_time_serieses.v[$nobs])[1, 1, :]
    T_les = @lift interior(obs.field_time_serieses.T[$nobs])[1, 1, :]
    e_les = @lift interior(obs.field_time_serieses.e[$nobs])[1, 1, :]

    ax_T = T_axs[j]
    ax_u = u_axs[j]
    #ax_e = e_axs[j]

    linewidth = 5
    color = :indigo
    lines!(ax_T, T_les, z; linewidth, color=(color, 0.8))
    lines!(ax_u, u_les, z; linewidth, color=(color, 0.8), label = "u, LES")
    lines!(ax_u, v_les, z; linewidth, color=(color, 0.5), label = "v, LES")
    #lines!(ax_e, e_les, z; linewidth, color=(color, 0.8))

    linewidth = 2
    colors = [:darkorange3, :black, :lightseagreen]
    for i = 1:3
        u_catke = @lift interior(u_catke_time_series[$n])[i, j, :]
        v_catke = @lift interior(v_catke_time_series[$n])[i, j, :]
        T_catke = @lift interior(T_catke_time_series[$n])[i, j, :]
        e_catke = @lift interior(e_catke_time_series[$n])[i, j, :]

        color = colors[i]
        label = labels[i]
        lines!(ax_T, T_catke, z; linewidth, color)
        lines!(ax_u, u_catke, z; linewidth, color, label = "u, $label")
        lines!(ax_u, v_catke, z; linewidth, color=(color, 0.5), label = "v, $label")
        #lines!(ax_e, e_catke, z; linewidth, color)
    end

end

times = u_catke_time_series.times
txt = @lift string("Five boundary layers after ", prettytime(times[$n]))
Label(fig[0, 2:5], txt)

Label(fig[3, 1], "Temperature \n (ᵒC)", tellheight=false)
Label(fig[4, 1], "Velocity \n (m s⁻¹)", tellheight=false)
#Label(fig[4, 1], "Turbulent \n kinetic \n energy", tellheight=false)
Legend(fig[4, 7], u_axs[1], tellheight=false, framevisible=false)

display(fig)

#=
record(fig, string("catke_", suite, ".mp4"), 10:800; framerate=16) do nn
    @info "Drawing frame $nn of $Nt..."
    n[] = nn
end
=#

max_stars_repo_name: glwagner/LocalOceanClosureCalibration

max_issues_repo_name: glwagner/LocalOceanClosureCalibration

max_forks_repo_name: glwagner/LocalOceanClosureCalibration


text:
@testset "Clique percolation" begin
    function setofsets(array_of_arrays)
        Set(map(BitSet, array_of_arrays))
    end

    function test_cliques(graph, expected)
        # Make test results insensitive to ordering
        Set(clique_percolation(graph)) == setofsets(expected)
    end

    g = Graph(5)
    add_edge!(g, 1, 2)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 1)
    add_edge!(g, 1, 4)
    add_edge!(g, 4, 5)
    add_edge!(g, 5, 1)
    @test test_cliques(g, Array[[1, 2, 3], [1, 4, 5]])
end


text:
using GaussianProcesses
using GaussianProcesses: update_mll_and_dmll!, update_cK!, FullCovariancePrecompute
using BenchmarkTools
using DataFrames
using JLD
using Random
using CSV
using LinearAlgebra

kerns = Dict(
    "se" => SEIso(0.3,0.3),
    "mat12" => Mat12Iso(0.3,0.3),
    "rq" => RQIso(0.3,0.3,0.3),
    "se+rq" => SEIso(0.3,0.3) + RQIso(0.3,0.3,0.3),
    "se+mat12" => SEIso(0.3,0.3) + Mat12Iso(0.3,0.3),
    "se*rq" => SEIso(0.3,0.3) * RQIso(0.3,0.3,0.3),
    "se*mat12" => SEIso(0.3,0.3) * Mat12Iso(0.3,0.3),
    "se+mat12+rq" => SEIso(0.3,0.3) + Mat12Iso(0.3,0.3) + RQIso(0.3,0.3,0.3),
    "(se+mat12)*rq" => (SEIso(0.3,0.3) + Mat12Iso(0.3,0.3)) * RQIso(0.3,0.3,0.3),
    # "(se+se2)*rq" => (SEIso(0.3,0.3) + SEIso(0.5,0.5)) * RQIso(0.3,0.3,0.3),
    "mask(se, [1])" => Masked(SEIso(0.3,0.3), [1]),
    "mask(se, [1])+mask(rq, [2:10])" =>  Masked(SEIso(0.3,0.3), [1]) +  Masked(RQIso(0.3,0.3,0.3), collect(2:10)),
    "fix(se, σ)" => fix(SEIso(0.3,0.3), :lσ)
    )

function benchmark_kernel(group, kern)
    XY_df = CSV.read("simdata.csv")
    Xt = XY_df[:,[:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :x10]]
    X = Matrix(transpose(Matrix(Xt)))
    Y = XY_df[:Y]
    @assert length(Y) == nt
    buffer = FullCovariancePrecompute(nt)
    gp = GP(X, Y, MeanConst(0.0), kern, 0.3)
    group["cK"] = @benchmarkable update_cK!($gp)
    group["mll_and_dmll"] = @benchmarkable update_mll_and_dmll!($gp, $buffer)
end

const d = 10        # Input observation dimension
const nt = 3000     # Number of training points
Random.seed!(1)
X = randn(d, nt) # Training data
Y = randn(nt)
XY_df = DataFrame(X')
XY_df[:Y] = Y
XY_df
CSV.write("simdata.csv", XY_df; writeheader=true)

SUITE = BenchmarkGroup()

for (label, k) in kerns
    SUITE[label] = BenchmarkGroup([label])
    benchmark_kernel(SUITE[label], k)
end
;

results = run(SUITE, verbose=false, seconds=10000, samples=20, evals=1)
;

knames = sort(collect(keys(kerns)))
times = [time(results[k]["mll_and_dmll"])/10^6 for k in knames]

df = DataFrame(kernel = knames, times=times)
print(df)

CSV.write("bench_results/GaussianProcesses_jl.csv", df; writeheader=false)


text:
include("../src/Regression.jl")
include("../src/Model.jl")

using BenchmarkTools
using Parameters
using Test
using Plots
using LsqFit

## test functions
test_f(x) = x

hill = Model.hill

## Test LinRange with `step size`
for i in Regression.rangeStepsize(0.2, -13.34, 0.1)
    @show i
end

## Test the time eclapse of `rangeStepsize`
# Result: same allocation
@time test_f.(Regression.rangeStepsize(0, -1300, 1))
@time test_f.(LinRange(0,-1300,1300))

## Test monotonous testing
x_inc = [1,2,3,4,5]
x_dec = [5,3,2,1,0.4]
x_non = [0,1,2,0,2]

@test Regression.is_monotonous(test_f, [4,0]; resolution=20) == true
@test Regression.is_monotonous(x_inc) == true
@test Regression.is_monotonous(x_dec) == true
@test Regression.is_monotonous(x_non) == false

## Compare functor type

# test model
@with_kw struct test_m
    p1
    p2
end
function (set::test_m)(x)
    @unpack p1, p2 = set
    return p1*exp(-x*p2)
end

@. exp_model(x, p) = p[1]*exp(-x*p[2])

# model input
f_ob = Regression.convert_functor_model(test_m) # convert functor to function
@test exp_model(3, [1,2]) == f_ob(3, [1,2])

# Create data
xdata = range(0, stop=10, length=20)
ydata = exp_model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))
p0_bounds = [1.2, 1.2] # we have to start inside the bounds

# Three methods of curve fit
@show fit_functor1 = Regression.curve_fit(f_ob, xdata, ydata,  p0_bounds)
@show fit_tutorial =LsqFit.curve_fit(exp_model, xdata, ydata,  p0_bounds)
@show fit_functor = Regression.curve_fit(test_m, xdata, ydata,  p0_bounds)

@test fit_functor1.param ≈ fit_tutorial.param ≈ fit_functor.param

plot(xdata,ydata, label="Real")
plot!(xdata, test_m(fit_functor1.param...).(xdata), label="Fitted with Regression.jl")
plot!(xdata, exp_model(xdata, fit_tutorial.param), label="Fitted with LsqFit.jl")

## Fit Hill Model
#` Don't use [0] to fit hill function
k, n = 0.5, 6.0
model_real = hill(k, n )

@. H(x, p ) = x^p[2] /(p[1]^p[2] + x^p[2])

xdata = range(0, stop=1, length=20)
ydata = model.(xdata) #+ 0.0001 * rand(length(xdata))

fit = Regression.curve_fit(Model.hill, xdata, ydata,[0.1,0.1])
fit_err = Regression.curve_fit(Model.hill, xdata, ydata,[0.,0.])
fit1 = LsqFit.curve_fit(H, xdata, ydata, [1.,1.])
fit2 = LsqFit.curve_fit(H, xdata, ydata, [1.,1.]; autodiff=:forwarddiff)
fit3 = LsqFit.curve_fit(H, xdata, ydata, [1.,1.]; autodiff=:finiteforward)

#` Plotting
plot(xdata, ydata, label="Noise Data")
plot!(xdata, Model.hill(fit.param...).(xdata) ,label="Fitted Curve")
plot!(xdata, Model.hill(fit_err.param...).(xdata) ,label="Error Curve with 0 init")
plot!(xdata, H(xdata, fit1.param) ,label="Fitted Curve:Default")
plot!(xdata, H(xdata, fit2.param) ,label="Fitted Curve:forwarddiff")
plot!(xdata, H(xdata, fit3.param) ,label="Fitted Curve:finiteforward")

@test k ≈ fit.param[1]
@test n ≈ fit.param[2]


text:
# Copyright 2015-2018 Andreas Noack
module TSVD

    export tsvd

    import Base: *, hcat, maximum, size
    import LinearAlgebra: axpy!, mul!

    using LinearAlgebra
    using LinearAlgebra: BlasComplex, BlasFloat, BlasInt, BlasReal
    using Adapt: adapt

    include("common.jl")
    include("eig.jl")
    include("svd.jl")

end