text:
@doc raw"""
    Centered()

Used in conjunction with `SparseVariationalApproximation`.
States that the `q` field of [`SparseVariationalApproximation`](@ref) is to be interpreted
directly as the approximate posterior over the pseudo-points.

This is also known as the "unwhitened" parametrization [1].

See also [`NonCentered`](@ref).

[1] - https://en.wikipedia.org/wiki/Whitening_transformation
"""
struct Centered end

@doc raw"""
    NonCentered()

Used in conjunction with `SparseVariationalApproximation`.
States that the `q` field of [`SparseVariationalApproximation`](@ref) is to be interpreted
as the approximate posterior over `cholesky(cov(u)).L \ (u - mean(u))`, where `u` are the
pseudo-points.

This is also known as the "whitened" parametrization [1].

See also [`Centered`](@ref).

[1] - https://en.wikipedia.org/wiki/Whitening_transformation
"""
struct NonCentered end

struct SparseVariationalApproximation{Parametrization,Tfz<:FiniteGP,Tq<:AbstractMvNormal}
    fz::Tfz
    q::Tq
end

@doc raw"""
    SparseVariationalApproximation(
        ::Parametrization, fz::FiniteGP, q::AbstractMvNormal
    ) where {Parametrization}

Produce a `SparseVariationalApproximation{Parametrization}`, which packages the prior over
the pseudo-points, `fz`, and the approximate posterior at the pseudo-points, `q`, together
into a single object.

The `Parametrization` determines the precise manner in which `q` and `fz` are interpreted.
Existing parametrizations include [`Centered`](@ref) and [`NonCentered`](@ref).
"""
function SparseVariationalApproximation(
    ::Parametrization, fz::Tfz, q::Tq
) where {Parametrization,Tfz<:FiniteGP,Tq<:AbstractMvNormal}
    return SparseVariationalApproximation{Parametrization,Tfz,Tq}(fz, q)
end

"""
    SparseVariationalApproximation(fz::FiniteGP, q::AbstractMvNormal)

Packages the prior over the pseudo-points `fz`, and the approximate posterior at the
pseudo-points, which is `mean(fz) + cholesky(cov(fz)).L * ε`, `ε ∼ q`.

Shorthand for
```julia
SparseVariationalApproximation(NonCentered(), fz, q)
```
"""
function SparseVariationalApproximation(fz::FiniteGP, q::AbstractMvNormal)
    return SparseVariationalApproximation(NonCentered(), fz, q)
end

@doc raw"""
    posterior(sva::SparseVariationalApproximation{Centered})

Compute the approximate posterior [1] over the process `f =
sva.fz.f`, given inducing inputs `z = sva.fz.x` and a variational
distribution over inducing points `sva.q` (which represents ``q(u)``
where `u = f(z)`). The approximate posterior at test points ``x^*``
where ``f^* = f(x^*)`` is then given by:

```math
q(f^*) = \int p(f | u) q(u) du
```
which can be found in closed form.

[1] - Hensman, James, Alexander Matthews, and Zoubin Ghahramani. "Scalable
variational Gaussian process classification." Artificial Intelligence and
Statistics. PMLR, 2015.
"""
function AbstractGPs.posterior(sva::SparseVariationalApproximation{Centered})
    # m* = K*u Kuu⁻¹ (mean(q) - mean(fz))
    #    = K*u α
    # Centered: α = Kuu⁻¹ (m - mean(fz))
    # [NonCentered: α = Lk⁻ᵀ m]
    # V** = K** - K*u (Kuu⁻¹ - Kuu⁻¹ cov(q) Kuu⁻¹) Ku*
    #     = K** - K*u (Kuu⁻¹ - Kuu⁻¹ cov(q) Kuu⁻¹) Ku*
    #     = K** - (K*u Lk⁻ᵀ) (Lk⁻¹ Ku*) + (K*u Lk⁻ᵀ) Lk⁻¹ cov(q) Lk⁻ᵀ (Lk⁻¹ Ku*)
    #     = K** - A'A + A' Lk⁻¹ cov(q) Lk⁻ᵀ A
    #     = K** - A'A + A' Lk⁻¹ Lq Lqᵀ Lk⁻ᵀ A
    #     = K** - A'A + A' B B' A
    # A = Lk⁻¹ Ku*
    # Centered: B = Lk⁻¹ Lq
    # [NonCentered: B = Lq]
    q, fz = sva.q, sva.fz
    m, S = mean(q), _chol_cov(q)
    Kuu = _chol_cov(fz)
    B = chol_lower(Kuu) \ chol_lower(S)
    α = Kuu \ (m - mean(fz))
    data = (Kuu=Kuu, B=B, α=α)
    return ApproxPosteriorGP(sva, fz.f, data)
end

#
# NonCentered Parametrization.
#

@doc raw"""
    posterior(sva::SparseVariationalApproximation{NonCentered})

Compute the approximate posterior [1] over the process `f =
sva.fz.f`, given inducing inputs `z = sva.fz.x` and a variational
distribution over inducing points `sva.q` (which represents ``q(ε)``
where `ε = cholesky(cov(fz)).L \ (f(z) - mean(f(z)))`). The approximate posterior at test
points ``x^*`` where ``f^* = f(x^*)`` is then given by:

```math
q(f^*) = \int p(f | ε) q(ε) du
```
which can be found in closed form.

[1] - Hensman, James, Alexander Matthews, and Zoubin Ghahramani. "Scalable
variational Gaussian process classification." Artificial Intelligence and
Statistics. PMLR, 2015.
"""
function AbstractGPs.posterior(sva::SparseVariationalApproximation{NonCentered})
    # u = Lk v + mean(fz), v ~ q
    # m* = K*u Kuu⁻¹ Lk (mean(u) - mean(fz))
    #    = K*u (Lk Lkᵀ)⁻¹ Lk mean(q)
    #    = K*u Lk⁻ᵀ Lk⁻¹ Lk mean(q)
    #    = K*u Lk⁻ᵀ mean(q)
    #    = K*u α
    # NonCentered: α = Lk⁻ᵀ m
    # [Centered: α = Kuu⁻¹ (m - mean(fz))]
    # V** = K** - K*u (Kuu⁻¹ - Kuu⁻¹ Lk cov(q) Lkᵀ Kuu⁻¹) Ku*
    #     = K** - K*u (Kuu⁻¹ - (Lk Lkᵀ)⁻¹ Lk cov(q) Lkᵀ (Lk Lkᵀ)⁻¹) Ku*
    #     = K** - K*u (Kuu⁻¹ - Lk⁻ᵀ Lk⁻¹ Lk cov(q) Lkᵀ Lk⁻ᵀ Lk⁻¹) Ku*
    #     = K** - K*u (Kuu⁻¹ - Lk⁻ᵀ cov(q) Lk⁻¹) Ku*
    #     = K** - (K*u Lk⁻ᵀ) (Lk⁻¹ Ku*) - (K*u Lk⁻ᵀ) Lq Lqᵀ (Lk⁻¹ Ku*)
    #     = K** - A'A - (K*u Lk⁻ᵀ) Lq Lqᵀ (Lk⁻¹ Ku*)
    #     = K** - A'A - A' B B' A
    # A = Lk⁻¹ Ku*
    # NonCentered: B = Lq
    # [Centered: B = Lk⁻¹ Lq]
    q, fz = sva.q, sva.fz
    m = mean(q)
    Kuu = _chol_cov(fz)
    α = chol_lower(Kuu)' \ m
    Sv = _chol_cov(q)
    B = chol_lower(Sv)
    data = (Kuu=Kuu, B=B, α=α)
    return ApproxPosteriorGP(sva, fz.f, data)
end

function AbstractGPs.posterior(
    sva::SparseVariationalApproximation, fx::FiniteGP, ::AbstractVector{<:Real}
)
    @assert sva.fz.f === fx.f
    return posterior(sva)
end

#
# Various methods implementing the Internal AbstractGPs API.
# See AbstractGPs.jl API docs for more info.
#

function Statistics.mean(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation}, x::AbstractVector
)
    return mean(f.prior, x) + cov(f.prior, x, inducing_points(f)) * f.data.α
end

# A = Lk⁻¹ Ku* is the projection matrix used in computing the predictive variance of the SparseVariationalApproximation posterior.
function _A_and_Kuf(f, x)
    Kuf = cov(f.prior, inducing_points(f), x)
    A = chol_lower(f.data.Kuu) \ Kuf
    return A, Kuf
end

_A(f, x) = first(_A_and_Kuf(f, x))

function Statistics.cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation}, x::AbstractVector
)
    A = _A(f, x)
    return cov(f.prior, x) - At_A(A) + At_A(f.data.B' * A)
end

function Statistics.var(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation}, x::AbstractVector
)
    A = _A(f, x)
    return var(f.prior, x) - diag_At_A(A) + diag_At_A(f.data.B' * A)
end

function StatsBase.mean_and_cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation}, x::AbstractVector
)
    A, Kuf = _A_and_Kuf(f, x)
    μ = mean(f.prior, x) + Kuf' * f.data.α
    Σ = cov(f.prior, x) - At_A(A) + At_A(f.data.B' * A)
    return μ, Σ
end

function StatsBase.mean_and_var(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation}, x::AbstractVector
)
    A, Kuf = _A_and_Kuf(f, x)
    μ = mean(f.prior, x) + Kuf' * f.data.α
    Σ_diag = var(f.prior, x) - diag_At_A(A) + diag_At_A(f.data.B' * A)
    return μ, Σ_diag
end

function Statistics.cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation},
    x::AbstractVector,
    y::AbstractVector,
)
    B = f.data.B
    Ax = _A(f, x)
    Ay = _A(f, y)
    return cov(f.prior, x, y) - Ax'Ay + Ax' * B * B' * Ay
end

#
# Misc utility.
#

inducing_points(f::ApproxPosteriorGP{<:SparseVariationalApproximation}) = f.approx.fz.x

#
# elbo
#

"""
    elbo(
        sva::SparseVariationalApproximation,
        fx::FiniteGP,
        y::AbstractVector{<:Real};
        num_data=length(y),
        quadrature=DefaultQuadrature(),
    )

Compute the Evidence Lower BOund from [1] for the process `f = fx.f ==
svgp.fz.f` where `y` are observations of `fx`, pseudo-inputs are given by `z =
svgp.fz.x` and `q(u)` is a variational distribution over inducing points `u =
f(z)`.

`quadrature` selects which method is used to calculate the expected loglikelihood in
the ELBO. The options are: `DefaultQuadrature()`, `Analytic()`, `GaussHermite()` and
`MonteCarlo()`. For likelihoods with an analytic solution, `DefaultQuadrature()` uses this
exact solution. If there is no such solution, `DefaultQuadrature()` either uses
`GaussHermite()` or `MonteCarlo()`, depending on the likelihood.

N.B. the likelihood is assumed to be Gaussian with observation noise `fx.Σy`.
Further, `fx.Σy` must be isotropic - i.e. `fx.Σy = α * I`.

[1] - Hensman, James, Alexander Matthews, and Zoubin Ghahramani. "Scalable
variational Gaussian process classification." Artificial Intelligence and
Statistics. PMLR, 2015.
"""
function AbstractGPs.elbo(
    sva::SparseVariationalApproximation,
    fx::FiniteGP{<:AbstractGP,<:AbstractVector,<:Diagonal{<:Real,<:Fill}},
    y::AbstractVector{<:Real};
    num_data=length(y),
    quadrature=DefaultQuadrature(),
)
    @assert sva.fz.f === fx.f
    return _elbo(quadrature, sva, fx, y, GaussianLikelihood(fx.Σy[1]), num_data)
end

function AbstractGPs.elbo(
    ::SparseVariationalApproximation, ::FiniteGP, ::AbstractVector; kwargs...
)
    return error(
        "The observation noise fx.Σy must be homoscedastic.\n",
        "To avoid this error, construct fx using: f = GP(kernel); fx = f(x, σ²)",
        ", where σ² is a positive Real.",
    )
end

"""
    elbo(
        sva::SparseVariationalApproximation,
        lfx::LatentFiniteGP,
        y::AbstractVector;
        num_data=length(y),
        quadrature=DefaultQuadrature(),
    )

Compute the ELBO for a LatentGP with a possibly non-conjugate likelihood.
"""
function AbstractGPs.elbo(
    sva::SparseVariationalApproximation,
    lfx::LatentFiniteGP,
    y::AbstractVector;
    num_data=length(y),
    quadrature=DefaultQuadrature(),
)
    @assert sva.fz.f === lfx.fx.f
    return _elbo(quadrature, sva, lfx.fx, y, lfx.lik, num_data)
end

# Compute the common elements of the ELBO
function _elbo(
    quadrature::QuadratureMethod,
    sva::SparseVariationalApproximation,
    fx::FiniteGP,
    y::AbstractVector,
    lik,
    num_data::Integer,
)
    @assert sva.fz.f === fx.f

    f_post = posterior(sva)
    q_f = marginals(f_post(fx.x))
    variational_exp = expected_loglik(quadrature, y, q_f, lik)

    n_batch = length(y)
    scale = num_data / n_batch
    return sum(variational_exp) * scale - _prior_kl(sva)
end

_prior_kl(sva::SparseVariationalApproximation{Centered}) = kldivergence(sva.q, sva.fz)

function _prior_kl(sva::SparseVariationalApproximation{NonCentered})
    m_ε = mean(sva.q)
    C_ε = _cov(sva.q)

    # trace_term = tr(C_ε)  # does not work due to PDMat / Zygote issues
    L = chol_lower(_chol_cov(sva.q))
    trace_term = sum(L .^ 2)  # TODO remove AD workaround

    return (trace_term + m_ε'm_ε - length(m_ε) - logdet(C_ε)) / 2
end


text:
##############################################################################
##
## FixedEffect
##
## The categoricalarray may have pools that are never referred. Note that the pool does not appear in FixedEffect anyway.
##
##############################################################################

struct FixedEffect{R <: AbstractVector{<:Integer}, I <: AbstractVector{<: Real}}
    refs::R                 # refs of the original CategoricalVector
    interaction::I          # the continuous interaction
    n::Int                  # Number of potential values (= maximum(refs))
    function FixedEffect{R, I}(refs, interaction, n) where {R <: AbstractVector{<:Integer}, I <: AbstractVector{<: Real}}
        isa(refs, Vector) && maximum(refs) > n && error("Categorical Vector used to construct Fixed Effect is malformed. Some elements of refs do not refer to pool")
        new(refs, interaction, n)
    end
end

function FixedEffect(args...; interaction::AbstractVector = Ones{Float64}(length(args[1])))
    if length(args) != 1
        FixedEffect(group(args...); interaction = interaction)
    else
        x = args[1]
        if !isa(x, CategoricalVector)
            FixedEffect(categorical(x); interaction = interaction)
        else
            FixedEffect{typeof(x.refs), typeof(interaction)}(x.refs, interaction, length(x.pool))
        end
    end
end

function Base.show(io::IO, fe::FixedEffect)
    println(io, "Refs:        ", length(fe.refs), "-element ", typeof(fe.refs))
    println(io, "             ", Int.(fe.refs[1:min(5, length(fe.refs))]), "...")
    println(io, "Interaction: ", length(fe.interaction), "-element ", typeof(fe.interaction))
    println(io, "             ", fe.interaction[1:min(5, length(fe.interaction))], "...")
end

Base.ismissing(fe::FixedEffect) = any(fe.refs .== 0) | ismissing(fe.interaction)
Base.length(fe::FixedEffect) = length(fe.refs)
Base.eltype(fe::FixedEffect) = eltype(I)

##############################################################################
##
## group transform multiple CategoricalVector into one
## Output is a CategoricalVector where pool is type Int64, equal to ranking of group
## Missing in some row mean result has Missing on this row
##
##############################################################################

function group(args...)
    v = categorical(args[1])
    if length(args) == 1
        refs = v.refs
    else
        refs, ngroups = convert(Vector{UInt}, v.refs), length(levels(v))
        for j = 2:length(args)
            v = categorical(args[j])
            refs, ngroups = pool_combine!(refs, v, ngroups)
        end
    end
    factorize!(refs)
end
function pool_combine!(refs::Array{T, N}, dv::CategoricalVector, ngroups::Integer) where {T, N}
    for i in 1:length(refs)
        # if previous one is NA or this one is NA, set to NA
        refs[i] = (dv.refs[i] == 0 || refs[i] == zero(T)) ? zero(T) : refs[i] + (dv.refs[i] - 1) * ngroups
    end
    return refs, ngroups * length(levels(dv))
end

function factorize!(refs::Vector{T}) where {T}
    uu = sort!(unique(refs))
    has_missing = uu[1] == 0
    dict = Dict{T, Int}(zip(uu, (1-has_missing):(length(uu)-has_missing)))
    newrefs = zeros(UInt32, length(refs))
    for i in 1:length(refs)
         newrefs[i] = dict[refs[i]]
    end
    Tout = has_missing ? Union{Int, Missing} : Int
    CategoricalArray{Tout, 1}(newrefs, CategoricalPool(collect(1:(length(uu)-has_missing))))
end

##############################################################################
##
## Find connected components
##
##############################################################################
# Return a vector of sets that contains the indices of each unique value
function refsrev(fe::FixedEffect)
    out = Vector{Int}[Int[] for _ in 1:fe.n]
    for i in eachindex(fe.refs)
        push!(out[fe.refs[i]], i)
    end
    return out
 end

# Returns a vector of all components
# A component is a vector that, for each fixed effect,
# contains all the refs that are included in the component.
function components(fes::AbstractVector{<:FixedEffect})
    refs_vec = Vector{UInt32}[fe.refs for fe in fes]
    refsrev_vec = Vector{Vector{Int}}[refsrev(fe) for fe in fes]
    visited = falses(length(refs_vec[1]))
    out = Vector{Set{Int}}[]
    for i in eachindex(visited)
        if !visited[i]
            # obs not visited yet, so create new component
            component_vec = Set{UInt32}[Set{UInt32}() for _ in 1:length(refsrev_vec)]
            # visit all obs in the same components
            tovisit = Set{Int}(i)
            while !isempty(tovisit)
                for (component, refs, refsrev) in zip(component_vec, refs_vec, refsrev_vec)
                    ref = refs[i]
                    # if group is not in component yet
                    if ref ∉ component
                        # add group to the component
                        push!(component, ref)
                        # visit other observations in same group
                        union!(tovisit, refsrev[ref])
                    end
                end
                # mark obs as visited
                i = pop!(tovisit)
                visited[i] = true
            end
            push!(out, component_vec)
        end
    end
    return out
end

##############################################################################
##
## normalize! a vector of fixedeffect coefficients using connected components
##
##############################################################################

function normalize!(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect}; kwargs...)
    # The solution is generally not unique. Find connected components and scale accordingly
    idx = findall(fe -> isa(fe.interaction, Ones), fes)
    length(idx) >= 2 && rescale!(view(fecoefs, idx), view(fes, idx))
    return fecoefs
end

function rescale!(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect})
    for component_vec in components(fes)
        m = 0.0
        # demean all fixed effects except the first
        for j in length(fecoefs):(-1):2
            fecoef, component = fecoefs[j], component_vec[j]
            mj = 0.0
            for k in component
                mj += fecoef[k]
            end
            mj = mj / length(component)
            for k in component
                fecoef[k] -= mj
            end
            m += mj
        end
        # rescale the first fixed effects
        fecoef, component = fecoefs[1], component_vec[1]
        for k in component
            fecoef[k] += m
        end
    end
end

function full(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect})
    [fecoef[fe.refs] for (fecoef, fe) in zip(fecoefs, fes)]
end

ext: jl


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems.Serialization
using Test
using SparseArrays

@bench @testset "Serialization" begin
    file = "sparse_matrix.txt"

    Random.seed!(0)
    A = sprand(40, 20, 0.05)
    write_sparse_matrix(file, A)

    B = read_sparse_matrix(file)

    rm(file)

    @test A == B
end


text:
using BlockDiagonalMatrices
using SparseArrays
using LinearAlgebra
using TimerOutputs

struct LocalDephasingModel <: Model
    params::ModelParameters
    Jx::BlockDiagonal
    Jy::BlockDiagonal
    Jz::BlockDiagonal
    Jx2::BlockDiagonal
    Jy2::BlockDiagonal
    Jz2::BlockDiagonal
    inefficient_measurement::BlockDiagonal
    second_term::SuperOperator
    M0::BlockDiagonal
    dM::BlockDiagonal
    measurement::Array{Eigen}

    function LocalDephasingModel(modelparams::ModelParameters, liouvillianfile::Union{String,Nothing}=nothing)
        Nj = modelparams.Nj
        dt = modelparams.dt
        Gamma = modelparams.Gamma
        kind = modelparams.kind
        ω = modelparams.omega
        η = modelparams.eta

        # Spin operators
        (Jx, Jy, Jz) = map(blockdiagonal, jspin(Nj))

        # TODO: Find better name
        inefficient_measurement = sqrt((1 - η) * dt * Gamma) * Jy

        second_term = dt * (kind / 2) * (isnothing(liouvillianfile) ?
                                            initliouvillian(Nj) :
                                            initliouvillian(Nj, liouvillianfile))

        dropzeros!(second_term)
        second_term = SuperOperator(second_term)

        Jx2 = Jx^2
        Jy2 = Jy^2
        Jz2 = Jz^2
        H = ω * Jz
        dH = Jz

        # Kraus-like operator, trajectory-independent part
        M0 = (I - 1im * H * dt -
              0.25 * dt * kind * Nj * I - # The Id comes from the squares of sigmaz_j
              (Gamma / 2) * Jy2 * dt)

        # Derivative of the Kraus-like operator wrt to ω
        dM = -1im * dH * dt

        measurement = map(x -> eigen(Matrix(x)), Jy.blocks)

        new(modelparams, Jx, Jy, Jz, Jx2, Jy2, Jz2, inefficient_measurement, second_term, M0, dM, measurement)
    end
end

    get_time(m::Model) = get_time(m.params)

function initliouvillian(Nj::Integer)
    sys = piqs.Dicke(Nj)
    sys.dephasing = 4.

    liouvillian = tosparse(sys.liouvillian())
    return liouvillian + Nj * I
end

function initliouvillian(Nj::Integer, filename::String)
    liouvillian = sparse_fromfile(filename)
    return liouvillian + Nj * I
end

function measure_current(state::BlockDiagonalState, model::LocalDephasingModel)
    @inline dW() = sqrt(model.params.dt) * randn() # Define the Wiener increment
    # Homodyne current (Eq. 35)
    # dy = 2 sqrt(Gamma * eta) * tr(ρ * Jy) * dt + dW
    mul!(state._tmp1, model.Jy, state.ρ)
    return 2 * sqrt(model.params.Gamma * model.params.eta) * real(tr(state._tmp1)) * model.params.dt + dW()
end

function get_kraus(model::LocalDephasingModel, dy::Real)
    p = model.params
    kraus = similar(model.M0)
    # Kraus operator Eq. (36)
    @inbounds for i in eachindex(kraus.blocks)
        kraus.blocks[i] = model.M0.blocks[i] + sqrt(p.eta * p.Gamma) * model.Jy.blocks[i] * dy +
        p.eta * (p.Gamma / 2) * model.Jy2.blocks[i] * (dy^2 - p.dt)
    end
    return kraus
end

function updatestate!(state::BlockDiagonalState, model::LocalDephasingModel, dy::Real)
    # Get Kraus operator for measured current dy
    @timeit_debug "get_kraus" M = get_kraus(model, dy)

    # Non-allocating code for
    # new_ρ = Mpre * Mpost * ρ + second_term * ρ

    @timeit_debug "rho" begin
        mul!(state._tmp1, state.ρ, M')
        mul!(state._new_ρ, M, state._tmp1)

            @timeit_debug "superop" apply_superop!(state._tmp1, model.second_term, state.ρ)

        if model.params.eta < 1.0
            mul!(state._tmp2, model.inefficient_measurement, state.ρ)
            mul!(state._tmp1, state._tmp2, model.inefficient_measurement', 1.0, 1.0)
        end

        # TODO: Replace with broadcasting once implemented
        @inbounds for i in eachindex(state._new_ρ.blocks)
            state._new_ρ.blocks[i] .+= state._tmp1.blocks[i]
        end

        zchop!(state._new_ρ) # Round off elements smaller than 1e-14
        tr_ρ = tr(state._new_ρ)
        # Evolve the unnormalized derivative wrt ω
    end
    # Non-allocating code for
    # τ = (Mpre * (Mpost * τ  +  dMpost * ρ) + dMpre * Mpost * ρ +
    #      second_term * τ )/ tr_ρ;
    @timeit_debug "tau" begin
        mul!(state._tmp1, state.ρ, model.dM')
        mul!(state._tmp1, state.τ, M', 1., 1.)

        @timeit_debug "superop" apply_superop!(state._tmp2, model.second_term, state.τ)
        mul!(state._tmp2, M, state._tmp1, 1., 1.)

        if model.params.eta < 1.0
            mul!(state._tmp1, model.inefficient_measurement, state.τ)
            mul!(state._tmp2, state._tmp1, model.inefficient_measurement', 1.0, 1.0)
        end

        mul!(state._tmp1, state.ρ, M')
        mul!(state._tmp2, model.dM, state._tmp1, 1., 1.)

        # TODO: Use broadcasting when it is implemented
        @inbounds for i in eachindex(state.τ.blocks)
            state.τ.blocks[i] .= state._tmp2.blocks[i] ./ tr_ρ
        end
    end
    zchop!(state.τ) # Round off elements smaller than 1e-14

    tr_τ = tr(state.τ)
    # Now we can renormalize ρ and its derivative wrt ω
    # TODO: Use broadcasting when it is implemented
    @timeit_debug "norms" begin
        @inbounds for i = eachindex(state.ρ.blocks)
            state.ρ.blocks[i] .= state._new_ρ.blocks[i] ./ tr_ρ
        end
        @inbounds for i = eachindex(state.dρ.blocks)
            state.dρ.blocks[i] .= state.τ.blocks[i] .- tr_τ .* state.ρ.blocks[i]
        end
    end
    return tr_ρ, tr_τ
end

expectation_value!(state::State, op::AbstractArray) = real(tr(mul!(state._tmp1, op, density_matrix(state))))

coherentspinstate(model::LocalDephasingModel) = blockdiag_css(model.params.Nj)


text:
module AxPlot

using PyPlot
using PyCall
using LinearAlgebra
using Flux: Tracker
using Requires

include("main.jl")
include("utils.jl")

function __init__()
    @require InferGMM="e4654f38-ff99-11e8-3938-77f3ff3d0ead" include("gmm.jl")
end

end

max_stars_repo_name: ornithos/AxPlot

max_issues_repo_name: ornithos/AxPlot

max_forks_repo_name: ornithos/AxPlot


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct LCFilter <: Filter
        lf::Float64
        rf::Float64
        cf::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of a LCL filter outside the converter

# Arguments
- `lf::Float64`: filter inductance, validation range: `(0, nothing)`
- `rf::Float64`: filter resistance, validation range: `(0, nothing)`
- `cf::Float64`: filter capacitance, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the LCFilter model are:
	ir_filter: Real current out of the filter,
	ii_filter: Imaginary current out of the filter
- `n_states::Int`: LCFilter has two states
"""
mutable struct LCFilter <: Filter
    "filter inductance"
    lf::Float64
    "filter resistance"
    rf::Float64
    "filter capacitance"
    cf::Float64
    ext::Dict{String, Any}
    "The states of the LCFilter model are:
	ir_filter: Real current out of the filter,
	ii_filter: Imaginary current out of the filter"
    states::Vector{Symbol}
    "LCFilter has two states"
    n_states::Int
end

function LCFilter(lf, rf, cf, ext=Dict{String, Any}(), )
    LCFilter(lf, rf, cf, ext, [:ir_filter, :ii_filter], 2, )
end

function LCFilter(; lf, rf, cf, ext=Dict{String, Any}(), states=[:ir_filter, :ii_filter], n_states=2, )
    LCFilter(lf, rf, cf, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function LCFilter(::Nothing)
    LCFilter(;
        lf=0,
        rf=0,
        cf=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`LCFilter`](@ref) `lf`."""
get_lf(value::LCFilter) = value.lf
"""Get [`LCFilter`](@ref) `rf`."""
get_rf(value::LCFilter) = value.rf
"""Get [`LCFilter`](@ref) `cf`."""
get_cf(value::LCFilter) = value.cf
"""Get [`LCFilter`](@ref) `ext`."""
get_ext(value::LCFilter) = value.ext
"""Get [`LCFilter`](@ref) `states`."""
get_states(value::LCFilter) = value.states
"""Get [`LCFilter`](@ref) `n_states`."""
get_n_states(value::LCFilter) = value.n_states

"""Set [`LCFilter`](@ref) `lf`."""
set_lf!(value::LCFilter, val) = value.lf = val
"""Set [`LCFilter`](@ref) `rf`."""
set_rf!(value::LCFilter, val) = value.rf = val
"""Set [`LCFilter`](@ref) `cf`."""
set_cf!(value::LCFilter, val) = value.cf = val
"""Set [`LCFilter`](@ref) `ext`."""
set_ext!(value::LCFilter, val) = value.ext = val


text:
"""
Provides the basic functionality to interact with the Collatz conjecture. The
parameterisation uses the same (P,a,b) notation as Conway's generalisations.
Besides the function and reverse function, there is also functionality to retrieve
the hailstone sequence, the "stopping time"/"total stopping time", or tree-graph.

# Examples
```
julia> import Pkg; Pkg.add("Collatz"); using Collatz
```
"""
module Collatz
#TODO: ^ add jldoctest back to the above example when this is in general

export _ErrMsg, _CC, _KNOWN_CYCLES
export collatz_function, reverse_collatz_function, hailstone_sequence, stopping_time, tree_graph

"The four known cycles (besides 0->0), for the default parameterisation."
const _KNOWN_CYCLES = [[1, 4, 2], [-1, -2], [-5, -14, -7, -20, -10], [-17, -50, -25, -74, -37, -110, -55, -164, -82, -41, -122, -61, -182, -91, -272, -136, -68, -34]]
"The value up to which has been proven numerically, for the default parameterisation."
const __VERIFIED_MAXIMUM = 295147905179352825856
"The value down to which has been proven numerically, for the default parameterisation."
const __VERIFIED_MINIMUM = -272  #&TODO: Check the actual lowest bound.

"""
Error message constant.
"""
module _ErrMsg
const SANE_PARAMS_P = "'P' should not be 0 ~ violates modulo being non-zero."
const SANE_PARAMS_A = "'a' should not be 0 ~ violates the reversability."
end # module _ErrMsg
import ._ErrMsg

"""
Cycle Control: Descriptive flags to indicate when some event occurs in the
hailstone sequences, when set to verbose, or stopping time check. A module used
to wrap an enum to reference the values as `_CC.ABC` rather than `ABC::CC`

# Examples
```jldoctest
julia> string(_CC.STOPPING_TIME)
"STOPPING_TIME"
julia> string(_CC.TOTAL_STOPPING_TIME)
"TOTAL_STOPPING_TIME"
julia> string(_CC.CYCLE_INIT)
"CYCLE_INIT"
julia> string(_CC.CYCLE_LENGTH)
"CYCLE_LENGTH"
julia> string(_CC.MAX_STOP_OOB)
"MAX_STOP_OOB"
julia> string(_CC.ZERO_STOP)
"ZERO_STOP"
```
"""
module _CC
# The elements of an enum are accessed as string(ABC::CC) or just string(ABC) so
# wrap the enum in a module so they can be referenced as string(_CC.ABC)
@enum CC begin
    STOPPING_TIME
    TOTAL_STOPPING_TIME
    CYCLE_INIT
    CYCLE_LENGTH
    MAX_STOP_OOB  # ~ "out of bounds"
    ZERO_STOP
end # @enum CC
end # module _CC
import ._CC

"""
    __assert_sane_parameterisation(P, a, b)

Handles the sanity check for the parameterisation (P,a,b) required by both
the function and reverse function.

# Args
- `P::Integer`: Modulus used to devide n, iff n is equivalent to (0 mod P).
- `a::Integer`: Factor by which to multiply n.
- `b::Integer`: Value to add to the scaled value of n.
"""
function __assert_sane_parameterisation(P::Integer, a::Integer, b::Integer)
    # Sanity check (P,a,b) ~ P absolutely can't be 0. a "could" be zero
    # theoretically, although would violate the reversability (if ~a is 0 then a
    # value of "b" as the input to the reverse function would have a pre-emptive
    # value of every number not divisible by P). The function doesn't _have_ to
    # be reversable, but we are only interested in dealing with the class of
    # functions that exhibit behaviour consistant with the collatz function. If
    # _every_ input not divisable by P went straight to "b", it would simply
    # cause a cycle consisting of "b" and every b/P^z that is an integer. While
    # P in [-1, 1] could also be a reasonable check, as it makes every value
    # either a 1 or 2 length cycle, it's not strictly an illegal operation.
    # "b" being zero would cause behaviour not consistant with the collatz
    # function, but would not violate the reversability, so no check either.
    # " != 0" is redundant for python assertions.
    @assert P != zero(P) _ErrMsg.SANE_PARAMS_P
    @assert a != zero(a) _ErrMsg.SANE_PARAMS_A
end

"""
    collatz_function(n; P=2, a=3, b=1)

Returns the output of a single application of a Collatz-esque function.

# Args
- `n::Integer`: The value on which to perform the Collatz-esque function

# Kwargs
- `P::Integer=2`: Modulus used to devide n, iff n is equivalent to (0 mod P).
- `a::Integer=3`: Factor by which to multiply n.
- `b::Integer=1`: Value to add to the scaled value of n.

# Examples
```jldoctest
julia> collatz_function(5)
16
```
```jldoctest
julia> collatz_function(14, P=7)
2
```
```jldoctest
julia> collatz_function(15, P=7, a=5, b=3)
78
```
"""
function collatz_function(n::Integer; P::Integer=2, a::Integer=3, b::Integer=1)
    __assert_sane_parameterisation(P,a,b)
    if n%P == 0
        # In Julia, "Integer//Integer" creates a "Rational{Int64}" type, rather
        # than strict integer division (which is done instead with "÷"). See;
        # https://docs.julialang.org/en/v1/manual/mathematical-operations/
        # https://docs.julialang.org/en/v1/manual/complex-and-rational-numbers/
        # So we can either "n÷P", or "numerator(n//P)"
        return n÷P
    else
        return (a*n+b)
    end
end

"""
    reverse_collatz_function(n; P=2, a=3, b=1)

Returns the output of a single application of a Collatz-esque reverse function.

# Args
- `n::Integer`: The value on which to perform the reverse Collatz function

# Kwargs
- `P::Integer=2`: Modulus used to devide n, iff n is equivalent to (0 mod P).
- `a::Integer=3`: Factor by which to multiply n.
- `b::Integer=1`: Value to add to the scaled value of n.

# Examples
```jldoctest
julia> print(reverse_collatz_function(1))
[2]
```
```jldoctest
julia> print(reverse_collatz_function(4))
[1, 8]
```
```jldoctest
julia> print(reverse_collatz_function(3, P=-3, a=-2, b=-5))
[-9, -4]
```
"""
function reverse_collatz_function(n::Integer; P::Integer=2, a::Integer=3, b::Integer=1)
    __assert_sane_parameterisation(P,a,b)
    # Every input can be reversed as the result of "n/P" division, which yields
    # "Pn"... {f(n) = an + b}≡{(f(n) - b)/a = n} ~ if n was such that the
    # muliplication step was taken instead of the division by the modulus, then
    # (f(n) - b)/a) must be an integer that is not in (0 mod P). Because we're
    # not placing restrictions on the parameters yet, although there is a better
    # way of shortcutting this for the default variables, we need to always
    # attempt (f(n) - b)/a)
    pre_values = [P*n]
    if (n-b)%a == 0 && (n-b)%(P*a) != 0
        push!(pre_values, ((n-b)÷a))
        sort!(pre_values)
    end
    return pre_values
end

"""
    __initial_value_outside_verified_range(x)

Checks if the initial value is greater than __VERIFIED_MAXIMUM or less than
__VERIFIED_MINIMUM. Only intended for the default parameterisation.

# Args
- `x::Integer`: The initial value to check if it is within range or not.
"""
function __initial_value_outside_verified_range(x::Integer)
    return (__VERIFIED_MAXIMUM < x) || (x < __VERIFIED_MINIMUM)
end

"""
    __stopping_time_terminus(n, total_stop)

Provides the appropriate lambda to use to check if iterations on an initial
value have reached either the stopping time, or total stopping time.

# Args
- `n::Integer`: The initial value to confirm against a stopping time check.
- `total_stop::Bool`: If false, the lambda will confirm that iterations
    of n have reached the oriented stopping time to reach a value closer
    to 0. If true, the lambda will simply check equality to 1.
"""
function __stopping_time_terminus(n::Integer, total_stop::Bool)
    # https://docs.julialang.org/en/v1/manual/functions/#man-anonymous-functions
    if total_stop
        return (function (x) x == 1 end)
    elseif n >= 0
        return (function (x) x < n && x > 0 end)
    else
        return (function (x) x > n && x < 0 end)
    end
end

#TODO: Make the hailstone calc arbitrary integer safe!
"""
    hailstone_sequence(initial_value; P=2, a=3, b=1, max_total_stopping_time=1000, total_stopping_time=true, verbose=true)

Returns a list of successive values obtained by iterating a Collatz-esque function.

Until either 1 is reached, or the total amount of iterations
exceeds max_total_stopping_time. Unless total_stopping_time is False,
which will terminate the hailstone at the "stopping time" value, i.e. the
first value less than the initial value. While the sequence has the
capability to determine that it has encountered a cycle, the cycle from "1"
wont be attempted or reported as part of a cycle, regardless of default or
custom parameterisation, as "1" is considered a "total stop".

# Args
- `initial_value::Integer`: The value to begin the hailstone sequence from.

# Kwargs
- `P::Integer=2`: Modulus used to devide n, iff n is equivalent to (0 mod P).
- `a::Integer=3`: Factor by which to multiply n.
- `b::Integer=1`: Value to add to the scaled value of n.
- `max_total_stopping_time::Integer=1000`: Maximum amount of times to iterate the
    function, if 1 is not reached.
- `total_stopping_time::Bool=true`: Whether or not to execute until the "total"
    stopping time (number of iterations to obtain 1) rather than the regular stopping
    time (number of iterations to reach a value less than the initial value).
- `verbose::Bool=true`: If set to verbose, the hailstone sequence will include
    control string sequences to provide information about how the sequence
    terminated, whether by reaching a stopping time or entering a cycle.

# Examples
```jldoctest
julia> print(hailstone_sequence(16, verbose=false))
[16, 8, 4, 2, 1]
```
```jldoctest
julia> print(hailstone_sequence(16))
Any[16, 8, 4, 2, 1, Any[Collatz._CC.TOTAL_STOPPING_TIME, 4]]
```
```jldoctest
julia> print(hailstone_sequence(761, P=5, a=2, b=3, verbose=false))
[761, 1525, 305, 61, 125, 25, 5, 1]
```
# Example cycle!
```jldoctest
julia> print(hailstone_sequence(-56))
Any[-56, -28, Collatz._CC.CYCLE_INIT, Any[-14, -7, -20, -10, -5], Any[Collatz._CC.CYCLE_LENGTH, 5]]
```
```jldoctest
julia> print(hailstone_sequence(-56, verbose=false))
[-56, -28, -14, -7, -20, -10, -5, -14]
```
"""
function hailstone_sequence(initial_value::Integer; P::Integer=2, a::Integer=3, b::Integer=1, max_total_stopping_time::Integer=1000, total_stopping_time::Bool=true, verbose::Bool=true)
    # Call out the collatz_function before any magic returns to trap bad values.
    _ = collatz_function(initial_value,P=P,a=a,b=b)
    # 0 is always an immediate stop.
    if initial_value == 0; if verbose; return [[_CC.ZERO_STOP, 0]]; else; return [0]; end; end
    # 1 is always an immediate stop, with 0 stopping time.
    if initial_value == 1; if verbose; return [[_CC.TOTAL_STOPPING_TIME, 0]]; else; return [1]; end; end
    terminate = __stopping_time_terminus(initial_value, total_stopping_time)
    # Start the hailstone sequence.
    _max_total_stopping_time = max(max_total_stopping_time, 1)
    hailstone = [initial_value] # Initialises as Vector{Integer}
    if verbose; hailstone = Vector{Any}(hailstone); end
    cyclic = (function (x) x in hailstone end)
    for k in 1:_max_total_stopping_time
        _next = collatz_function(last(hailstone),P=P,a=a,b=b)
        # Check if the next hailstone is either the stopping time, total
        # stopping time, the same as the initial value, or stuck at zero.
        if terminate(_next)
            push!(hailstone, _next)
            if verbose
                if _next == 1
                    m = _CC.TOTAL_STOPPING_TIME
                else
                    m = _CC.STOPPING_TIME
                end
                push!(hailstone, [m, length(hailstone)-1])
            end
            break
        end
        if cyclic(_next)
            cycle_init = 0
            for j in 0:(length(hailstone)-1)
                if hailstone[length(hailstone)-j] == _next
                    # Because Julia is 1-based indexed, add 1 to the J
                    cycle_init = j+1
                    break
                end
            end
            if verbose
                hailstone = append!(hailstone[1:(length(hailstone)-cycle_init)], [_CC.CYCLE_INIT, hailstone[(length(hailstone)-cycle_init+1):length(hailstone)], [_CC.CYCLE_LENGTH, cycle_init]])
            else
                push!(hailstone, _next)
            end
            break
        end
        if _next == 0
            push!(hailstone, 0)
            if verbose
                push!(hailstone, [_CC.ZERO_STOP, -(length(hailstone)-1)])
            end
            break
        end
        push!(hailstone, _next)
    end
    # Julia has no for ~ else yet https://github.com/JuliaLang/julia/issues/1289
    if length(hailstone) >= max_total_stopping_time
        if verbose
            push!(hailstone, [_CC.MAX_STOP_OOB, _max_total_stopping_time])
        end
    end
    return hailstone
end

"""
    stopping_time(initial_value; P=2, a=3, b=1, max_stopping_time=1000, total_stopping_time=false)

Returns the stopping time, the amount of iterations required to reach a
value less than the initial value, or nothing if max_stopping_time is exceeded.

Alternatively, if total_stopping_time is True, then it will instead count
the amount of iterations to reach 1. If the sequence does not stop, but
instead ends in a cycle, the result will be infinity. If (P,a,b) are such
that it is possible to get stuck on zero, the result will be the negative of
what would otherwise be the "total stopping time" to reach 1, where 0 is
considered a "total stop" that should not occur as it does form a cycle of
length 1.

# Args
- `initial_value::Integer`: The value for which to find the stopping time.

# Kwargs
- `P::Integer=2`: Modulus used to devide n, iff n is equivalent to (0 mod P).
- `a::Integer=3`: Factor by which to multiply n.
- `b::Integer=1`: Value to add to the scaled value of n.
- `max_stopping_time::Integer=1000`: Maximum amount of times to iterate the function, if the stopping
    time is not reached. IF the max_stopping_time is reached, the function will return nothing.
- `total_stopping_time::Bool=false`: Whether or not to execute until the "total"
    stopping time (number of iterations to obtain 1) rather than the regular
    stopping time (number of iterations to reach a value less than the initial value).

# Examples
```jldoctest
julia> stopping_time(5)
3
julia> stopping_time(27)
96
julia> stopping_time(21, P=5, a=2, b=3, total_stopping_time=true)
Inf
```
"""
function stopping_time(initial_value::Integer; P::Integer=2, a::Integer=3, b::Integer=1, max_stopping_time::Integer=1000, total_stopping_time::Bool=false)
    # The information is contained in the verbose form of a hailstone sequence.
    # Although the "max_~_time" for hailstones is name for "total stopping" time
    # and the "max_~_time" for this "stopping time" function is _not_ "total",
    # they are handled the same way, as the default for "total_stopping_time"
    # for hailstones is true, but for this, is false. Thus the naming difference
    end_msg = last(hailstone_sequence(initial_value, P=P, a=a, b=b, verbose=true, max_total_stopping_time=max_stopping_time, total_stopping_time=total_stopping_time))
    # For total/regular/zero stopping time, the value is already the same as
    # that present, for cycles we report infinity instead of the cycle length,
    # and for max stop out of bounds, we report nothing instead of the max stop cap
    return get(Dict(_CC.TOTAL_STOPPING_TIME => end_msg[2],
                    _CC.STOPPING_TIME => end_msg[2],
                    _CC.CYCLE_LENGTH => Inf, # infinity
                    _CC.ZERO_STOP => end_msg[2],
                    _CC.MAX_STOP_OOB => nothing,
                    ), end_msg[1], nothing)
end

# TODO: Determine why the output of print(tree_graph(1, 3)) in the example in the below
# yields a "Dict{Any, Any}(Collatz._CC.CYCLE_INIT => 1, 8 => Dict{Any, Any}())" in jldoctest
# but a "Dict{Any, Any}(8 => Dict{Any, Any}(), Collatz._CC.CYCLE_INIT => 1)" from REPL?
# Both examples are covered in tests, which doesn't appear to care about order, but jldoctest
# is just a stringy parser?? so can't jldoctest these two.
"""
    tree_graph(initial_value, max_orbit_distance; P=2, a=3, b=1, __cycle_prevention=nothing)

Returns nested dictionaries that model the directed tree graph up to a
maximum nesting of max_orbit_distance, with the initial_value as the root.

# Args
- `initial_value::Integer`: The root value of the directed tree graph.
- `max_orbit_distance::Integer`: Maximum amount of times to iterate the reverse
    function. There is no natural termination to populating the tree graph, equivalent
    to the termination of hailstone sequences or stopping time attempts, so this is not
    an optional argument like max_stopping_time / max_total_stopping_time, as it is the
    intended target of orbits to obtain, rather than a limit to avoid uncapped computation.

# Kwargs
- `P::Integer=2`: Modulus used to devide n, iff n is equivalent to (0 mod P).
- `a::Integer=3`: Factor by which to multiply n.
- `b::Integer=1`: Value to add to the scaled value of n.

# Internal Kwargs
- `__cycle_prevention::Union{Set{Integer},Nothing}=nothing`: Used to prevent cycles
    from precipitatingby keeping track of all values added across previous nest depths.

# Examples
```
julia> print(tree_graph(1, 3))
Dict{Int64, Dict{Any, Any}}(1 => Dict(2 => Dict{Any, Any}(4 => Dict{Any, Any}(Collatz._CC.CYCLE_INIT => 1, 8 => Dict{Any, Any}()))))
```
```
julia> print(tree_graph(4, 3))
Dict{Int64, Dict{Any, Any}}(4 => Dict(8 => Dict{Any, Any}(16 => Dict{Any, Any}(5 => Dict{Any, Any}(), 32 => Dict{Any, Any}())), 1 => Dict{Any, Any}(2 => Dict{Any, Any}(Collatz._CC.CYCLE_INIT => 4))))
```
"""
function tree_graph(initial_value::Integer, max_orbit_distance::Integer; P::Integer=2, a::Integer=3, b::Integer=1, __cycle_prevention::Union{Set{Integer},Nothing}=nothing)
    # Call out the reverse_collatz_function before any magic returns to trap bad values.
    _ = reverse_collatz_function(initial_value,P=P,a=a,b=b)
    tgraph = Dict(initial_value=>Dict())
    if max(0, max_orbit_distance) == 0; return tgraph; end
    # Handle cycle prevention for recursive calls ~
    # Shouldn't use a mutable object initialiser for a default.
    if __cycle_prevention == nothing; __cycle_prevention = Set{Integer}(); end
    push!(__cycle_prevention, initial_value)
    for branch_value in reverse_collatz_function(initial_value, P=P, a=a, b=b)
        if branch_value in __cycle_prevention
            tgraph[initial_value][_CC.CYCLE_INIT] = branch_value
        else
            tgraph[initial_value][branch_value] = tree_graph(branch_value, max_orbit_distance-1, P=P, a=a, b=b, __cycle_prevention=__cycle_prevention)[branch_value]
        end
    end
    return tgraph
end

end # module

max_stars_repo_name: Skenvy/Collatz

max_issues_repo_name: Skenvy/Collatz

max_forks_repo_name: Skenvy/Collatz


text:
using DataConvenience
using Test

include("canonicalcor.jl")
include("janitor.jl")
include("read-csv-in-chunks.jl")
include("test-fsort-dataframes.jl")
include("missing.jl")
include("sample.jl")

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

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
using SeisProcessing
using Test

nt=250;
dt=0.004;
nx=30;
dx=5;

param = Dict(:nt=>nt, :dt=>dt,:nx1=>nx,:dx1=>dx,
              :p1=>[0.0003, -0.0005, 0.000], :p2=>[0,0,0],:p3=>[0,0,0],:p4=>[0,0,0],
              :tau=>[0.3, 0.4,0.5], :amp=>[ -1.0, 0.7,1],:f0=>20);

d=SeisLinearEvents(;param...);

h = collect(0:nx-1)*dx;

# 2- Define linear moveouts p
pmin = -2*dt/dx;
pmax = 2*dt/dx;
np = 60
dp = (pmax-pmin)/(np-1)
p = collect(pmin:dp:pmax);

# 3- Transform to Radon domain

m = SeisRadonFreqInv(d, order="linear", dt=dt, h=h, p=p, flow=2, fhigh=80,
                     mu=0.00001 )

# 4- Transform back to data domain
d2 = SeisRadonFreqFor(m, nt, order="linear", dt=dt, h=h, p=p, flow=2, fhigh=80)

# 5- Calculate and print modeling error (%)
e2 = 100*(sum((d-d2).^2))/sum(d.^2);
println("Relative error = ", e2, " %")
@test e2 < 0.1


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using DifferentialEquations, DiffEqCallbacks

############ Model-Exchange ############

# Read next time event from fmu and provide it to the integrator
function time_choice(c::fmi2Component, integrator)
    eventInfo = fmi2NewDiscreteStates(c)
    fmi2EnterContinuousTimeMode(c)
    if Bool(eventInfo.nextEventTimeDefined)
        eventInfo.nextEventTime
    else
        Inf
    end
end

# Handles events and returns the values and nominals of the changed continuous states.
function handleEvents(c::fmi2Component, enterEventMode::Bool, exitInContinuousMode::Bool)

    if enterEventMode
        fmi2EnterEventMode(c)
    end

    eventInfo = fmi2NewDiscreteStates(c)

    valuesOfContinuousStatesChanged = eventInfo.valuesOfContinuousStatesChanged
    nominalsOfContinuousStatesChanged = eventInfo.nominalsOfContinuousStatesChanged

    #set inputs here
    #fmiSetReal(myFMU, InputRef, Value)

    while eventInfo.newDiscreteStatesNeeded == fmi2True
        # update discrete states
        eventInfo = fmi2NewDiscreteStates(c)
        valuesOfContinuousStatesChanged = eventInfo.valuesOfContinuousStatesChanged
        nominalsOfContinuousStatesChanged = eventInfo.nominalsOfContinuousStatesChanged

        if eventInfo.terminateSimulation == fmi2True
            @error "Event info returned error!"
        end
    end

    if exitInContinuousMode
        fmi2EnterContinuousTimeMode(c)
    end

    return valuesOfContinuousStatesChanged, nominalsOfContinuousStatesChanged

end

# Returns the event indicators for an FMU.
function condition(c::fmi2Component, out, x, t, integrator) # Event when event_f(u,t) == 0

    fmi2SetTime(c, t)
    fmi2SetContinuousStates(c, x)
    indicators = fmi2GetEventIndicators(c)

    copy!(out, indicators)

end

# Handles the upcoming events.
function affectFMU!(c::fmi2Component, integrator, idx)
    # Event found - handle it
    continuousStatesChanged, nominalsChanged = handleEvents(c, true, Bool(sign(idx)))

    if continuousStatesChanged == fmi2True
        integrator.u = fmi2GetContinuousStates(c)
    end

    if nominalsChanged == fmi2True
        x_nom = fmi2GetNominalsOfContinuousStates(c)
    end
    timeEventCb = PresetTimeCallback(2.0, (integrator) -> affect!(c, integrator, 0))
end

# Does one step in the simulation.
function stepCompleted(c::fmi2Component, x, t, integrator)

     (status, enterEventMode, terminateSimulation) = fmi2CompletedIntegratorStep(c, fmi2True)
     if enterEventMode == fmi2True
        affect!(c, integrator, 0)
     end

end

# Returns the state derivatives of the FMU.
function fx(c::fmi2Component, x, p, t)
    fmi2SetTime(c, t)
    fmi2SetContinuousStates(c, x)
    dx = fmi2GetDerivatives(c)
end

# save FMU values
function saveValues(c::fmi2Component, recordValues, u, t, integrator)
    Tuple(fmiGetReal(c, recordValues)...,)
end

"""
Source: FMISpec2.0.2[p.90 ff]: 3.2.4 Pseudocode Example

Simulates a FMU instance for the given simulation time interval.
State- and Time-Events are handled correctly.

Returns a tuple of type (ODESolution, DiffEqCallbacks.SavedValues).
If keyword `recordValues` is not set, a tuple of type (ODESolution, nothing) is returned for consitency.
"""
function fmi2SimulateME(c::fmi2Component, t_start::Real = 0.0, t_stop::Real = 1.0;
    solver = nothing,
    customFx = nothing,
    recordValues::fmi2ValueReferenceFormat = nothing,
    saveat = [],
    setup = true)

    recordValues = prepareValueReference(c, recordValues)
    solution = nothing
    callbacks = []
    savedValues = nothing

    if length(recordValues) > 0
        savedValues = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(recordValues))...})

        savingCB = SavingCallback((u,t,integrator) -> saveValues(c, recordValues, u, t, integrator),
                                  savedValues,
                                  saveat=saveat)
        push!(callbacks, savingCB)
    end

    if solver == nothing
        solver = Tsit5()
    end

    if customFx == nothing
        customFx = (x, p, t) -> fx(c, x, p, t)
    end

    if setup
        fmi2Reset(c)
        fmi2SetupExperiment(c, t_start, t_stop)
        fmi2EnterInitializationMode(c)
        fmi2ExitInitializationMode(c)
    end

    eventHandling = c.fmu.modelDescription.numberOfEventIndicators > 0

    # First evaluation of the FMU
    x0 = fmi2GetContinuousStates(c)
    x0_nom = fmi2GetNominalsOfContinuousStates(c)

    fmi2SetContinuousStates(c, x0)

    handleEvents(c, false, false)

    # Get states of handling initial Events
    x0 = fmi2GetContinuousStates(c)
    x0_nom = fmi2GetNominalsOfContinuousStates(c)

    p = []
    problem = ODEProblem(customFx, x0, (t_start, t_stop), p,)

    if eventHandling

        eventInfo = fmi2NewDiscreteStates(c)
        fmi2EnterContinuousTimeMode(c)

        timeEvents = (eventInfo.nextEventTimeDefined == fmi2True)

        eventCb = VectorContinuousCallback((out, x, t, integrator) -> condition(c, out, x, t, integrator),
                                           (integrator, idx) -> affectFMU!(c, integrator, idx),
                                           Int64(c.fmu.modelDescription.numberOfEventIndicators);
                                           rootfind = DiffEqBase.RightRootFind)
        push!(callbacks, eventCb)

        stepCb = FunctionCallingCallback((x, t, integrator) -> stepCompleted(c, x, t, integrator);
                                         func_everystep = true,
                                         func_start = true)
        push!(callbacks, stepCb)

        if timeEvents
            timeEventCb = IterativeCallback((integrator) -> time_choice(c, integrator),
                                            (integrator) -> affectFMU!(c, integrator, 0), Float64; initial_affect = true)

            push!(callbacks, timeEventCb)
        end
    end

    if length(callbacks) > 0
        solution = solve(problem, solver, callback = CallbackSet(callbacks...), saveat = saveat)
    else
        solution = solve(problem, solver, saveat = saveat)
    end

    return solution, savedValues
end

############ Co-Simulation ############

"""
Starts a simulation of the Co-Simulation FMU instance.

Returns a tuple of (success::Bool, DiffEqCallbacks.SavedValues) with success = `true` or `false`.
If keyword `recordValues` is not set, a tuple of type (success::Bool, nothing) is returned for consitency.

ToDo: Improove Documentation.
"""
function fmi2SimulateCS(c::fmi2Component, t_start::Real, t_stop::Real;
                        recordValues::fmi2ValueReferenceFormat = nothing,
                        saveat = [],
                        setup = true)

    recordValues = prepareValueReference(c, recordValues)
    variableSteps = c.fmu.modelDescription.CScanHandleVariableCommunicationStepSize

    success = false
    savedValues = nothing

    # default setup
    if length(saveat) == 0
        saveat = LinRange(t_start, t_stop, 100)
    end
    dt = (t_stop - t_start) / length(saveat)

    # setup if no variable steps
    if variableSteps == false
        if length(saveat) >= 2
            dt = saveat[2] - saveat[1]
        end
    end

    if setup
        fmi2Reset(c)
        fmi2SetupExperiment(c, t_start, t_stop)
        fmi2EnterInitializationMode(c)
        fmi2ExitInitializationMode(c)
    end

    t = t_start

    record = length(recordValues) > 0

    #numDigits = length(string(round(Integer, 1/dt)))

    if record
        savedValues = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(recordValues))...})

        i = 1

        values = (fmi2GetReal(c, recordValues)...,)
        DiffEqCallbacks.copyat_or_push!(savedValues.t, i, t)
        DiffEqCallbacks.copyat_or_push!(savedValues.saveval, i, values, Val{false})

        while t < t_stop
            if variableSteps
                if length(saveat) > i
                    dt = saveat[i+1] - saveat[i]
                else
                    dt = t_stop - saveat[i]
                end
            end

            fmi2DoStep(c, t, dt)
            t = t + dt #round(t + dt, digits=numDigits)
            i += 1

            values = (fmi2GetReal(c, recordValues)...,)
            DiffEqCallbacks.copyat_or_push!(savedValues.t, i, t)
            DiffEqCallbacks.copyat_or_push!(savedValues.saveval, i, values, Val{false})
        end

        success = true
    else
        i = 1
        while t < t_stop
            if variableSteps
                if length(saveat) > i
                    dt = saveat[i+1] - saveat[i]
                else
                    dt = t_stop - saveat[i]
                end
            end

            fmi2DoStep(c, t, dt)
            t = t + dt #round(t + dt, digits=numDigits)
            i += 1
        end

        success = true
    end

    success, savedValues
end

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

"""
Starts a simulation of the fmu instance for the matching fmu type, if both types are available, CS is preferred.

Returns:
    - a tuple of (success::Bool, DiffEqCallbacks.SavedValues) with success = `true` or `false` for CS-FMUs
    - a tuple of (ODESolution, DiffEqCallbacks.SavedValues) for ME-FMUs
    - if keyword `recordValues` is not set, a tuple of type (..., nothing)

ToDo: Improove Documentation.
"""
function fmi2Simulate(c::fmi2Component, t_start::Real = 0.0, t_stop::Real = 1.0;
                      recordValues::fmi2ValueReferenceFormat = nothing,
                      saveat = [],
                      setup = true)

    if fmi2IsCoSimulation(c.fmu)
        return fmi2SimulateCS(c, t_start, t_stop; recordValues=recordValues, saveat=saveat, setup=setup)
    elseif fmi2IsModelExchange(c.fmu)
        return fmi2SimulateME(c, t_start, t_stop; recordValues=recordValues, saveat=saveat, setup=setup)
    else
        error(unknownFMUType)
    end
end


text:
using Base.Test

using ErlPort.ErlTerms


text:
using SpecialFunctions, Roots
using Parameters, Distributions
using Plots, Random
const sf = SpecialFunctions

export DirichletPDFParams, ValencePDFParams
export DIRICHLET_TYPE, VALENCE_TYPE
export get_scaled_θ
export plot_input_pdfs, int_xtotx, xtotx
export get_input_pdf_func
export input_pdf_map

VALENCE_TYPE = 1
DIRICHLET_TYPE = 2

"""
    Abstract type for any PDF parametrisation.
"""
abstract type AbstractPDFParams end

"""
    struct DirichletPDFParams <: AbstractPDFParams

Full Dirichlet specification of input PDF parameters.
"""
@with_kw struct DirichletPDFParams <: AbstractPDFParams
    param_type::Integer = DIRICHLET_TYPE
    seed::Integer = 0
    weights::Vector{Float64} = [1, 1, 1, 1, 1, 1, 1, 1, 1]
    θ::Vector{Float64} = rand(MersenneTwister(seed), Dirichlet(weights))
    K_u::Float64
    λ_u::Float64 = (θ[1] * (K_u + 1)) / (2 - θ[1])
    K_d::Float64
    λ_d::Float64 = (θ[2] * (K_d + 1)) / (1 - θ[2])
    λ_g1::Float64
    λ_g2::Float64
    K_g::Float64
    λ_q::Float64
end

"""
    struct ValencePDFParams <: AbstractPDFParams

Valence specification of input PDF parameters.
"""
@with_kw struct ValencePDFParams <: AbstractPDFParams
    param_type::Integer = VALENCE_TYPE
    λ_u::Float64
    K_u::Float64
    λ_d::Float64
    K_d::Float64
    λ_g1::Float64
    λ_g2::Float64
    K_g::Float64
    λ_q::Float64
    seed::Integer = 0
    weights::Vector{Float64} = [1, 1, 1, 1, 1, 1, 1]
    θ::Vector{Float64} = get_dirichlet_samples(λ_u, K_u, λ_d, K_d, seed, weights)
end

"""
    get_scaled_θ(λ_u, K_u, λ_d, K_d, θ_tmp)

Given a set of Dirichlet samples, θ_tmp, scale
according to the valence params. Relevant for
`ValencePDFParams`.
"""
function get_scaled_θ(λ_u::Float64, K_u::Float64, λ_d::Float64, K_d::Float64,
                      θ_tmp::Vector{Float64})

    A_u = 2 / sf.beta(λ_u, K_u + 1)
    A_d = 1 / sf.beta(λ_d, K_d + 1)

    I_u = A_u * sf.beta(λ_u + 1, K_u + 1)
    I_d = A_d * sf.beta(λ_d + 1, K_d + 1)

    remaining = 1 - (I_u + I_d)

    return θ_tmp * remaining

end

"""
    get_dirichlet_samples(λ_u, K_u, λ_d, K_d, seed, weights)

Given valance shape parameters and weights, get samples θ
that are correctly scaled. Relevant for `ValencePDFParams`
"""
function get_dirichlet_samples(λ_u::Float64, K_u::Float64, λ_d::Float64,
                               K_d::Float64, seed::Integer, weights::Vector{Float64})

    Random.seed!(seed);
    θ_tmp = rand(Dirichlet(weights))

    θ = get_scaled_θ(λ_u, K_u, λ_d, K_d, θ_tmp)

    return θ

end

"""
    x_uv_x(x, λ_u, K_u)

Momentum density of u valence component.
Beta function
    A_u x^λ_u (1-x)^K_u
A_u is set by λ_u and K_u.
"""
function x_uv_x(x::Float64, λ_u::Float64, K_u::Float64)

    A_u = 2 / sf.beta(λ_u, K_u + 1)

    return A_u * x^λ_u * (1 - x)^K_u

end

"""
    x_dv_x(x, λ_d, K_d)

Momentum density of d valence component.
Beta function
    A_d x^λ_d (1 - x)^K_d
A_d is set by λ_d and K_d.
"""
function x_dv_x(x::Float64, λ_d::Float64, K_d::Float64)

    A_d = 1 / sf.beta(λ_d, K_d + 1)

    return A_d * x^λ_d * (1 - x)^K_d

end

"""
        x_g_x(x, λ_g1, λ_g2, K_g, w1, w2)

Momentum density of gluon component.
    A_g1 x^λ_g1 (1 - x)^K_g + A_g2 x^λ_g2
Amplitudes are set by weights `w1` and `w2`.
"""
function x_g_x(x::Float64, λ_g1::Float64, λ_g2::Float64, K_g::Float64,
               w1::Float64, w2::Float64)

    A_g1 = w1 / sf.beta(λ_g1 + 1, K_g + 1)

    x_g1_x = A_g1 * x^λ_g1 * (1 - x)^K_g

    A_g2 = w2 / sf.beta(λ_g2 + 1, 5 + 1)

    x_g2_x = A_g2 * x^λ_g2 * (1 - x)^5

    return x_g1_x + x_g2_x

end

"""
    x_q_x(x, λ_q, w)

Momentum density of non-valence
quark component.
    A_q x^λ_q
Amplitude is set by weight `w`
"""
function x_q_x(x::Float64, λ_q::Float64, w::Float64)

    A_q = (w / 2) / sf.beta(λ_q + 1, 5 + 1)

    return A_q * x^λ_q * (1 - x)^5

end

"""
    x_total_x(x, λ_u, K_u, λ_d, K_d, λ_g1, λ_g2, K_g, λ_q, θ)

Total momentum density.
"""
function xtotx(x::Float64, λ_u::Float64, K_u::Float64, λ_d::Float64, K_d::Float64,
               λ_g1::Float64, λ_g2::Float64, K_g::Float64, λ_q::Float64, θ::Vector{Float64})

    xuvx = x_uv_x(x, λ_u, K_u)

    xdvx = x_dv_x(x, λ_d, K_d)

    # For ValencePDFParams
    if length(θ) == 7

        xgx = x_g_x(x, λ_g1, λ_g2, K_g, θ[1], θ[2])

        xubarx = x_q_x(x, λ_q, θ[3])

        xdbarx = x_q_x(x, λ_q, θ[4])

        xsx = x_q_x(x, λ_q, θ[5])

        xcx = x_q_x(x, λ_q, θ[6])

        xbx = x_q_x(x, λ_q, θ[7])

    # For DirichletPDFParams
    elseif length(θ) == 9

        xgx = x_g_x(x, λ_g1, λ_g2, K_g, θ[3], θ[4])

        xubarx = x_q_x(x, λ_q, θ[5])

        xdbarx = x_q_x(x, λ_q, θ[6])

        xsx = x_q_x(x, λ_q, θ[7])

        xcx = x_q_x(x, λ_q, θ[8])

        xbx = x_q_x(x, λ_q, θ[9])

    else

        @error "length(θ) must be 7 or 9 for a valid PDF parametrisation"

    end

    return xuvx + xdvx + xgx + xubarx + xdbarx + xsx + xcx + xbx

end

"""
    x_total_x(x, pdf_params)

Total momentum density.
"""
function xtotx(x::Float64, pdf_params::AbstractPDFParams)

    pdf = pdf_params

    return xtotx(x, pdf.λ_u, pdf.K_u, pdf.λ_d, pdf.K_d, pdf.λ_g1,
                 pdf.λ_g2, pdf.K_g, pdf.λ_q, pdf.θ)
end

"""
    int_xtotx(λ_u, K_u, λ_d, K_d, λ_g1, λ_g2, K_g, λ_q, θ)

Total integrated momentum density. Should equal 1.
"""
function int_xtotx(λ_u::Float64, K_u::Float64, λ_d::Float64, K_d::Float64,
                   λ_g1::Float64, λ_g2::Float64, K_g::Float64, λ_q::Float64, θ::Array{Float64})

    A_u = 2 / sf.beta(λ_u, K_u + 1)
    A_d = 1 / sf.beta(λ_d, K_d + 1)

    # ValencePDFParams
    if length(θ) == 7

        A_g1 = θ[1] / sf.beta(λ_g1 + 1, K_g + 1)
        A_g2 = θ[2] / sf.beta(λ_g2 + 1, 5 + 1)

        A_ubar = (θ[3] / 2) / sf.beta(λ_q + 1, 5 + 1)
        A_dbar = (θ[4] / 2) / sf.beta(λ_q + 1, 5 + 1)

        A_s = (θ[5] / 2) / sf.beta(λ_q + 1, 5 + 1)
        A_c = (θ[6] / 2) / sf.beta(λ_q + 1, 5 + 1)
        A_b = (θ[7] / 2) / sf.beta(λ_q + 1, 5 + 1)

    # DirichletPDFParams
    elseif length(θ) == 9

        A_g1 = θ[3] / sf.beta(λ_g1 + 1, K_g + 1)
        A_g2 = θ[4] / sf.beta(λ_g2 + 1, 5 + 1)

        A_ubar = (θ[5] / 2) / sf.beta(λ_q + 1, 5 + 1)
        A_dbar = (θ[6] / 2) / sf.beta(λ_q + 1, 5 + 1)

        A_s = (θ[7] / 2) / sf.beta(λ_q + 1, 5 + 1)
        A_c = (θ[8] / 2) / sf.beta(λ_q + 1, 5 + 1)
        A_b = (θ[9] / 2) / sf.beta(λ_q + 1, 5 + 1)

    else

         @error "length(θ) must be 7 or 9 for a valid PDF parametrisation"

    end

    I_u = A_u * sf.beta(λ_u + 1, K_u + 1)
    I_d = A_d * sf.beta(λ_d + 1, K_d + 1)
    I_g = A_g1 * sf.beta(λ_g1 + 1, K_g + 1) + A_g2 * sf.beta(λ_g2 + 1, 5 + 1)
    I_q = 2 * (A_ubar + A_dbar + A_s + A_c + A_b) * sf.beta(λ_q + 1, 5 + 1)

    I_tot = I_u + I_d + I_g + I_q

    I_tot
end

"""
    int_xtotx(pdf_params)

Total integrated momentum density. Should equal 1.
"""
function int_xtotx(pdf_params::AbstractPDFParams)

    pdf = pdf_params

    result = int_xtotx(pdf.λ_u, pdf.K_u, pdf.λ_d, pdf.K_d,
                       pdf.λ_g1, pdf.λ_g2, pdf.K_g, pdf.λ_q, pdf.θ)

    return result
end

"""
    plot_input_pdfs(pdf_params, xmin, xmax, nx)

Plot the input PDFs defined by hyper_params over
the given x range.
"""
function plot_input_pdfs end

function plot_input_pdfs(pdf_params::ValencePDFParams, xmin::Float64=1.0e-2,
                         xmax::Float64=1.0, nx::Integer=1000)

    x_grid = range(xmin, stop=xmax, length=nx)
    pdf = pdf_params

    p = plot(x_grid, [x_uv_x(x, pdf.λ_u, pdf.K_u) for x in x_grid], label="x uv(x)", lw=3)
    p = plot!(x_grid, [x_dv_x(x, pdf.λ_d, pdf.K_d) for x in x_grid], label="x dv(x)", lw=3)

    p = plot!(x_grid, [x_g_x(x, pdf.λ_g1, pdf.λ_g2, pdf.K_g, pdf.θ[1], pdf.θ[2])
                       for x in x_grid], label="x g(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[3]) for x in x_grid], label="x ubar(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[4]) for x in x_grid], label="x dbar(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[5]) for x in x_grid], label="x s(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[6]) for x in x_grid], label="x c(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[7]) for x in x_grid], label="x b(x)", lw=3)

    p = plot!(xlabel="x")
    p = plot!(xaxis=:log, yaxis=:log, legend=:outertopright)
    p = ylims!(1e-8, 30)

    return p
end

function plot_input_pdfs(pdf_params::DirichletPDFParams, xmin::Float64=1.0e-2,
                         xmax::Float64=1.0, nx::Integer=1000)

    x_grid = range(xmin, stop=xmax, length=nx)
    pdf = pdf_params

    p = plot(x_grid, [x_uv_x(x, pdf.λ_u, pdf.K_u) for x in x_grid], label="x uv(x)", lw=3)
    p = plot!(x_grid, [x_dv_x(x, pdf.λ_d, pdf.K_d) for x in x_grid], label="x dv(x)", lw=3)

    p = plot!(x_grid, [x_g_x(x, pdf.λ_g1, pdf.λ_g2, pdf.K_g, pdf.θ[3], pdf.θ[4])
                       for x in x_grid], label="x g(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[5]) for x in x_grid], label="x ubar(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[6]) for x in x_grid], label="x dbar(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[7]) for x in x_grid], label="x s(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[8]) for x in x_grid], label="x c(x)", lw=3)
    p = plot!(x_grid, [x_q_x(x, pdf.λ_q, pdf.θ[9]) for x in x_grid], label="x b(x)", lw=3)

    p = plot!(xlabel="x")
    p = plot!(xaxis=:log, yaxis=:log, legend=:outertopright)
    p = ylims!(1e-8, 30)

    return p
end

"""
    get_input_pdf_func(pdf_params)

Get the function to input the PDFs into QCDNUM.
"""
function get_input_pdf_func end

function get_input_pdf_func(pdf_params::ValencePDFParams)::Function

    pdf = pdf_params

    func = function _input_pdfs(i, x)::Float64
        i = i[]
        x = x[]

        f = 0.0

        # gluon
        if (i == 0)
            f = x_g_x(x, pdf.λ_g1, pdf.λ_g2, pdf.K_g, pdf.θ[1], pdf.θ[2])
        end

        # u valence
        if (i == 1)
            f = x_uv_x(x, pdf.λ_u, pdf.K_u)
        end

        # d valence
        if (i == 2)
            f = x_dv_x(x, pdf.λ_d, pdf.K_d)
        end

        # ubar
        if (i == 3)
            f = x_q_x(x, pdf.λ_q, pdf.θ[3])
        end

        # dbar
        if (i == 4)
            f = x_q_x(x, pdf.λ_q, pdf.θ[4])
        end

        # s and sbar
        if (i == 5) || (i == 6)
            f = x_q_x(x, pdf.λ_q, pdf.θ[5])
        end

        # c and cbar
        if (i == 7) || (i == 8)
            f = x_q_x(x, pdf.λ_q, pdf.θ[6])
        end

        # d and dbar
        if (i == 9) || (i == 10)
            f = x_q_x(x, pdf.λ_q, pdf.θ[7])
        end

        return f
    end

    return func
end

function get_input_pdf_func(pdf_params::DirichletPDFParams)::Function

    pdf = pdf_params

    func = function _input_pdfs(i, x)::Float64
        i = i[]
        x = x[]

        f = 0.0

        # gluon
        if (i == 0)
            f = x_g_x(x, pdf.λ_g1, pdf.λ_g2, pdf.K_g, pdf.θ[3], pdf.θ[4])
        end

        # u valence
        if (i == 1)
            f = x_uv_x(x, pdf.λ_u, pdf.K_u)
        end

        # d valence
        if (i == 2)
            f = x_dv_x(x, pdf.λ_d, pdf.K_d)
        end

        # ubar
        if (i == 3)
            f = x_q_x(x, pdf.λ_q, pdf.θ[5])
        end

        # dbar
        if (i == 4)
            f = x_q_x(x, pdf.λ_q, pdf.θ[6])
        end

        # s and sbar
        if (i == 5) || (i == 6)
            f = x_q_x(x, pdf.λ_q, pdf.θ[7])
        end

        # c and cbar
        if (i == 7) || (i == 8)
            f = x_q_x(x, pdf.λ_q, pdf.θ[8])
        end

        # d and dbar
        if (i == 9) || (i == 10)
            f = x_q_x(x, pdf.λ_q, pdf.θ[9])
        end

        return f
    end

    return func
end

"""
    input_pdf_map

The relevant mapping for use with QCDNUM and `get_input_pdf_func`.
"""
#                         tb  bb  cb  sb  ub  db  g   d   u   s   c   b   t
input_pdf_map = Float64.([0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0., # 1 # U valence
                          0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0., # 2 # D valence
                          0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., # 3 # u sea
                          0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., # 4 # d sea
                          0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., # 5 # s
                          0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., # 6 # sbar
                          0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., # 7 # c
                          0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., # 8 # cbar
                          0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., # 9 # b
                          0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., # 10 # bbar
                          0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., # 11
                          0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]); # 12


text:
types = {
    Int8, Uint8, Int16, Uint16, Int32, Uint32, Int64, Uint64, Float32, Float64,
    Rational{Int8}, Rational{Uint8}, Rational{Int16}, Rational{Uint16},
    Rational{Int32}, Rational{Uint32}, Rational{Int64}, Rational{Uint64}
}
vals = [
    typemin(Int64),
    -integer(maxintfloat(Float64))+(-4:1),
    typemin(Int32),
    -integer(maxintfloat(Float32))+(-4:1),
    -2:2,
    integer(maxintfloat(Float32))+(-1:4),
    typemax(Int32),
    integer(maxintfloat(Float64))+(-1:4),
    typemax(Int64),
]

for T=types, S=types, x=vals
    a = convert(T,x)
    b = convert(S,x)
    #println("$(typeof(a)) $a")
    #println("$(typeof(b)) $b")
    @test !isequal(a,b) || hash(a)==hash(b)
    # for y=vals
    #     println("T=$T; S=$S; x=$x; y=$y")
    #     c = convert(T,x//y)
    #     d = convert(S,x//y)
    #     @test !isequal(a,b) || hash(a)==hash(b)
    # end
end

@test hash(RopeString("1","2")) == hash("12")
@test hash(SubString("--hello--",3,7)) == hash("hello")
@test hash(:(X.x)) == hash(:(X.x))
@test hash(:(X.x)) != hash(:(X.y))

@test hash([1,2]) == hash(sub([1,2,3,4],1:2))

# make sure >>> is used
@test bitmix(2, -3) != bitmix(2, -4)
@test bitmix(-3, 2) != bitmix(-4, 2)

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
using Plots
using StatsPlots
using RDatasets
using StatsBase
using DataFrames
using DataFramesMeta
cars = dataset("datasets", "mtcars")

describe(cars)

# Splitting by engine type
# AND
@linq where(:Cyl .== 4, :Gear .== 5)(cars)

@linq where(.|(:Cyl .== 4, :Cyl .== 6))(cars)
@linq where(:Model .== "Honda Civic")(cars)

@linq select(:Model, :MPG)(where(:Cyl .!= 4)(cars))

@linq orderby(:MPG)(select(:Model, :MPG)(where(:Cyl .== 4)(cars)))

@linq orderby(-:MPG)(select(:Model, :MPG)(where(:Cyl .== 4)(cars)))

@linq by(:Cyl; meanMPG=mean(:MPG))(transform(; KPL=:MPG .* 0.425144)(cars))

unstack(cars, :Model, :Cyl, :MPG)

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
using Test, Random, CUDA, Knet, AutoGrad, Statistics, Knet.Layers21

@testset "layers21/batchnorm" begin

    x, s, b, s2, b2, s3, b3 = (Param(CUDA.randn(Float64,x...)) for x in ((5,4,3,2),(1,1,3,1),(1,1,3,1),(5,4,3,1),(5,4,3,1),(5,1,1,1),(5,1,1,1)))

    # @gcheck needs Inference and Training to compute the same function
    # For this they need to use the same mean/var
    # @gcheck sums the output of non-scalar functions by default
    # This has the effect of multiplying scale with a 0 mean array
    # So we'll use getindex(y,1) instead
    function bnormtest(x,s,b; dims=(1,2,4), o...)
        m = mean(value(x); dims)
        v = var(value(x); mean=m,corrected=false,dims)
        bn = BatchNorm(; dims, mean=m, var=v, bias=b, scale=s, o...)
        getindex(bn(x), 1)
    end

    @test bnormtest(x,s,b) ≈ value(@diff(bnormtest(x,s,b)))
    @test @gcheck bnormtest(x,s,b)
    @test @gcheck bnormtest(x,s2,b2; dims=4)
    @test @gcheck bnormtest(x,s3,b3; dims=(2,3,4))
    @test @gcheck bnormtest(x,s,b; epsilon = 0)
    @test @gcheck bnormtest(x,s,b; momentum = 0.99)
end


text:
mutable struct A0 <: Model
    f0::Int
end

@testset "clean!" begin
    a = A0(0)
    @test clean!(a) == ""
    function M.clean!(a::A0)
        warn = ""
        if a.f0 < 0
            warn *= "Field a is negative, resetting to 5."
            a.f0 = 5
        end
        return warn
    end
    a = A0(-2)
    @test clean!(a) == "Field a is negative, resetting to 5."
    @test a.f0 == 5
end

# No type, no default
@mlj_model mutable struct A1
    a
end

@testset "@mlj-1" begin
    a = A1()
    @test ismissing(a.a)
    a.a = 5
    @test a.a == 5
end

# No type, with default
@mlj_model mutable struct A1b
    a = 5
end

@testset "@mlj-2" begin
    a = A1b()
    @test a.a == 5
    a.a = "hello"
    @test a.a == "hello"
end

# If a type is given but no default value is given, then the macro tries to fill
# a default value; either 0 if it's a Number type, or an empty string and otherwise fails.
@mlj_model mutable struct A1c
    a::Int
end

@testset "@mlj-3" begin
    a = A1c()
    @test a.a == 0
    a = A1c(a=7)
    @test a.a == 7
    @test_throws InexactError A1c(a=5.3)
    @test_throws MethodError A1c(a="hello")
end

# Type is given and default is given
@mlj_model mutable struct A1d
    a::Int = 5
end

@testset "@mlj-4" begin
    a = A1d()
    @test a.a == 5
    a = A1d(a=7)
    @test a.a == 7
end

# No type is given but a default and constraint
@mlj_model mutable struct A1e
    a = 5::(_ > 0)
end

@testset "@mlj-5" begin
    a = A1e()
    @test a.a == 5
    a = A1e(a=7)
    @test a.a == 7
    @test @test_logs (:warn, "Constraint `model.a > 0` failed; using default: a=5.") A1e(a=-1).a==5
    a = A1e(a=7.5)
    @test a.a == 7.5
end

# Type is given with default and constraint
@mlj_model mutable struct A1f
    a::Int = 5::(_ > 0)
end

@testset "@mlj-6" begin
    a = A1f()
    @test a.a == 5
    a = A1f(a=7)
    @test a.a == 7
    @test_throws InexactError A1f(a=7.5)
    @test @test_logs (:warn, "Constraint `model.a > 0` failed; using default: a=5.") A1f(a=-1).a==5
end

abstract type FooBar end
@mlj_model mutable struct B1a <: FooBar
    a::Symbol = :auto::(_ in (:auto, :semi))
end

@testset "@mlj-7" begin
    b = B1a()
    @test b.a == :auto
    b = B1a(a=:semi)
    @test b.a == :semi
    @test @test_logs (:warn, "Constraint `model.a in (:auto, :semi)` failed; using default: a=:auto.") B1a(a=:autos).a == :auto
    @test_throws MethodError B1a(b="blah")
end

# == dependence on other types

@mlj_model mutable struct B1b
    a::SemiMetric = Euclidean()::(_ isa Metric)
end

@mlj_model mutable struct B1c
    a::SemiMetric = Euclidean()
end

@testset "@mlj-dist" begin
    @test B1b().a isa Euclidean
    @test @test_logs (:warn, "Constraint `model.a isa Metric` failed; using default: a=Euclidean().") B1b(a=BhattacharyyaDist()).a isa Euclidean
    @test B1c().a isa Euclidean
end

# Implicit defaults
@mlj_model mutable struct Ca
    a::String
end

@mlj_model mutable struct Cb
    a::Any
end

@mlj_model mutable struct Cc
    a::Union{Nothing,Int}
end

@mlj_model mutable struct Cd
    a::Union{Missing,Int}
end

@testset "@mlj-10" begin
    @test Ca().a == ""
    @test Cb().a === missing
    @test Cc().a === nothing
    @test Cd().a === missing
end


text:
isperfect(n::Integer) = n == sum([n % i == 0 ? i : 0 for i = 1:(n - 1)])
perfects(n::Integer) = filter(isperfect, 1:n)

@show perfects(10000)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module UKrig

using SpecialFunctions: besselk, gamma
import DynamicPolynomials: @polyvar, monomials
using FixedPolynomials
using LinearAlgebra
using Statistics
using Random

export Mnu, Gnu, generate_Gnu_krigY, generate_Mnu_krigY

tν𝒦t(t,ν) = t^ν * besselk(ν, t)

function Mnu(t, ν)::Float64
	pt, pν, p0, p1 = promote(t, ν, Float64(0), Float64(1)) # note besselk always returns a Float64 apparently
	return (pt==p0) ? p1 : tν𝒦t(√(2pν)*pt,pν) * 2^(1-pν) / gamma(pν)
end

# the const on the principle irregular term
scν(ν)       = - (2ν)^ν * gamma(ν + 1//2) * gamma(1-ν) / gamma(2ν+1) / sqrt(π)
scν(ν::Int)  = - 2 * (-2ν)^ν * gamma(ν + 1//2) / gamma(ν) / gamma(2ν+3) / sqrt(π)

function Gnu(t::T, ν::Int) where T<:Real
	if t==0
		return T(0)
	end
	return scν(ν) * t^(2ν) * log(t)
end

# Gnu(t::T, ν) where T<:Real = scν(ν) * t^(2ν)
function Gnu(t::T, ν) where T<:Real
	return scν(ν) * t^(2ν)
end

"""
`generate_Gnu_krigY(;fdata, xdata, ν, σg, σε) -> (x::Array->krigY(x), x::Array->fp(x'), b, c)`
"""
function generate_Gnu_krigY(Ydata,Xdata1,Xdata2,Xdata3, ν)
	m   =  7
	n₁  = length(Xdata1)
	#kriging matrix F
	@polyvar x y z
	p=1+x+y+x^2+y^2+x*y+z
	F=[t(x=>vx,y=>vy,z=>vz) for (vx, vy,vz) in zip(Xdata1,Xdata2,Xdata3), t in p]
	#distance matrix and Gnu
 	dist_mat= sqrt.(norm.( (Xdata1.-Xdata1').^2  + (Xdata2.-Xdata2').^2) )
	Gnu_mat=Gnu.(dist_mat, ν) |>Symmetric
	#generate the nullspace of F
	Mᵀ = nullspace(F')
	M = transpose(Mᵀ)
	#pre-compute matrix
	#This one takes so much time
 	Mat₁=Symmetric(M*Gnu_mat*Mᵀ)
	Mat₂=Symmetric(M*Mᵀ)

	Ydata_M=M*Ydata

    function KrigY(σg,σε,X1,X2,X3)
		G₁₁ = (σg^2) .* Gnu_mat
		Ξ   = [
			G₁₁ .+ σε^2*I(n₁)    	 	F
			F'       zeros(m,m)
		]
		cb = Ξ \ vcat(Ydata, zeros(m))
		c  = cb[1:length(Ydata)]
		b  = cb[length(Ydata)+1:end]
		K= (σg^2) .* Gnu.( sqrt.(norm.( (X1.-Xdata1').^2  + (X2.-Xdata2').^2) )      , ν)
		FF = [t(x=>vx,y=>vy,z=>vz) for (vx,vy,vz) in zip(X1,X1,X3), t in p]
		return K*c.+FF*b
	end

	function loglike(σg,σε)

		Σ₁=σg^2*Mat₁+σε^2*Mat₂
		ch=cholesky(Σ₁)
 		return -1/2*sum( (ch.L\(Ydata_M)).^2)- sum(log.(diag(ch.L)))

	end
	return KrigY, loglike
end

end


text:
using PtFEM, Test

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 3, 2, 9, :z, Quadrilateral(4, 2), true),
  :properties => [
     100.0 0.3;
    1000.0 0.45
    ],
  :etype => [1, 2, 1, 2, 1, 2],
  :x_coords => [0.0, 4.0, 10.0, 30.0],
  :y_coords => [0.0, -4.0, -10.0],
  :support => [
    (1, [0 1]),
    (2, [0 1]),
    (3, [0 0]),
    (6, [0 0]),
    (9, [0 0]),
    (10, [0 1]),
    (11, [0 1]),
    (12, [0 0])
    ],
  :loaded_nodes => [
    (1, [0.0  -2.6667]),
    (4, [0.0 -23.3333]),
    (7, [0.0 -24.0])
    ]
)

@time m, dis_df, fm_df = p51(data)

@test m.loads[0:12] ≈ [0.0,-0.03176060466236135,-0.003231272451306888,
  0.0013949965318819561,-0.039904996293397346,0.0011648123636934574,
  -0.0024979956898815428,0.0017035565573309804,-0.006045921706318328,
  0.0013301584221295814,-0.000442142378941997,0.0025875517277645756,
  0.00030906084483465454]


text:
###############
# Branch Type #
###############

type Branch <: GitHubType
    name::Nullable{GitHubString}
    label::Nullable{GitHubString}
    ref::Nullable{GitHubString}
    sha::Nullable{GitHubString}
    commit::Nullable{Commit}
    user::Nullable{Owner}
    repo::Nullable{Repo}
    _links::Nullable{Dict}
    protection::Nullable{Dict}
end

Branch(data::Dict) = json2github(Branch, data)
Branch(name::AbstractString) = Branch(Dict("name" => name))

namefield(branch::Branch) = isnull(branch.name) ? branch.ref : branch.name

###############
# API Methods #
###############

function branches(repo; options...)
    results, page_data = gh_get_paged_json("/repos/$(name(repo))/branches"; options...)
    return map(Branch, results), page_data
end

function branch(repo, branch_obj; options...)
    result = gh_get_json("/repos/$(name(repo))/branches/$(name(branch_obj))"; options...)
    return Branch(result)
end


text:
using OctreeBH
using Test

include("test_ngbs.jl")
include("test_gravtree.jl")
include("test_balance.jl")


text:
struct BayesianEVA{T} <: fittedEVA{T}
    "Extreme value model definition"
    model::T
    "MCMC outputs"
    sim::Mamba.Chains
end

function getdistribution(fm::BayesianEVA)::Array{<:ContinuousUnivariateDistribution,2}

    v = fm.sim.value[:,:,1]

    V = Extremes.slicematrix(v, dims=2)

    D = Extremes.getdistribution.(fm.model, V)

    d = Extremes.unslicematrix(D, dims=2)

    return d

end

"""
    quantile(fm::BayesianEVA,p::Real)::Real

Compute the quantile of level `p` from the fitted Bayesian model `fm`.

If the model is stationary, then a quantile is returned for each MCMC steps.

If the model is non-stationary, a matrix of quantiles is returned, where each row
corresponds to a MCMC step and each column to a covariate.

"""
function quantile(fm::BayesianEVA,p::Real)::Array{<:Real}

    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    θ = slicematrix(fm.sim.value[:,:,1], dims=2)

    q = quantile.(fm.model, θ, p)

    if !(typeof(q) <: Vector{<:Real})
        q = unslicematrix(q, dims=2)
    end

    return q

end

"""
    returnlevel(fm::BayesianEVA{BlockMaxima}, returnPeriod::Real)::ReturnLevel

Compute the return level corresponding to the return period `returnPeriod` from the fitted model `fm`.

"""
function returnlevel(fm::BayesianEVA{BlockMaxima}, returnPeriod::Real)::ReturnLevel

      @assert returnPeriod > zero(returnPeriod) "the return period should be positive."

      # quantile level
      p = 1-1/returnPeriod

      Q = quantile(fm, p)

      return ReturnLevel(BlockMaximaModel(fm), returnPeriod, Q)

end

function cint(rl::ReturnLevel{BayesianEVA{BlockMaxima}}, confidencelevel::Real=.95)::Vector{Vector{Real}}

      @assert rl.returnperiod > zero(rl.returnperiod) "the return period should be positive."
      @assert zero(confidencelevel)<confidencelevel<one(confidencelevel) "the confidence level should be in (0,1)."

      # quantile level
      p = 1-1/rl.returnperiod

      α = (1 - confidencelevel)

      Q = Chains(rl.value)

      ci = Mamba.hpd(Q, alpha = α)

      return slicematrix(ci.value[:,:,1], dims=2)

end

"""
    returnlevel(fm::BayesianEVA{ThresholdExceedance}, threshold::Real, nobservation::Int,
        nobsperblock::Int, returnPeriod::Real)::ReturnLevel

Compute the return level corresponding to the return period `returnPeriod` from the fitted model `fm`.

The threshold should be a scalar. A varying threshold is not yet implemented.

"""
function returnlevel(fm::BayesianEVA{ThresholdExceedance}, threshold::Real, nobservation::Int,
    nobsperblock::Int, returnPeriod::Real)::ReturnLevel

    @assert returnPeriod > zero(returnPeriod) "the return period should be positive."

    # Exceedance probability
    ζ = length(fm.model.data.value)/nobservation

    # Appropriate quantile level given the probability exceedance and the number of obs per year
    p = 1-1/(returnPeriod * nobsperblock * ζ)

    Q = quantile(fm, p)

    return ReturnLevel(PeakOverThreshold(fm, threshold, nobservation, nobsperblock),
        returnPeriod, threshold .+ Q)

end

function cint(rl::ReturnLevel{BayesianEVA{ThresholdExceedance}}, confidencelevel::Real=.95)::Vector{Vector{Real}}

    @assert rl.returnperiod > zero(rl.returnperiod) "the return period should be positive."
    @assert zero(confidencelevel)<confidencelevel<one(confidencelevel) "the confidence level should be in (0,1)."

    α = 1 - confidencelevel

    Q = Chains(rl.value)

    ci = Mamba.hpd(Q, alpha = α)

    return slicematrix(ci.value[:,:,1], dims=2)

end

"""
    showfittedEVA(io::IO, obj::BayesianEVA; prefix::String = "")

Displays a BayesianEVA with the prefix `prefix` before every line.
"""
function showfittedEVA(io::IO, obj::BayesianEVA; prefix::String = "")

    println(io, prefix, "BayesianEVA")
    println(io, prefix, "model :")
    showEVA(io, obj.model, prefix = prefix*"\t")
    println(io)
    println(io, prefix, "sim :")
    showChain(io, obj.sim, prefix = prefix*"\t")

end

"""
    showChain(io::IO, obj::Mamba.Chains; prefix::String = "")

Displays a Mamba.Chains with the prefix `prefix` before every line.
"""
function showChain(io::IO, chain::Mamba.Chains; prefix::String = "")

    println(io, prefix, "Mamba.Chains")
    println(io, prefix, "Iterations :\t\t", chain.range[1], ":", chain.range[end])
    println(io, prefix, "Thinning interval :\t", step(chain.range))
    println(io, prefix, "Chains :\t\t", length(chain.chains))
    println(io, prefix, "Samples per chain :\t", size(chain.value, 1))
    println(io, prefix, "Value :\t\t\t", typeof(chain.value), "[", size(chain.value, 1),
        ",", size(chain.value, 2), ",", size(chain.value, 3),"]")

end

"""
    transform(fm::BayesianEVA{BlockMaxima})::BayesianEVA

Transform the fitted model for the original covariate scales.
"""
function transform(fm::BayesianEVA{BlockMaxima})::BayesianEVA

    locationcovstd = fm.model.location.covariate
    logscalecovstd = fm.model.logscale.covariate
    shapecovstd = fm.model.shape.covariate

    locationcov = Extremes.reconstruct.(locationcovstd)
    logscalecov = Extremes.reconstruct.(logscalecovstd)
    shapecov = Extremes.reconstruct.(shapecovstd)

    # Model on the original covariate scale
    model = BlockMaxima(fm.model.data, locationcov = locationcov, logscalecov = logscalecov, shapecov = shapecov)

    # Transformation of the parameter estimates
    z = fm.sim.value[:,:,1]

    x = deepcopy(z)
    ind = Extremes.paramindex(fm.model)

    for (var, par) in zip([locationcovstd, logscalecovstd, shapecovstd],[:μ, :ϕ, :ξ])
        if !isempty(var)
            a = getfield.(var, :scale)
            b = getfield.(var, :offset)

            for i=1:length(a)
                x[:,ind[par][1]] = x[:,ind[par][1]] - z[:,ind[par][1+i]] * b[i]/a[i]
                x[:,ind[par][1+i]] = z[:,ind[par][1+i]]/a[i]
            end
        end
    end

    sim = fm.sim
    sim.value[:,:,1] = x

    # Contruction of the fittedEVA structure
    return BayesianEVA(model, sim)

end

"""
    transform(fm::BayesianEVA{ThresholdExceedance})

Transform the fitted model for the original covariate scales.
"""
function transform(fm::BayesianEVA{ThresholdExceedance})

    logscalecovstd = fm.model.logscale.covariate
    shapecovstd = fm.model.shape.covariate

    logscalecov = Extremes.reconstruct.(logscalecovstd)
    shapecov = Extremes.reconstruct.(shapecovstd)

    # Model on the original covariate scale
    model = ThresholdExceedance(fm.model.data, logscalecov = logscalecov, shapecov = shapecov)

    # Transformation of the parameter estimates
    z = fm.sim.value[:,:,1]

    x = deepcopy(z)
    ind = Extremes.paramindex(fm.model)

    for (var, par) in zip([logscalecovstd, shapecovstd],[:ϕ, :ξ])
        if !isempty(var)
            a = getfield.(var, :scale)
            b = getfield.(var, :offset)

            for i=1:length(a)
                x[:,ind[par][1]] = x[:,ind[par][1]] - z[:,ind[par][1+i]] * b[i]/a[i]
                x[:,ind[par][1+i]] = z[:,ind[par][1+i]]/a[i]
            end
        end
    end

    sim = fm.sim
    sim.value[:,:,1] = x

    # Contruction of the fittedEVA structure
    return BayesianEVA(model, sim)

end

"""
    parametervar(fm::BayesianEVA)::Array{Float64, 2}

Compute the covariance parameters estimate of the fitted model `fm`.

"""
function parametervar(fm::BayesianEVA)::Array{Float64, 2}

    x = fm.sim.value[:,:,1]

    C = cov(x)

    return C

end

function cint(fm::BayesianEVA, confidencelevel::Real=.95)::Array{Array{Float64,1},1}

    @assert 0<confidencelevel<1 "the confidence level should be between 0 and 1."

    α = 1-confidencelevel

    # Chain summary
    ci = hpd(fm.sim[:,:,1], alpha = α)

    return slicematrix(ci.value[:,:,1], dims=2)

end

"""
    findposteriormode(fm::BayesianEVA)::Vector{<:Real}

Find the maximum a posteriori probability (MAP) estimate.
"""
function findposteriormode(fm::BayesianEVA)::Vector{<:Real}

    θ = Extremes.slicematrix(fm.sim.value[:,:,1], dims=2)

    ll = Extremes.loglike.(fm.model,θ)

    ind = argmax(ll)

    θ̂ = θ[ind]

    return θ̂

end


text:
function Restriction(nz::Int,nx::Int,ext::Int,igz::Array{Int,1},igx::Array{Int,1})

# This function creates the reciver restriction/sampling operator/matrix.
#
# INPUTS:     nz        - Number of grid points in the z-direction
#             nx        - Number of grid points in the x-direction
#             igz       - z-indices for receiver locations
#             igx       - x-indices for receiver locations
#
# OUTPUTS:    R         - Restriction operator

    if length(igx)!=length(igz)
        error("The lengths of 'igx' and 'igz' must be the same.")
    end

    I = zeros(Int,length(igz))
    J = zeros(Int,length(igz))
    V = zeros(Int,length(igz))
    for i = 1:length(igz)
        I[i] = i
        J[i] = (ext+igx[i]-1)*nz + ext + igz[i]
        V[i] = 1
    end

    R = sparse(I,J,V,length(igz),nz*nx)

    return R

end


text:
using ClimaCore: Geometry, Domains, Meshes, Topologies
using Test

@testset "interior faces iterator" begin
    @testset "1 element across each panel" begin
        topology = Topologies.Topology2D(
            Meshes.EquiangularCubedSphere(Domains.SphereDomain(5.0), 1),
        )
        @test length(Topologies.interior_faces(topology)) == 12
    end
    @testset "3 elements across each panel" begin
        topology = Topologies.Topology2D(
            Meshes.EquiangularCubedSphere(Domains.SphereDomain(5.0), 3),
        )
        @test length(Topologies.interior_faces(topology)) == 6 * 2 * 6 + 12 * 3
    end
end

@testset "boundaries" begin
    @testset "1 element across each panel" begin
        topology = Topologies.Topology2D(
            Meshes.EquiangularCubedSphere(Domains.SphereDomain(5.0), 1),
        )
        @test isempty(Topologies.boundary_tags(topology))
    end
    @testset "3 elements across each panel" begin
        topology = Topologies.Topology2D(
            Meshes.EquiangularCubedSphere(Domains.SphereDomain(5.0), 3),
        )
        @test isempty(Topologies.boundary_tags(topology))
    end
end

@testset "boundaries" begin
    @testset "1 element across each panel" begin
        topology = Topologies.Topology2D(
            Meshes.EquiangularCubedSphere(Domains.SphereDomain(5.0), 1),
        )
        length(Topologies.vertices(topology)) == 8
        for uvert in Topologies.vertices(topology)
            @test length(uvert) == 3
        end
    end
    @testset "3 elements across each panel" begin
        topology = Topologies.Topology2D(
            Meshes.EquiangularCubedSphere(Domains.SphereDomain(5.0), 3),
        )
        length(Topologies.vertices(topology)) == 8 + 12 * 2 + 6 * 4
        for uvert in Topologies.vertices(topology)
            @test length(uvert) in (3, 4)
        end
    end
end


text:
mutable struct OBOPParameters
    gama
    max_generations
    population_size
    elite_size
    mutation_size
    function OBOPParameters(gama::Float64)
        # Max and Min population size
        population_max::Int64 = 1000
        population_min::Int64 = 50

        # Stop criterium
        max_generations::Int64 = convert(Int64, round(gama^(log(population_min) / log(gama) - population_max), digits=0))

        # Set elite size
        percentual_elite = 0.10 + ((1.0 / max_generations) * (0.25 - 0.10))
        elite_size = convert(Int64, round(population_max * percentual_elite, digits=0))

        # Set mutation size
        percentual_mutation = 0.05 + (1.0 - (1.0 / max_generations)) * (0.20 - 0.05)
        mutation_size = convert(Int64, round(population_max * percentual_mutation, digits=0))

        new(gama, max_generations, population_max, elite_size, mutation_size)
    end
end

function update_parameters!(generation::Int64, parameters::OBOPParameters)
    # population size: cooling annealing
    parameters.population_size = convert(Int64, round(( parameters.population_size * parameters.gama), digits=0))

    # Update elite size
    percentual_elite = 0.10 + ((generation / parameters.max_generations) * (0.25 - 0.10))
    parameters.elite_size = convert(Int64, round(parameters.population_size * percentual_elite, digits=0))

    # Update mutation size
    percentual_mutation = 0.05 + (1.0 - (generation / parameters.max_generations)) * (0.20 - 0.05)
    parameters.mutation_size = convert(Int64, round(parameters.population_size * percentual_mutation, digits=0))
end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
"""
    checkmethod_plot(method)

"""
function checkmethod_plot(
        method::LearningMethod;
        model = mockmodel(method),
        sample = mocksample(method),
        devicefn = cpu,
        context = Training())
    @testset "Plotting interface" begin
        x, y = DLPipelines.encode(method, context, sample)
        ŷ = DLPipelines._predictx(method, model, x, devicefn)

        @testset "plotsample!" begin
            @test_nowarn plotsample(method, sample; resolution = (200, 200))
        end
        @testset "plotxy!" begin
            @test_nowarn plotxy(method, x, y; resolution = (200, 200))
        end
        @testset "plotprediction!" begin
            @test_nowarn plotprediction(method, x, ŷ, y; resolution = (200, 200))
        end
    end
end


text:
export make_resident, evict

## memory

function make_resident(ctx::ZeContext, dev::ZeDevice, buf::AbstractBuffer, size=sizeof(buf))
    zeContextMakeMemoryResident(ctx, dev, buf, size)
end

function evict(ctx::ZeContext, dev::ZeDevice, buf::AbstractBuffer, size=sizeof(buf))
    zeContextEvictMemory(ctx, dev, buf, size)
end


text:
using ErlPort.Exceptions
using ErlPort.ErlTerms.Encode

include("Testing.jl")

# tests for supporting functions
testcase() do
    @test charint4pack(1) == b"\x00\x00\x00\x01"
    @test charint4pack(1000) == b"\x00\x00\x03\xe8"
    @test charint4pack(72000) == b"\x00\x01\x19@"
    @test charint4pack(1000000000) == b";\x9a\xca\x00"
    @test charint4pack(4294967295) == b"\xff\xff\xff\xff"
    @test charint2pack(1) == b"\x00\x01"
    @test charint2pack(1000) == b"\x03\xe8"
    @test charint2pack(10000) == b"'\x10"
    @test charint2pack(50092) == b"\xc3\xac"
    @test charint2pack(65535) == b"\xff\xff"
end

# tests for atoms
testcase() do
    @test encode(:erlang) == b"\x83d\0\6erlang"
    @test encode(:test) == b"\x83d\0\4test"
    @test encode(symbol("")) == b"\x83d\0\0"
    @test encodeterm(:erlang) == b"d\0\6erlang"
    @test encodeterm(:test) == b"d\0\4test"
    @test encodeterm(symbol("")) == b"d\0\0"
end

# tests for bools
testcase() do
    @test encode(true) == b"\x83d\0\4true"
    @test encode(:true) == b"\x83d\0\4true"
    @test encode(false) == b"\x83d\0\5false"
    @test encode(:false) == b"\x83d\0\5false"
    @test encodeterm(true) == b"d\0\4true"
    @test encodeterm(:true) == b"d\0\4true"
    @test encodeterm(false) == b"d\0\5false"
    @test encodeterm(:false) == b"d\0\5false"
end

# tests for nothing
testcase() do
    @test encode(nothing) == b"\x83d\0\7nothing"
    @test encode(:nothing) == b"\x83d\0\7nothing"
    @test encodeterm(nothing) == b"d\0\7nothing"
    @test encodeterm(:nothing) == b"d\0\7nothing"
end

# tests for list
testcase() do
end

# tests for string
testcase() do
end

# tests for tuples
testcase() do
end

# tests for improper list
testcase() do
end

# tests for short int
testcase() do
    @test encode(0) == b"\x83a\0"
    @test encode(1) == b"\x83a\1"
    @test encode(255) == b"\x83a\xff"
end

# tests for int
testcase() do
    @test encode(-1) == b"\x83b\xff\xff\xff\xff"
    @test encode(256) == b"\x83b\0\0\1\0"
    @test encode(-2147483648) == b"\x83b\x80\0\0\0"
    @test encode(2147483647) == b"\x83b\x7f\xff\xff\xff"
end

# tests for long int
testcase() do
    @test encode(-2147483649) == b"\x83n\4\1\1\0\0\x80"
    @test encode(2147483648) == b"\x83n\4\0\0\0\0\x80"
    @test encode(BigInt(2)^2040) == vcat(b"\x83o\0\0\1\0\0", zeros(255), b"\1")
    @test encode(-BigInt(2)^2040) == vcat(b"\x83o\0\0\1\0\1", zeros(255), b"\1")
end

# tests for int errors
# XXX add some!

# tests for float
testcase() do
end

# tests for opaque objects
testcase() do
end

# tests for compressed terms
testcase() do
end