text:
"""
    LinearMixedModel

Linear mixed-effects model representation

## Fields

* `formula`: the formula for the model
* `allterms`: a vector of random-effects terms, the fixed-effects terms and the response
* `sqrtwts`: vector of square roots of the case weights.  Can be empty.
* `A`: an `nt × nt` symmetric `BlockMatrix` of matrices representing `hcat(Z,X,y)'hcat(Z,X,y)`
* `L`: a `nt × nt` `BlockMatrix` - the lower Cholesky factor of `Λ'AΛ+I`
* `optsum`: an [`OptSummary`](@ref) object

## Properties

* `θ` or `theta`: the covariance parameter vector used to form λ
* `β` or `beta`: the fixed-effects coefficient vector
* `λ` or `lambda`: a vector of lower triangular matrices repeated on the diagonal blocks of `Λ`
* `σ` or `sigma`: current value of the standard deviation of the per-observation noise
* `b`: random effects on the original scale, as a vector of matrices
* `reterms`: a `Vector{ReMat{T}}` of random-effects terms.
* `feterms`: a `Vector{FeMat{T}}` of the fixed-effects model matrix and the response
* `u`: random effects on the orthogonal scale, as a vector of matrices
* `lowerbd`: lower bounds on the elements of θ
* `X`: the fixed-effects model matrix
* `y`: the response vector
"""
struct LinearMixedModel{T<:AbstractFloat} <: MixedModel{T}
    formula::FormulaTerm
    allterms::Vector{Union{ReMat{T}, FeMat{T}}}
    sqrtwts::Vector{T}
    A::BlockMatrix{T}            # cross-product blocks
    L::BlockMatrix{T}
    optsum::OptSummary{T}
end
LinearMixedModel(f::FormulaTerm, tbl; contrasts = Dict{Symbol,Any}(), wts = []) =
    LinearMixedModel(
        f::FormulaTerm,
        Tables.columntable(tbl),
        contrasts = contrasts,
        wts = wts,
    )
function LinearMixedModel(
    f::FormulaTerm,
    tbl::Tables.ColumnTable;
    contrasts = Dict{Symbol,Any}(),
    wts = [],
)
    # TODO: perform missing_omit() after apply_schema() when improved
    # missing support is in a StatsModels release
    tbl, _ = StatsModels.missing_omit(tbl, f)
    form = apply_schema(f, schema(f, tbl, contrasts), LinearMixedModel)
    # tbl, _ = StatsModels.missing_omit(tbl, form)

    y, Xs = modelcols(form, tbl)

    y = reshape(float(y), (:, 1)) # y as a floating-point matrix
    T = eltype(y)

    reterms = ReMat{T}[]
    feterms = FeMat{T}[]
    for (i, x) in enumerate(Xs)
        if isa(x, ReMat{T})
            push!(reterms, x)
        else
            cnames = coefnames(form.rhs[i])
            push!(feterms, FeMat(x, isa(cnames, String) ? [cnames] : collect(cnames)))
        end
    end
    push!(feterms, FeMat(y, [""]))

    # detect and combine RE terms with the same grouping var
    if length(reterms) > 1
        reterms = amalgamate(reterms)
    end

    sort!(reterms, by = nranef, rev = true)

    allterms = convert(Vector{Union{ReMat{T},FeMat{T}}}, vcat(reterms, feterms))
    A, L = createAL(allterms)
    lbd = foldl(vcat, lowerbd(c) for c in reterms)
    θ = foldl(vcat, getθ(c) for c in reterms)
    optsum = OptSummary(θ, lbd, :LN_BOBYQA, ftol_rel = T(1.0e-12), ftol_abs = T(1.0e-8))
    fill!(optsum.xtol_abs, 1.0e-10)
    LinearMixedModel(form, allterms, sqrt.(convert(Vector{T}, wts)), A, L, optsum)
end

fit(
    ::Type{LinearMixedModel},
    f::FormulaTerm,
    tbl;
    wts = [],
    contrasts = Dict{Symbol,Any}(),
    verbose::Bool = false,
    REML::Bool = false,
) = fit(
    LinearMixedModel,
    f,
    Tables.columntable(tbl),
    wts = wts,
    contrasts = contrasts,
    verbose = verbose,
    REML = REML,
)

fit(
    ::Type{LinearMixedModel},
    f::FormulaTerm,
    tbl::Tables.ColumnTable;
    wts = wts,
    contrasts = contrasts,
    verbose = verbose,
    REML = REML,
) = fit!(
    LinearMixedModel(f, tbl, contrasts = contrasts, wts = wts),
    verbose = verbose,
    REML = REML,
)

fit(
    ::Type{MixedModel},
    f::FormulaTerm,
    tbl;
    wts = [],
    contrasts = Dict{Symbol,Any}(),
    verbose::Bool = false,
    REML::Bool = false,
) = fit(
    LinearMixedModel,
    f,
    tbl,
    wts = wts,
    contrasts = contrasts,
    verbose = verbose,
    REML = REML,
)

fit(
    ::Type{MixedModel},
    f::FormulaTerm,
    tbl,
    d::Normal,
    l::IdentityLink;
    wts = [],
    contrasts = Dict{Symbol,Any}(),
    verbose::Bool = false,
    REML::Bool = false,
    offset = [],
    fast::Bool = false,
    nAGQ::Integer = 1,
) = fit(
    LinearMixedModel,
    f,
    tbl,
    wts = wts,
    contrasts = contrasts,
    verbose = verbose,
    REML = REML,
)

function StatsBase.coef(m::LinearMixedModel{T}) where {T}
    piv = fetrm(m).piv
    invpermute!(fixef!(similar(piv, T), m), piv)
end

βs(m::LinearMixedModel) = NamedTuple{(Symbol.(coefnames(m))...,)}(coef(m))

function StatsBase.coefnames(m::LinearMixedModel)
    Xtrm = fetrm(m)
    invpermute!(copy(Xtrm.cnames), Xtrm.piv)
end

function StatsBase.coeftable(m::LinearMixedModel)
    co = coef(m)
    se = stderror!(similar(co), m)
    z = co ./ se
    pvalue = ccdf.(Chisq(1), abs2.(z))
    names = coefnames(m)

    CoefTable(
        hcat(co, se, z, pvalue),
        ["Estimate", "Std.Error", "z value", "P(>|z|)"],
        names,
        4, # pvalcol
        3, # teststatcol
    )
end

"""
    cond(m::MixedModel)

Return a vector of condition numbers of the λ matrices for the random-effects terms
"""
LinearAlgebra.cond(m::MixedModel) = cond.(m.λ)

"""
    condVar(m::LinearMixedModel)

Return the conditional variances matrices of the random effects.

The random effects are returned by `ranef` as a vector of length `k`,
where `k` is the number of random effects terms.  The `i`th element
is a matrix of size `vᵢ × ℓᵢ`  where `vᵢ` is the size of the
vector-valued random effects for each of the `ℓᵢ` levels of the grouping
factor.  Technically those values are the modes of the conditional
distribution of the random effects given the observed data.

This function returns an array of `k` three dimensional arrays,
where the `i`th array is of size `vᵢ × vᵢ × ℓᵢ`.  These are the
diagonal blocks from the conditional variance-covariance matrix,

    s² Λ(Λ'Z'ZΛ + I)⁻¹Λ'
"""
function condVar(m::LinearMixedModel{T}) where {T}
    retrms = m.reterms
    t1 = first(retrms)
    L11 = m.L[Block(1, 1)]
    if !isone(length(retrms)) || !isa(L11, Diagonal{T,Vector{T}})
        throw(ArgumentError("code for multiple or vector-valued r.e. not yet written"))
    end
    ll = first(t1.λ)
    Ld = L11.diag
    Array{T,3}[reshape(abs2.(ll ./ Ld) .* varest(m), (1, 1, length(Ld)))]
end

function createAL(allterms::Vector{Union{ReMat{T},FeMat{T}}}) where {T}
    k = length(allterms)
    sz = [isa(t, ReMat) ? size(t, 2) : rank(t) for t in allterms]
    A = BlockArray(undef_blocks, AbstractMatrix{T}, sz, sz)
    L = BlockArray(undef_blocks, AbstractMatrix{T}, sz, sz)
    for j = 1:k
        for i = j:k
            Lij = L[Block(i, j)] = densify(allterms[i]' * allterms[j])
            A[Block(i, j)] = deepcopy(isa(Lij, BlockedSparse) ? Lij.cscmat : Lij)
        end
    end
    nretrm = sum(Base.Fix2(isa, ReMat), allterms)
    for i = 2:nretrm      # check for fill-in due to non-nested grouping factors
        ci = allterms[i]
        for j = 1:(i-1)
            cj = allterms[j]
            if !isnested(cj, ci)
                for l = i:k
                    L[Block(l, i)] = Matrix(L[Block(l, i)])
                end
                break
            end
        end
    end
    A, L
end

StatsBase.deviance(m::MixedModel) = objective(m)

GLM.dispersion(m::LinearMixedModel, sqr::Bool = false) = sqr ? varest(m) : sdest(m)

GLM.dispersion_parameter(m::LinearMixedModel) = true

StatsBase.dof(m::LinearMixedModel) = size(m)[2] + nθ(m) + 1

function StatsBase.dof_residual(m::LinearMixedModel)::Int
    # nobs - rank(FE) - 1 (dispersion)
    # this differs from lme4 by not including nθ
    # a better estimate would be a number somewhere between the number of
    # variance components and the number of conditional modes
    # nobs, rank FE, num conditional modes, num grouping vars
    nobs(m) - size(m)[2] - 1
end

"""
    feind(m::MixedModel)

An internal utility to return the index in `m.allterms` of the fixed-effects term.
"""
feind(m::MixedModel) = findfirst(Base.Fix2(isa, FeMat), m.allterms)

"""
    feL(m::MixedModel)

Return the lower Cholesky factor for the fixed-effects parameters, as an `LowerTriangular`
`p × p` matrix.
"""
function feL(m::LinearMixedModel)
    k = feind(m)
    LowerTriangular(m.L.blocks[k, k])
end

"""
    fetrm(m::LinearMixedModel)

Return the fixed-effects term from `m.allterms`
"""
fetrm(m) = m.allterms[feind(m)]

"""
    fit!(m::LinearMixedModel[; verbose::Bool=false, REML::Bool=false])

Optimize the objective of a `LinearMixedModel`.  When `verbose` is `true` the values of the
objective and the parameters are printed on stdout at each function evaluation.
"""
function fit!(m::LinearMixedModel{T}; verbose::Bool = false, REML::Bool = false) where {T}
    optsum = m.optsum
    opt = Opt(optsum)
    feval = 0
    optsum.REML = REML
    function obj(x, g)
        isempty(g) || error("gradient not defined")
        feval += 1
        val = objective(updateL!(setθ!(m, x)))
        feval == 1 && (optsum.finitial = val)
        verbose && println("f_", feval, ": ", round(val, digits = 5), " ", x)
        val
    end
    NLopt.min_objective!(opt, obj)
    fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial))
    ## check if small non-negative parameter values can be set to zero
    xmin_ = copy(xmin)
    lb = optsum.lowerbd
    for i in eachindex(xmin_)
        if iszero(lb[i]) && zero(T) < xmin_[i] < T(0.001)
            xmin_[i] = zero(T)
        end
    end
    if xmin_ ≠ xmin
        if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
            fmin = zeroobj
            copyto!(xmin, xmin_)
        end
    end
    ## ensure that the parameter values saved in m are xmin
    updateL!(setθ!(m, xmin))

    optsum.feval = feval
    optsum.final = xmin
    optsum.fmin = fmin
    optsum.returnvalue = ret
    ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited")
    if ret ∈ [:FAILURE, :INVALID_ARGS, :OUT_OF_MEMORY, :FORCED_STOP, :MAXFEVAL_REACHED]
        @warn("NLopt optimization failure: $ret")
    end
    m
end

function fitted!(v::AbstractArray{T}, m::LinearMixedModel{T}) where {T}
    ## FIXME: Create and use `effects(m) -> β, b` w/o calculating β twice
    Xtrm = fetrm(m)
    vv = mul!(vec(v), Xtrm, fixef!(similar(Xtrm.piv, T), m))
    for (rt, bb) in zip(m.reterms, ranef(m))
        unscaledre!(vv, rt, bb)
    end
    v
end

StatsBase.fitted(m::LinearMixedModel{T}) where {T} = fitted!(Vector{T}(undef, nobs(m)), m)

"""
    fixef!(v::Vector{T}, m::LinearMixedModel{T})

Overwrite `v` with the pivoted fixed-effects coefficients of model `m`

For full-rank models the length of `v` must be the rank of `X`.  For rank-deficient models
the length of `v` can be the rank of `X` or the number of columns of `X`.  In the latter
case the calculated coefficients are padded with -0.0 out to the number of columns.
"""
function fixef!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
    Xtrm = fetrm(m)
    if isfullrank(Xtrm)
        ldiv!(feL(m)', copyto!(v, m.L.blocks[end, end-1]))
    else
        ldiv!(
            feL(m)',
            view(copyto!(fill!(v, -zero(T)), m.L.blocks[end, end-1]), 1:(Xtrm.rank)),
        )
    end
    v
end

"""
    fixef(m::MixedModel)

Return the fixed-effects parameter vector estimate of `m`.

In the rank-deficient case the truncated parameter vector, of length `rank(m)` is returned.
This is unlike `coef` which always returns a vector whose length matches the number of
columns in `X`.
"""
fixef(m::LinearMixedModel{T}) where {T} = fixef!(Vector{T}(undef, fetrm(m).rank), m)

"""
    fixefnames(m::MixedModel)

Return a (permuted and truncated in the rank-deficient case) vector of coefficient names.
"""
function fixefnames(m::LinearMixedModel{T}) where {T}
    Xtrm = fetrm(m)
    Xtrm.cnames[1:Xtrm.rank]
end

"""
    fnames(m::MixedModel)

Return the names of the grouping factors for the random-effects terms.
"""
fnames(m::MixedModel) = ((tr.trm.sym for tr in m.reterms)...,)

"""
    getθ(m::LinearMixedModel)

Return the current covariance parameter vector.
"""
getθ(m::LinearMixedModel{T}) where {T} = foldl(vcat, getθ.(m.reterms))

function getθ!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
    k = 0
    for t in m.reterms
        nt = nθ(t)
        getθ!(view(v, (k+1):(k+nt)), t)
        k += nt
    end
    v
end

function Base.getproperty(m::LinearMixedModel{T}, s::Symbol) where {T}
    if s == :θ || s == :theta
        getθ(m)
    elseif s == :β || s == :beta
        coef(m)
    elseif s == :βs || s == :betas
        βs(m)
    elseif s == :λ || s == :lambda
        getproperty.(m.reterms, :λ)
    elseif s == :σ || s == :sigma
        sdest(m)
    elseif s == :σs || s == :sigmas
        σs(m)
    elseif s == :σρs || s == :sigmarhos
        σρs(m)
    elseif s == :b
        ranef(m)
    elseif s == :feterms
        convert(Vector{FeMat{T}}, filter(Base.Fix2(isa, FeMat), getfield(m, :allterms)))
    elseif s == :objective
        objective(m)
    elseif s == :corr
        vcov(m, corr=true)
    elseif s == :vcov
        vcov(m, corr=false)
    elseif s == :PCA
        NamedTuple{fnames(m)}(PCA.(m.reterms))
    elseif s == :pvalues
        ccdf.(Chisq(1), abs2.(coef(m) ./ stderror(m)))
    elseif s == :reterms
        convert(Vector{ReMat{T}}, filter(Base.Fix2(isa, ReMat), getfield(m, :allterms)))
    elseif s == :stderror
        stderror(m)
    elseif s == :u
        ranef(m, uscale = true)
    elseif s == :lowerbd
        m.optsum.lowerbd
    elseif s == :X
        modelmatrix(m)
    elseif s == :y
        vec(last(m.feterms).x)
    elseif s == :rePCA
        rePCA(m)
    else
        getfield(m, s)
    end
end

"""
    issingular(m::LinearMixedModel, θ=m.θ)

Test whether the model `m` is singular if the parameter vector is `θ`.
"""
issingular(m::LinearMixedModel, θ=m.θ) = any(isapprox.(lowerbd(m), θ))

function StatsBase.leverage(m::LinearMixedModel{T}) where {T}
    # This can be done more efficiently but reusing existing tools is easier.
    # The i'th leverage value is obtained by replacing the response with the i'th
    # basis vector, updating A and L, then taking the sum of squared values of the
    # last row of L, excluding the last position.
    yorig = copy(m.y)
    l = length(m.allterms)
    value = map(eachindex(yorig)) do i
        fill!(m.y, zero(T))
        m.y[i] = one(T)
        reevaluateAend!(m)
        updateL!(m)
        sum(j -> sum(abs2, m.L[Block(l, j)]), 1:(l-1))
    end
    copyto!(m.y, yorig)
    updateL!(reevaluateAend!(m))
    value
end

function StatsBase.loglikelihood(m::LinearMixedModel)
    if m.optsum.REML
        throw(ArgumentError("loglikelihood not available for models fit by REML"))
    end
    -objective(m) / 2
end

lowerbd(m::LinearMixedModel) = m.optsum.lowerbd

function StatsBase.modelmatrix(m::LinearMixedModel)
    fetrm = first(m.feterms)
    if fetrm.rank == size(fetrm, 2)
        fetrm.x
    else
        fetrm.x[:, invperm(fetrm.piv)]
    end
end

nθ(m::LinearMixedModel) = sum(nθ, m.allterms)

StatsBase.nobs(m::LinearMixedModel) = first(size(m))

"""
    objective(m::LinearMixedModel)

Return negative twice the log-likelihood of model `m`
"""
function objective(m::LinearMixedModel)
    wts = m.sqrtwts
    logdet(m) + ssqdenom(m) * (1 + log2π + log(varest(m))) -
    (isempty(wts) ? 0 : 2 * sum(log, wts))
end

StatsBase.predict(m::LinearMixedModel) = fitted(m)

Base.propertynames(m::LinearMixedModel, private = false) = (
    :formula,
    :sqrtwts,
    :A,
    :L,
    :optsum,
    :θ,
    :theta,
    :β,
    :beta,
    :λ,
    :lambda,
    :stderror,
    :σ,
    :sigma,
    :σs,
    :sigmas,
    :b,
    :u,
    :lowerbd,
    :X,
    :y,
    :corr,
    :vcov,
    :PCA,
    :rePCA,
    :reterms,
    :feterms,
    :objective,
    :pvalues,
)

"""
    pwrss(m::LinearMixedModel)

The penalized, weighted residual sum-of-squares.
"""
pwrss(m::LinearMixedModel) = abs2(sqrtpwrss(m))

"""
    ranef!(v::Vector{Matrix{T}}, m::MixedModel{T}, β, uscale::Bool) where {T}

Overwrite `v` with the conditional modes of the random effects for `m`.

If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise
on the original scale
"""
function ranef!(
    v::Vector,
    m::LinearMixedModel{T},
    β::AbstractArray{T},
    uscale::Bool,
) where {T}
    (k = length(v)) == length(m.reterms) || throw(DimensionMismatch(""))
    L = m.L
    for j = 1:k
        mul!(
            vec(copyto!(v[j], L[Block(length(m.allterms), j)])),
            L[Block(k + 1, j)]',
            β,
            -one(T),
            one(T),
        )
    end
    for i = k:-1:1
        Lii = L[Block(i, i)]
        vi = vec(v[i])
        ldiv!(adjoint(isa(Lii, Diagonal) ? Lii : LowerTriangular(Lii)), vi)
        for j = 1:(i-1)
            mul!(vec(v[j]), L[Block(i, j)]', vi, -one(T), one(T))
        end
    end
    if !uscale
        for (t, vv) in zip(m.reterms, v)
            lmul!(t.λ, vv)
        end
    end
    v
end

ranef!(v::Vector, m::LinearMixedModel, uscale::Bool) = ranef!(v, m, fixef(m), uscale)

"""
    ranef(m::LinearMixedModel; uscale=false, named=true)

Return, as a `Vector{Vector{T}}` (`Vector{NamedVector{T}}` if `named=true`),
the conditional modes of the random effects in model `m`.

If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise on
the original scale.
"""
function ranef(m::LinearMixedModel{T}; uscale = false, named = false) where {T}
    reterms = m.reterms
    v = [Matrix{T}(undef, size(t.z, 1), nlevs(t)) for t in reterms]
    ranef!(v, m, uscale)
    named || return v
    vnmd = map(NamedArray, v)
    for (trm, vnm) in zip(reterms, vnmd)
        setnames!(vnm, trm.cnames, 1)
        setnames!(vnm, string.(trm.trm.contrasts.levels), 2)
    end
    vnmd
end

LinearAlgebra.rank(m::LinearMixedModel) = first(m.feterms).rank

"""
    rePCA(m::LinearMixedModel; corr::Bool=true)

Return a named tuple of the normalized cumulative variance of a principal components
analysis of the random effects covariance matrices or correlation
matrices when `corr` is `true`.

The normalized cumulative variance is the proportion of the variance for the first
principal component, the first two principal components, etc.  The last element is
always 1.0 representing the complete proportion of the variance.
"""
function rePCA(m::LinearMixedModel; corr::Bool=true)
    pca = PCA.(m.reterms, corr=corr)
    NamedTuple{fnames(m)}(getproperty.(pca,:cumvar))
end

"""
    PCA(m::LinearMixedModel; corr::Bool=true)

Return a named tuple of the principal components analysis of the random effects
covariance matrices or correlation matrices when `corr` is `true`.
"""

function PCA(m::LinearMixedModel; corr::Bool=true)
    NamedTuple{fnames(m)}(PCA.(m.reterms, corr=corr))
end

"""
    reevaluateAend!(m::LinearMixedModel)

Reevaluate the last column of `m.A` from `m.feterms`.  This function should be called
after updating the response, `m.feterms[end]`.
"""
function reevaluateAend!(m::LinearMixedModel)
    A = m.A
    trmn = reweight!(last(m.allterms), m.sqrtwts)
    nblk = length(m.allterms)
    for (j, trm) in enumerate(m.allterms)
        mul!(A[Block(nblk, j)], trmn', trm)
    end
    m
end

"""
    refit!(m::LinearMixedModel[, y::Vector])

Refit the model `m` after installing response `y`.

If `y` is omitted the current response vector is used.
"""
refit!(m::LinearMixedModel) = fit!(reevaluateAend!(m))

function refit!(m::LinearMixedModel, y)
    resp = last(m.feterms)
    length(y) == size(resp, 1) || throw(DimensionMismatch(""))
    copyto!(resp, y)
    refit!(m)
end

StatsBase.residuals(m::LinearMixedModel) = response(m) .- fitted(m)

StatsBase.response(m::LinearMixedModel) = vec(last(m.feterms).x)

function reweight!(m::LinearMixedModel, weights)
    sqrtwts = map!(sqrt, m.sqrtwts, weights)
    reweight!.(m.allterms, Ref(sqrtwts))
    updateA!(m)
    updateL!(m)
end

"""
    sdest(m::LinearMixedModel)

Return the estimate of σ, the standard deviation of the per-observation noise.
"""
sdest(m::LinearMixedModel) = √varest(m)

"""
    setθ!(m::LinearMixedModel, v)

Install `v` as the θ parameters in `m`.
"""
function setθ!(m::LinearMixedModel, v)
    offset = 0
    for trm in m.reterms
        k = nθ(trm)
        setθ!(trm, view(v, (1:k) .+ offset))
        offset += k
    end
    m
end

Base.setproperty!(m::LinearMixedModel, s::Symbol, y) =
    s == :θ ? setθ!(m, y) : setfield!(m, s, y)

function Base.show(io::IO, m::LinearMixedModel)
    if m.optsum.feval < 0
        @warn("Model has not been fit")
        return nothing
    end
    n, p, q, k = size(m)
    REML = m.optsum.REML
    println(io, "Linear mixed model fit by ", REML ? "REML" : "maximum likelihood")
    println(io, " ", m.formula)
    oo = objective(m)
    if REML
        println(io, " REML criterion at convergence: ", oo)
    else
        nums = showoff([-oo / 2, oo, aic(m), bic(m)])
        fieldwd = max(maximum(textwidth.(nums)) + 1, 11)
        for label in [" logLik", "-2 logLik", "AIC", "BIC"]
            print(io, rpad(lpad(label, (fieldwd + textwidth(label)) >> 1), fieldwd))
        end
        println(io)
        print.(Ref(io), lpad.(nums, fieldwd))
        println(io)
    end
    println(io)

    show(io, VarCorr(m))

    print(io, " Number of obs: $n; levels of grouping factors: ")
    join(io, nlevs.(m.reterms), ", ")
    println(io)
    println(io, "\n  Fixed-effects parameters:")
    show(io, coeftable(m))
end

function σs(m::LinearMixedModel)
    σ = sdest(m)
    NamedTuple{fnames(m)}(((σs(t, σ) for t in m.reterms)...,))
end

function σρs(m::LinearMixedModel)
    σ = sdest(m)
    NamedTuple{fnames(m)}(((σρs(t, σ) for t in m.reterms)...,))
end

"""
    size(m::LinearMixedModel)

Returns the size of a mixed model as a tuple of length four:
the number of observations, the number of (non-singular) fixed-effects parameters,
the number of conditional modes (random effects), the number of grouping variables
"""
function Base.size(m::LinearMixedModel)
    n, p = size(fetrm(m))
    reterms = m.reterms
    n, p, sum(size.(reterms, 2)), length(reterms)
end

"""
    sqrtpwrss(m::LinearMixedModel)

Return the square root of the penalized, weighted residual sum-of-squares (pwrss).

This value is the contents of the `1 × 1` bottom right block of `m.L`
"""
sqrtpwrss(m::LinearMixedModel) = first(m.L.blocks[end, end])

"""
    ssqdenom(m::LinearMixedModel)

Return the denominator for penalized sums-of-squares.

For MLE, this value is either the number of observation for ML. For REML, this
value is the number of observations minus the rank of the fixed-effects matrix.
The difference is analagous to the use of n or n-1 in the denominator when
calculating the variance.
"""
function ssqdenom(m::LinearMixedModel)::Int
    nobs(m) - m.optsum.REML * first(m.feterms).rank
end

"""
    std(m::MixedModel)

Return the estimated standard deviations of the random effects as a `Vector{Vector{T}}`.
"""
function Statistics.std(m::LinearMixedModel)
    rl = rowlengths.(m.reterms)
    s = sdest(m)
    isfinite(s) ? rmul!(push!(rl, [1.0]), s) : rl
end

"""
    stderror!(v::AbstractVector, m::LinearMixedModel)

Overwrite `v` with the standard errors of the fixed-effects coefficients in `m`

The length of `v` should be the total number of coefficients (i.e. `length(coef(m))`).
When the model matrix is rank-deficient the coefficients forced to `-0.0` have an
undefined (i.e. `NaN`) standard error.
"""
function stderror!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
    L = feL(m)
    scr = Vector{T}(undef, size(L, 2))
    s = sdest(m)
    fill!(v, zero(T) / zero(T))  # initialize to appropriate NaN for rank-deficient case
    for i in eachindex(scr)
        fill!(scr, false)
        scr[i] = true
        v[i] = s * norm(ldiv!(L, scr))
    end
    invpermute!(v, fetrm(m).piv)
    v
end

function StatsBase.stderror(m::LinearMixedModel{T}) where {T}
    stderror!(similar(fetrm(m).piv, T), m)
end

"""
    updateA!(m::LinearMixedModel)

Update the cross-product array, `m.A`, from `m.allterms`

This is usually done after a reweight! operation.
"""
function updateA!(m::LinearMixedModel)
    allterms = m.allterms
    k = length(allterms)
    A = m.A
    for j = 1:k
        for i = j:k
            mul!(A[Block(i, j)], allterms[i]', allterms[j])
        end
    end
    m
end

"""
    updateL!(m::LinearMixedModel)

Update the blocked lower Cholesky factor, `m.L`, from `m.A` and `m.reterms` (used for λ only)

This is the crucial step in evaluating the objective, given a new parameter value.
"""
function updateL!(m::LinearMixedModel{T}) where {T}
    A = m.A
    L = m.L
    k = length(m.allterms)
    for j = 1:k                         # copy lower triangle of A to L
        for i = j:k
            copyto!(L[Block(i, j)], A[Block(i, j)])
        end
    end
    for (j, cj) in enumerate(m.reterms)  # pre- and post-multiply by Λ, add I to diagonal
        scaleinflate!(L[Block(j, j)], cj)
        for i = (j+1):k         # postmultiply column by Λ
            rmulΛ!(L[Block(i, j)], cj)
        end
        for jj = 1:(j-1)        # premultiply row by Λ'
            lmulΛ!(cj', L[Block(j, jj)])
        end
    end
    for j = 1:k                         # blocked Cholesky
        Ljj = L[Block(j, j)]
        LjjH = isa(Ljj, Diagonal) ? Ljj : Hermitian(Ljj, :L)
        for jj = 1:(j-1)
            rankUpdate!(LjjH, L[Block(j, jj)], -one(T))
        end
        cholUnblocked!(Ljj, Val{:L})
        LjjT = isa(Ljj, Diagonal) ? Ljj : LowerTriangular(Ljj)
        for i = (j+1):k
            Lij = L[Block(i, j)]
            for jj = 1:(j-1)
                mul!(Lij, L[Block(i, jj)], L[Block(j, jj)]', -one(T), one(T))
            end
            rdiv!(Lij, LjjT')
        end
    end
    m
end

"""
    varest(m::LinearMixedModel)

Returns the estimate of σ², the variance of the conditional distribution of Y given B.
"""
varest(m::LinearMixedModel) = pwrss(m) / ssqdenom(m)

"""
    vcov(m::LinearMixedModel)

Returns the variance-covariance matrix of the fixed effects.
If `corr=true`, then correlation of fixed effects is returned instead.
"""
function StatsBase.vcov(m::LinearMixedModel{T}; corr=false) where {T}
    Xtrm = fetrm(m)
    iperm = invperm(Xtrm.piv)
    p = length(iperm)
    r = Xtrm.rank
    Linv = inv(feL(m))
    permvcov = varest(m) * (Linv'Linv)
    if p == Xtrm.rank
        vv = permvcov[iperm, iperm]
    else
        covmat = fill(zero(T) / zero(T), (p, p))
        for j = 1:r, i = 1:r
            covmat[i, j] = permvcov[i, j]
        end
        vv = covmat[iperm, iperm]
    end

    corr ?  StatsBase.cov2cor!(vv, stderror(m)) : vv
end

"""
    zerocorr!(m::LinearMixedModel[, trmnms::Vector{Symbol}])

Rewrite the random effects specification for the grouping factors in `trmnms` to zero correlation parameter.

The default for `trmnms` is all the names of random-effects terms.

A random effects term is in the zero correlation parameter configuration when the off-diagonal elements of
λ are all zero - hence there are no correlation parameters in that term being estimated.
"""
function zerocorr!(m::LinearMixedModel{T}, trmns) where {T}
    reterms = m.reterms
    for trm in reterms
        if fname(trm) in trmns
            zerocorr!(trm)
        end
    end
    optsum = m.optsum
    optsum.lowerbd = foldl(vcat, lowerbd(c) for c in reterms)
    optsum.initial = foldl(vcat, getθ(c) for c in reterms)
    optsum.final = copy(optsum.initial)
    optsum.xtol_abs = fill!(copy(optsum.initial), 1.0e-10)
    optsum.initial_step = T[]

    # the model is no longer fitted
    optsum.feval = -1

    m
end

zerocorr!(m::LinearMixedModel) = zerocorr!(m, fnames(m))


text:
for N in [Float64, Float32]
    B = BallInf(zeros(N, 2), N(1))
    E = HPolyhedron([HalfSpace(N[1], N(0)), HalfSpace(N[-1], N(-1))])  # empty
    U = Universe{N}(3)
    ε = N(1//10)
    p = N(2)

    # constructors
    X = Bloating(B, ε, p)
    X⁻ = Bloating(B, -ε, p)
    Y = Bloating(E, ε, p)
    Z = Bloating(U, ε, p)
    @test_throws AssertionError Bloating(B, ε, N(9//10))

    # dimension
    @test dim(X) == dim(X⁻) == 2 && dim(Y) == 1 && dim(Z) == 3

    # isbounded
    @test isbounded(X) && isbounded(X⁻) && isbounded(Y) && !isbounded(Z) &&
        isboundedtype(typeof(X)) && !isboundedtype(typeof(Z))

    # isempty
    @test !isempty(X) && !isempty(X⁻) && isempty(Y) && !isempty(Z)

    # an_element
    for S ∈ [X, X⁻, Z]
        v = an_element(S)
        @test v isa AbstractVector{N} && length(v) == dim(S)
    end
    @test_throws ErrorException an_element(Y)

    # tests for different norms
    for p in N[1, 2, Inf]
        B = BallInf(zeros(N, 2), N(1))
        X = Bloating(B, ε, p)
        X⁻ = Bloating(B, -ε, p)
        E = HPolyhedron([HalfSpace(N[1], N(0)), HalfSpace(N[-1], N(-1))])
        Y = Bloating(E, ε, p)  # empty set
        Bp = Ballp(p, zeros(N, 2), ε)

        # support vector and support function
        d = N[1, 0]
        @test σ(d, X) == σ(d, B + Bp)
        @test ρ(d, X) == ρ(d, B + Bp) == N(1 + ε)
        @test ρ(d, X⁻) == N(1 - ε)
        if p != N(1)
            @test ρ(d, X⁻) == dot(d, σ(d, X⁻))
        end
        @test_throws ErrorException σ(N[1], Y)
        @test_throws ErrorException ρ(N[1], Y)
        d = N[1, 99//100]
        @test σ(d, X) == σ(d, B + Bp)
        @test ρ(d, X) == ρ(d, B + Bp)

        # subset
        if p == N(Inf)
            @test X⁻ ⊆ B ⊆ X
        elseif p == N(1)
            @test B ⊆ X
        end

        # bloating of negative bloating gives the identity
        if p != N(1)
            # need to overapproximate because constraints_list is missing
            B2 = box_approximation(Bloating(Bloating(B, ε, p), -ε, p))
            @test isequivalent(B, B2)
        end

        # tests for infinity norm
        if p == Inf
            B_bloated = BallInf(zeros(N, 2), N(11/10))
            @test constraints_list(X) == constraints_list(B_bloated)
        end
    end
end

ext: jl


text:
####
#### Viscous fluxes
####

@inline viscous_flux_ux(i, j, k, grid, u, ν_ccc) = ν_ccc * ℑxᶜᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * δxᶜᵃᵃ(i, j, k, grid, u)
@inline viscous_flux_uy(i, j, k, grid, u, ν_ffc) = ν_ffc * ℑyᵃᶠᵃ(i, j, k, grid, Ayᵃᵃᶜ) * δyᵃᶠᵃ(i, j, k, grid, u)
@inline viscous_flux_uz(i, j, k, grid, u, ν_fcf) = ν_fcf * ℑzᵃᵃᶠ(i, j, k, grid, Azᵃᵃᵃ)  * δzᵃᵃᶠ(i, j, k, grid, u)

@inline viscous_flux_vx(i, j, k, grid, v, ν_ffc) = ν_ffc * ℑxᶠᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * δxᶠᵃᵃ(i, j, k, grid, v)
@inline viscous_flux_vy(i, j, k, grid, v, ν_ccc) = ν_ccc * ℑyᵃᶜᵃ(i, j, k, grid, Ayᵃᵃᶜ) * δyᵃᶜᵃ(i, j, k, grid, v)
@inline viscous_flux_vz(i, j, k, grid, v, ν_cff) = ν_cff * ℑzᵃᵃᶠ(i, j, k, grid, Azᵃᵃᵃ)  * δzᵃᵃᶠ(i, j, k, grid, v)

@inline viscous_flux_wx(i, j, k, grid, w, ν_fcf) = ν_fcf * ℑxᶠᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * δxᶠᵃᵃ(i, j, k, grid, w)
@inline viscous_flux_wy(i, j, k, grid, w, ν_cff) = ν_cff * ℑyᵃᶠᵃ(i, j, k, grid, Ayᵃᵃᶜ) * δyᵃᶠᵃ(i, j, k, grid, w)
@inline viscous_flux_wz(i, j, k, grid, w, ν_ccc) = ν_ccc * ℑzᵃᵃᶜ(i, j, k, grid, Azᵃᵃᵃ)  * δzᵃᵃᶜ(i, j, k, grid, w)

####
#### Viscous dissipation operators
####

"""
    ∇ν∇u(i, j, k, grid::Grid, u::AbstractArray, ν::AbstractFloat)

Calculates viscous dissipation for the u-velocity via

    1/V * [δxᶠᵃᵃ(ν * ℑxᶜᵃᵃ(Ax) * δxᶜᵃᵃ(u)) + δyᵃᶜᵃ(ν * ℑyᵃᶠᵃ(Ay) * δyᵃᶠᵃ(u)) + δzᵃᵃᶜ(ν * ℑzᵃᵃᶠ(Az) * δzᵃᵃᶠ(u))]

which will end up at the location `fcc`.
"""
@inline function ∇ν∇u(i, j, k, grid, ν, u)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, viscous_flux_ux, u, ν) +
                                    δyᵃᶜᵃ(i, j, k, grid, viscous_flux_uy, u, ν) +
                                    δzᵃᵃᶜ(i, j, k, grid, viscous_flux_uz, u, ν))
end

"""
    ∇ν∇v(i, j, k, grid::Grid, v::AbstractArray, ν::AbstractFloat)

Calculates viscous dissipation for the v-velocity via

    1/Vᵛ * [δxᶜᵃᵃ(ν * ℑxᶠᵃᵃ(Ax) * δxᶠᵃᵃ(v)) + δyᵃᶠᵃ(ν * ℑyᵃᶜᵃ(Ay) * δyᵃᶜᵃ(v)) + δzᵃᵃᶜ(ν * ℑzᵃᵃᶠ(Az) * δzᵃᵃᶠ(v))]

which will end up at the location `cfc`.
"""
@inline function ∇ν∇v(i, j, k, grid, ν, v)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, viscous_flux_vx, v, ν) +
                                    δyᵃᶠᵃ(i, j, k, grid, viscous_flux_vy, v, ν) +
                                    δzᵃᵃᶜ(i, j, k, grid, viscous_flux_vz, v, ν))
end

"""
    ∇ν∇w(i, j, k, grid::Grid, w::AbstractArray, ν::AbstractFloat)

Calculates viscous dissipation for the w-velocity via

    1/Vʷ * [δxᶜᵃᵃ(ν * ℑxᶠᵃᵃ(Ax) * δxᶠᵃᵃ(w)) + δyᵃᶜᵃ(ν * ℑyᵃᶠᵃ(Ay) * δyᵃᶠᵃ(w)) + δzᵃᵃᶠ(ν * ℑzᵃᵃᶜ(Az) * δzᵃᵃᶜ(w))]

which will end up at the location `ccf`.
"""
@inline function ∇ν∇w(i, j, k, grid, ν, w)
    return 1/Vᵃᵃᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, viscous_flux_wx, w, ν) +
                                    δyᵃᶜᵃ(i, j, k, grid, viscous_flux_wy, w, ν) +
                                    δzᵃᵃᶠ(i, j, k, grid, viscous_flux_wz, w, ν))
end


text:
### A Pluto.jl notebook ###
# v0.11.10

using Markdown
using InteractiveUtils

# ╔═╡ b16a96ca-eba3-11ea-3a74-4fecb5dda838
begin
	using FileIO: load
	import LibSndFile
	using Plots
	using SampledSignals
	using FFTW
end

# ╔═╡ cf3ddc9c-ebab-11ea-1417-b7eb9e517b50
md"# Aligning video recordings with Julia"

# ╔═╡ 2c57c0fe-ed5b-11ea-2837-81321cf21495
md"Due to the restrictions imposed by COVID-19, dance teachers around the world are taking their classes online. As a dance student, it would be helpful to watch yourself against a recording of your teacher.

Something like this:"

# ╔═╡ 34a5c400-ebd9-11ea-04cf-3118b3902453
html"<iframe width=\"680\" height=\"383\" src=\"https://www.youtube-nocookie.com/embed/I8kuhv0FcTE\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>"

# ╔═╡ ae29c670-ebd9-11ea-3995-77063c4c3f4d
md"It is very unlikely the song in both videos starts at exactly the same. Manually syncing them is not easy.

On the other hand, computers can compare the sound waves of each video file and quickly determine how long to wait before playing one video or the other so they are both aligned."

# ╔═╡ eb4bd808-ebc9-11ea-3c71-e9850bb24137
md"## Sound waves"

# ╔═╡ f7a9d13a-ebac-11ea-2abd-79345ce6fbef
md"[Sound waves](https://pudding.cool/2018/02/waveforms/) are represented in digital audio by very rapidly sampling the distortions of the medium (e.g. air) through a microphone and storing the resulting data in a vector."

# ╔═╡ cfdabcf8-ebae-11ea-307a-7728476c9825
md"Jazmine, the teacher, sent Julia a [recording of a jazz song](https://freesound.org/s/37750/) ([cc0](https://creativecommons.org/publicdomain/zero/1.0/)), sampled 16,000 times per second:"

# ╔═╡ 49b7f3f4-eba7-11ea-14f1-31d09ab786bb
x = load("../data/jazz.ogg")

# ╔═╡ 0342bfec-ebc2-11ea-15a2-c769c8039d17
x_duration_seconds = length(x) / x.samplerate

# ╔═╡ c10ff5da-ebae-11ea-3205-6d1f6666ca5c
typeof(x), x.samplerate, length(x)

# ╔═╡ 6e4d19da-ebaf-11ea-113d-55eca1da0769
md"Each sample in $$x$$ corresponds to a microphone reading. The first 5 samples at $$1s$$ are:"

# ╔═╡ c7eccb4a-ebad-11ea-16f5-21299446a8a2
x[16_000:16_005, 1]

# ╔═╡ becebb5c-ebaf-11ea-3e5d-4f91c13b24a1
md"Plotting the audio wave:"

# ╔═╡ 3ab06516-eba5-11ea-2316-732fe86075f6
plot(domain(x), x[:, 1], xlabel="time (s)", ylabel="amplitude", title="Jazz song", linealpha=0.85, legend=false, linecolor=palette(:default)[2], fmt = :png, dpi=300)

# ╔═╡ c0e56d6e-ebb4-11ea-1250-d598737b4049
md"Julia, then danced against the same song and recorded herself near [the sea](https://freesound.org/s/9332/) ([cc-by-nc](https://creativecommons.org/licenses/by-nc/3.0/)):"

# ╔═╡ fdb90c24-ebb6-11ea-1a47-63001e2b8819
x2 = load("../data/jazz-waves.ogg")

# ╔═╡ 997f7a16-ed4a-11ea-3a29-891b059c3667
x2_duration_seconds = length(x2) / x2.samplerate

# ╔═╡ a3c0c930-ed4a-11ea-15e7-cb31b4e9675c
typeof(x2), x2.samplerate, length(x2)

# ╔═╡ 57466fbe-ed5a-11ea-2da0-199907c6bece
md"Let's look at Julia's audio wave:"

# ╔═╡ fba9972c-ebb4-11ea-2a8d-31350dd2855b
plot(domain(x2), x2[:, 1], xlabel="time (s)", ylabel="amplitude", title="Jazz song at the Sea", linealpha=0.85, legend=false, fmt = :png, dpi=300)

# ╔═╡ fa3508de-ed50-11ea-3357-19fd3a4b0441
md"We would like to know how much lag there is between the 2 different recordings of the same song: $$x$$ and $$x2$$."

# ╔═╡ b3997960-ebb0-11ea-2f91-af92117a47a2
md"## Using Cross Correlation to align audio signals"

# ╔═╡ fa4385a0-ebdc-11ea-2683-4107d5129e7e
md"Visually, we could superimpose both signals and slide each other over the time axis until they match as closely as possible.

For example, take Jazmine's recording ($$x$$) and place it next to Julia's ($$x2$$). They are clearly not aligned."

# ╔═╡ 9701b1a8-ed51-11ea-3992-fb42d17e5d0f
begin
	x_right_padding_1 = zeros(Float32, length(x2) - length(x), 1)
	x_padded_1 = cat(x, x_right_padding_1, dims=1)
	mixed_1 = cat(x2, x_padded_1, dims=2)
	plot(domain(mixed_1), mixed_1, xlabel="time (s)", ylabel="amplitude", title="At the sea + Jazz song", linealpha=0.85, legend=false, fmt = :png, dpi=300)
end

# ╔═╡ 07457180-ed51-11ea-148b-5b8c40e3470f
md"Then slide Julia's recording ($$x2$$) $$0.1$$ seconds to the right, then again $$0.1$$ seconds, …

 $$1$$ second later:"

# ╔═╡ 2c083f2a-ed51-11ea-1884-33902424a474
begin
	x_left_padding_2 = zeros(Float32, convert(Int, 1 * x.samplerate), 1)
	x_right_padding_2 = zeros(Float32, length(x2) - length(x) - length(x_left_padding_2), 1)
	x_padded_2 = cat(x_left_padding_2, x, x_right_padding_2, dims=1)
	mixed_2 = cat(x2, x_padded_2, dims=2)
	plot(domain(mixed_2), mixed_2, xlabel="time (s)", ylabel="amplitude", title="At the sea + Jazz song. Offset 1 second", linealpha=0.85, legend=false, fmt = :png, dpi=300)
end

# ╔═╡ 05c3c5a0-ed56-11ea-24fd-2d6628bf80ba
md" $$2$$ seconds later:"

# ╔═╡ 2133853e-ed56-11ea-06a7-9bf1a022192c
begin
	x_left_padding_3 = zeros(Float32, convert(Int, 2 * x.samplerate), 1)
	x_right_padding_3 = zeros(Float32, length(x2) - length(x) - length(x_left_padding_3), 1)
	x_padded_3 = cat(x_left_padding_3, x, x_right_padding_3, dims=1)
	mixed_3 = cat(x2, x_padded_3, dims=2)
	plot(domain(mixed_2), mixed_3, xlabel="time (s)", ylabel="amplitude", title="At the sea + Jazz song. Offset 2 seconds", linealpha=0.85, legend=false, fmt = :png, dpi=300)
end

# ╔═╡ 5871916a-ed56-11ea-2db0-132308b8e7aa
md" $$2.5$$ seconds later both $$x$$ and $$x2$$ are perfectly aligned."

# ╔═╡ 594b2b64-ed56-11ea-104b-cbf7510a12b3
begin
	x_left_padding_4 = zeros(Float32, convert(Int, 2.5 * x.samplerate), 1)
	x_right_padding_4 = zeros(Float32, length(x2) - length(x) - length(x_left_padding_4), 1)
	x_padded_4 = cat(x_left_padding_4, x, x_right_padding_4, dims=1)
	mixed_4 = cat(x2, x_padded_4, dims=2)
	plot(domain(mixed_2), mixed_4, xlabel="time (s)", ylabel="amplitude", title="Ocean Wave + Jazz. Offset 2.5 seconds", linealpha=0.85, legend=false, fmt = :png, dpi=300)
end

# ╔═╡ 07e0ed2c-ed51-11ea-0dd4-cfca2a31506b
md"This sliding process is essentially what a cross-correlation measures. It tells us how similar 2 signals are when displaced relative to each other."

# ╔═╡ 1997e6ea-ebb1-11ea-0875-8ff96ec060dc
md"Mathematically, the [cross-correlation](https://en.wikipedia.org/wiki/Cross-correlation) of two temporal series $$u$$ and $$v$$ is defined as:"

# ╔═╡ 16b55988-ebb1-11ea-3b2b-6f432c3eb7d1
md"$$(u \star v)(t) \triangleq\ \int_{-\infty}^\infty \overline{u(\tau)} v(t + \tau) \, d\tau$$"

# ╔═╡ d4b1188e-ebb9-11ea-3487-bf00755dbd76
md"Where $$\overline{u(\tau)}$$ is the [complex conjugate](https://en.wikipedia.org/wiki/Complex_conjugate) of $$u(\tau)$$. And $$t$$ is the lag."

# ╔═╡ d8b9598a-ebb9-11ea-3f08-7d5e64c688b5
md"The cross-correlation ($$\star$$) of two functions is equivalent to the convolution ($$*$$) of one of the signals reversed in the time axis:

$$(u \star v)(t) =  u(−t) ∗ g (t)$$

Moreover, the Fourier transform ($$\mathcal{F}$$) of a convolution [is equivalent](https://en.wikipedia.org/wiki/Convolution_theorem) to multiplying the Fourier tranforms of the functions:

$$\mathcal{F}\{u * v\} = \mathcal{F}\{u\} \cdot \mathcal{F}\{v\}$$

$$(u \star v)(t) = \mathcal{F}^{-1}\big\{\mathcal{F}\{u(-t)\}\cdot\mathcal{F}\{v(t)\}\big\}$$

As such we can [implement the cross-correlation](https://dsp.stackexchange.com/questions/736/how-do-i-implement-cross-correlation-to-prove-two-audio-files-are-similar) with the following steps:"

# ╔═╡ 2786714e-ebba-11ea-2f24-9f3640dd4b29
md"1. zero padding the signals and reversing one of them:"

# ╔═╡ 608c49d2-ebba-11ea-00ba-c541b8589037
begin
	x_max_size = abs(length(x2) - length(x))

	x_pad_size = length(x) + 2 * length(x_right_padding_1)
	x_pad = SampledSignals.SampleBuf(zeros(Float32, x_pad_size, 1), x.samplerate)
	x_padded = vcat(x, x_pad)
	x_padded_reversed = reverse(x_padded, dims=1)

	x2_pad_size = length(x2)
	x2_pad = SampledSignals.SampleBuf(zeros(Float32, x2_pad_size, 1), x2.samplerate)
	x2_padded = vcat(x2, x2_pad)

	# Listen to this audio at 12 seconds. It is the same Jazz song (x) but reversed and with empty audio at the start
	x_padded_reversed
end

# ╔═╡ 6d193a4c-ebbb-11ea-2ba6-c1341160a063
md"2. Performing a convolution between both signals with their Fourier transforms: "

# ╔═╡ 5cbb311c-ebc2-11ea-3b5e-910d2b940555
convolution = irfft(rfft(x2_padded) .* rfft(x_padded_reversed), length(x_padded))

# ╔═╡ 7deba858-ebbb-11ea-2602-cffad18a927a
plot(domain(x2_padded), convolution[:, 1], xlabel="time (s)", ylabel="convolution(t)", legend=false, fmt = :png, dpi=300)

# ╔═╡ 25a1d500-ebbd-11ea-1c37-35974b7f928d
md"The point where the signals are most correlated corresponds to the peak of the convolution. This is the lag between $$x$$ and $$x2$$:"

# ╔═╡ 979fa4d2-ebbf-11ea-1139-b167cf216ac7
lag_index, peak_value = argmax(convolution[:, 1]), maximum(convolution[:, 1])

# ╔═╡ 7213820c-ebc3-11ea-3c87-d1dcb117d446
md"The lag is exactly 2.5 seconds. We can now use this value to align the audio/video files for Jazmine and Julia."

# ╔═╡ 8863e55c-ebc2-11ea-2495-a962207eacc0
lag_seconds = x_duration_seconds * lag_index / length(x)

# ╔═╡ Cell order:
# ╟─cf3ddc9c-ebab-11ea-1417-b7eb9e517b50
# ╟─2c57c0fe-ed5b-11ea-2837-81321cf21495
# ╟─34a5c400-ebd9-11ea-04cf-3118b3902453
# ╟─ae29c670-ebd9-11ea-3995-77063c4c3f4d
# ╟─eb4bd808-ebc9-11ea-3c71-e9850bb24137
# ╟─f7a9d13a-ebac-11ea-2abd-79345ce6fbef
# ╟─cfdabcf8-ebae-11ea-307a-7728476c9825
# ╟─b16a96ca-eba3-11ea-3a74-4fecb5dda838
# ╠═49b7f3f4-eba7-11ea-14f1-31d09ab786bb
# ╠═0342bfec-ebc2-11ea-15a2-c769c8039d17
# ╠═c10ff5da-ebae-11ea-3205-6d1f6666ca5c
# ╟─6e4d19da-ebaf-11ea-113d-55eca1da0769
# ╠═c7eccb4a-ebad-11ea-16f5-21299446a8a2
# ╟─becebb5c-ebaf-11ea-3e5d-4f91c13b24a1
# ╠═3ab06516-eba5-11ea-2316-732fe86075f6
# ╟─c0e56d6e-ebb4-11ea-1250-d598737b4049
# ╠═fdb90c24-ebb6-11ea-1a47-63001e2b8819
# ╠═997f7a16-ed4a-11ea-3a29-891b059c3667
# ╠═a3c0c930-ed4a-11ea-15e7-cb31b4e9675c
# ╟─57466fbe-ed5a-11ea-2da0-199907c6bece
# ╠═fba9972c-ebb4-11ea-2a8d-31350dd2855b
# ╟─fa3508de-ed50-11ea-3357-19fd3a4b0441
# ╟─b3997960-ebb0-11ea-2f91-af92117a47a2
# ╟─fa4385a0-ebdc-11ea-2683-4107d5129e7e
# ╟─9701b1a8-ed51-11ea-3992-fb42d17e5d0f
# ╟─07457180-ed51-11ea-148b-5b8c40e3470f
# ╟─2c083f2a-ed51-11ea-1884-33902424a474
# ╟─05c3c5a0-ed56-11ea-24fd-2d6628bf80ba
# ╟─2133853e-ed56-11ea-06a7-9bf1a022192c
# ╟─5871916a-ed56-11ea-2db0-132308b8e7aa
# ╟─594b2b64-ed56-11ea-104b-cbf7510a12b3
# ╟─07e0ed2c-ed51-11ea-0dd4-cfca2a31506b
# ╟─1997e6ea-ebb1-11ea-0875-8ff96ec060dc
# ╟─16b55988-ebb1-11ea-3b2b-6f432c3eb7d1
# ╟─d4b1188e-ebb9-11ea-3487-bf00755dbd76
# ╟─d8b9598a-ebb9-11ea-3f08-7d5e64c688b5
# ╟─2786714e-ebba-11ea-2f24-9f3640dd4b29
# ╠═608c49d2-ebba-11ea-00ba-c541b8589037
# ╟─6d193a4c-ebbb-11ea-2ba6-c1341160a063
# ╠═5cbb311c-ebc2-11ea-3b5e-910d2b940555
# ╠═7deba858-ebbb-11ea-2602-cffad18a927a
# ╟─25a1d500-ebbd-11ea-1c37-35974b7f928d
# ╠═979fa4d2-ebbf-11ea-1139-b167cf216ac7
# ╟─7213820c-ebc3-11ea-3c87-d1dcb117d446
# ╠═8863e55c-ebc2-11ea-2495-a962207eacc0


text:
function complete(D::MaskedEDM{T}, loss::TripletEmbeddings.AbstractLoss, d::Int; kwargs...) where T <: Real
    @assert d ≥ 1 "d must be ≥ 1"

    triplets = Triplets(D)
    X = Embedding(d, size(D, 1))
    _ = fit!(loss, triplets, X; kwargs...)

    D̂ = EuclideanDistanceMatrix(X)
    s = mean(filter(!isnan, D./D̂))

    return s * D̂
end


text:
@doc raw"""
    solveNFE(problem,saveat)

# Arguments:
- `problem :: ProbOutput1D`: Output of probNFE using Input1D structure
- `problem :: ProbOutput2D`: Output of probNFE using Input2D structure
- `saveat  :: AbstractVector`: Vector containing the instants where the solution is saved at

Returns a structure containing the solution to the NFE saved at saveat instants.
The structure has `x`, `y` and `saveat` as fields to help plotting the solution.

# Examples
```julia-repl
julia> # 2D example
julia> saveat = [5.0,20.0]
julia> sol    = solveNFE(prob,saveat) # Deterministic solution saved at t=5 and t=20
julia> sol(20.0) # Solution at t=20
julia> x = sol.x # Spatial vector in x direction
julia> y = sol.y # Spatial vector in y direction
julia> using Plots
julia> plot(x,y,sol(20.0),st=:surface,title="Solution at t=20") # Surface plot of sol(20.0)
```
If the neural field is in 1D/2D, the solution is a vector/matrix.

-------------------------------------------------------------------------------------------
    solveNFE(problem,saveat,ϵ,np,ξ=0.1)

Solve stochastic version of an NFE for np trajectories, noise level ϵ and correlation ξ.

# Arguments:
- `problem :: ProbOutput1D`: Output of probNFE with Input1D
- `problem :: ProbOutput2D`: Output of probNFE with Input2D
- `saveat  :: AbstractVector`: Vector containing the instants where the solution is saved at
- `ϵ       :: Number`: Level of additive noise
- `np      :: Integer`: Number of simulations
- `ξ       :: AbstractFloat=0.1`: Spatial correlation parameter

Returns a structure containing the mean solution and the trajectories of the
NFE saved at saveat instants.

# Examples
```julia-repl
julia> # Stochastic solution saved at t=5,20, with noise level 0.01 simulated 100 times.
julia> sol_sto = solveNFE(prob,saveat,0.01,100)
julia> sol_sto(20.0,4) # 4th path at t=20
julia> sol_sto(5.0)    # Mean solution at t=20
```
"""
function solveNFE(prob::ProbOutput1D,saveat) # Deterministic method for 1D
    P      = prob.Plan
    Pinv   = prob.PlanInv
    krings = prob.krings
    α      = prob.α
    rings  = prob.Ω.rings1D
    N      = prob.Ω.N
    t      = prob.Ω.t
    x      = prob.Ω.x
    dt     = prob.Ω.dt
    V      = prob.V0 # Initial condition in natural space (needed if saveat contains t[1])
    v      = prob.v0 # Initial condition in Fourier space
    sv     = prob.sv # Firing Rate in Fourier space (all delays)
    hN     = N÷2+1   # Half of dim N (due to the output of rfft)

    # Pre-allocate vectors and matrices
    Vj  = Vector{Float64}(undef,N*length(saveat)) # Store the solution at saveat
    svu = Vector{ComplexF64}(undef,hN) # Firing Rate in Fourier space (delay u)
    I   = similar(V)
    î   = similar(svu) # rfft of I
    a   = similar(svu) # rfft of the integral operator A(X,t)
    j   = 1 # Index for saveat

    @inbounds for i = 1:length(t) # Time loop
        ti = t[i]
        # Store solution V at saveat instants
        if ti in saveat
            Vj[(j-1)*N+1:N*j] .= V
            j += 1
        end

        I.= [prob.I(i,ti) for i in x] # Discretise I
        mul!(î,P,I)                   # Compute rfft of I

        # Compute the integral operator's at frequency domain at t_i
        @. a = @views krings[1:hN]*sv[1:hN] # Init a
        @inbounds for u = 2:rings # Rings loop
            @. a += @views krings[1+hN*(u-1):hN*u]*sv[1+hN*(u-1):hN*u]
        end

        # Euler explicit scheme: v_i+1 = v_i + (Δt/α)*(î_i - v_i + a_i)
        @. v = v + (dt/α)*(î - v + a)

        V .= Pinv * v # Perform inverse Fourier transform to compute V in natural domain

        # Update sv
        # Every delay block moves one block down. The last one is deleted
        sv[hN+1:end] .= @view sv[1:hN*(rings-1)]
        mul!(svu,P,prob.S.(V)) # Apply the Real Fourier Transform to V
        sv[1:hN] .= svu        # Update sv first delay ring
    end #end time loop

    # Build our solution structure output
    sol = SolveOutDet1D(Vj,x,t,saveat,N)

    return sol

end

# Stochastic method for 1D
function solveNFE(prob::ProbOutput1D,saveat,ϵ,np,ξ=0.1)
    P      = prob.Plan
    Pinv   = prob.PlanInv
    krings = prob.krings
    α      = prob.α
    rings  = prob.Ω.rings1D
    N      = prob.Ω.N
    t      = prob.Ω.t
    x      = prob.Ω.x
    dt     = prob.Ω.dt
    hN     = N÷2+1   # Half of dim N (due to the output of rfft)

    # Pre-allocate matrices
    meanVj = Vector{Float64}(undef,N*length(saveat)) # Store the mean solution at saveat
    Vj  = Matrix{Float64}(undef,N*length(saveat),np) # Store the solution at saveat at path p
    V   = Vector{Float64}(undef,N)
    sv  = Vector{ComplexF64}(undef,hN*rings) # Firing Rate in Fourier space (all delays)
    svu = Vector{ComplexF64}(undef,hN) # Firing Rate in Fourier space (delay u)
    I   = similar(V)
    v   = similar(svu)
    î   = similar(svu) # rfft of I
    a   = similar(svu) # Fourier coefficients of the integral operator A(X,t)
    λ   = [exp(-(i^2*ξ^2)/(8*pi)) for i = 1:hN] # Correlation matrix

    # Pre-compute constant term of the stochastic part
    # fft gives the fourier coeff mult by N. Need to scale up noise, to have same magnitude
    # √dt due to the explicit Euler-Maruyama
    # Shift λ due to the also shifted output of fft
    c = N*sqrt(dt)*ϵ*fftshift(λ) # ϵ is the level of additive noise

    # Trajectories loop
    @inbounds for p = 1:np
        copy!(V,prob.V0)  # Initial condition V(X,0) = V0 (Natural domain)
        copy!(v,prob.v0)  # Initial condition v(x,0) = v0 (Fourier domain)
        copy!(sv,prob.sv) # Initialise with the initial values of sv

        j = 1 # Index for tj
        @inbounds for i = 1:length(t) # Time loop
            # w ~ N(0,1) + i*N(0,1)  stochastic term
            w = randn(Float64,hN) + im*randn(Float64,hN)
            ti = t[i]
            # Store solution V in saveat instants
            if ti in saveat
                Vj[(j-1)*N+1:N*j,p] .= V
                j += 1
            end

            I.= [prob.I(i,ti) for i in x] # Discretise I
            mul!(î,P,I)                   # rfft I

            # Compute the integral operator's at frequency domain at t_i
            @. a = @views krings[1:hN]*sv[1:hN] # Init a
            @inbounds for u = 2:rings # Ring loop
                @. a += @views krings[1+hN*(u-1):hN*u]*sv[1+hN*(u-1):hN*u]
            end

            # Euler-Maruyama explicit scheme: v_i+1 = v_i + (Δt/α)*(î_i - v_i + a_i) + √(Δt)*ϵ*λ*w
            @. v = v + (dt/α)*(î - v + a) + c*w

            V .= Pinv * v # Perform inverse Fourier transform to compute V in natural domain

            # Update sv
            # Every delay block moves one block down. The last one is deleted
            sv[hN+1:end] .= @view sv[1:hN*(rings-1)]
            mul!(svu,P,prob.S.(V)) # Apply the Real Fourier Transform to V
            sv[1:hN,:].= svu       # Update sv first delay ring
        end #end time loop

    end #end trajectories loop

    # Compute the mean solution across all trajectories
    @inbounds for i = 1:N*length(saveat)
        meanVj[i] = mean(@view Vj[i,:])
    end

    # Build our solution structure output
    sol = SolveOutSto1D(Vj,meanVj,x,t,saveat,N)

    return sol

end


text:
using Molly
using Aqua
import BioStructures # Imported to avoid clashing names
using CUDA
using FiniteDifferences
using ForwardDiff
using Zygote

using Base.Threads
using DelimitedFiles
using LinearAlgebra
using Statistics
using Test

@warn "This file does not include all the tests for Molly.jl due to CI time limits, " *
        "see the test directory for more"

# Allow testing of particular components
const GROUP = get(ENV, "GROUP", "All")
if GROUP == "Zygote"
    @warn "Only running Zygote tests as GROUP is set to Zygote"
end

# GLMakie doesn't work on CI
run_visualize_tests = !haskey(ENV, "CI")
if run_visualize_tests
    using GLMakie
    @info "The visualization tests will be run as this is not CI"
else
    @warn "The visualization tests will not be run as this is CI"
end

run_parallel_tests = nthreads() > 1
if run_parallel_tests
    @info "The parallel tests will be run as Julia is running on $(nthreads()) threads"
else
    @warn "The parallel tests will not be run as Julia is running on 1 thread"
end

run_gpu_tests = CUDA.functional()
if run_gpu_tests
    @info "The GPU tests will be run as a CUDA-enabled device is available"
else
    @warn "The GPU tests will not be run as a CUDA-enabled device is not available"
end

CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU

data_dir = normpath(@__DIR__, "..", "data")

temp_fp_pdb = tempname(cleanup=true) * ".pdb"
temp_fp_viz = tempname(cleanup=true) * ".mp4"

if GROUP == "All"
    # Some failures due to dependencies but there is an unbound args error
    Aqua.test_all(
        Molly;
        ambiguities=(recursive=false),
        unbound_args=false,
        undefined_exports=false,
    )

    include("basic.jl")
    include("simulation.jl")
    include("protein.jl")
    include("agent.jl")
end

if GROUP == "All" || GROUP == "Zygote"
    include("zygote.jl")
end


text:
module AMLPipelineBase

include("abstracttypes.jl")
using .AbsTypes
export Machine, Computer, Workflow, Learner, Transformer
export fit, fit!, transform, transform!, fit_transform, fit_transform!

include("utils.jl")
using .Utils
export holdout, kfold, score, infer_eltype,
       nested_dict_to_tuples,
       nested_dict_set!,
       nested_dict_merge,
       create_transformer,
       mergedict, getiris,
       skipmean,skipmedian,skipstd,
       aggregatorclskipmissing,
       getiris, getprofb,
       find_catnum_columns,
       train_test_split

include("baselinemodels.jl")
using .BaselineModels
export Baseline, Identity

include("basefilters.jl")
using .BaseFilters
export Imputer, OneHotEncoder, Wrapper

include("featureselector.jl")
using .FeatureSelectors
export FeatureSelector, CatFeatureSelector, NumFeatureSelector, CatNumDiscriminator

include("decisiontree.jl")
using .DecisionTreeLearners
export PrunedTree, RandomForest, Adaboost

include("ensemble.jl")
using .EnsembleMethods
export VoteEnsemble, StackEnsemble, BestLearner

include("crossvalidator.jl")
using .CrossValidators
export crossvalidate, pipe_performance

include("naremover.jl")
using .NARemovers
export NARemover

include("pipelines.jl")
using .Pipelines
export @pipeline, @pipelinex, @pipelinez
export |>, +, |, *, >>
export Pipeline, ComboPipeline

end # module


text:
amap = Dict(
rs_st1_1 => [rs_state1=>rs_state1],
rs_st1_2 => [rs_state1=>rs_state2],
rs_st2_1 => [rs_state2=>rs_state1],
rs_st2_2 => [rs_state2=>rs_state2],
rs_st13_2 => [rs_state1=>rs_state2,rs_state3=>rs_state2],
rs_st1_3 => [rs_state1=>rs_state3],
rs_st2_3 => [rs_state2=>rs_state3],
rs_st3_3 => [rs_state3=>rs_state3],
rs_st123_3 => [rs_state1=>rs_state3,rs_state2=>rs_state3,rs_state3=>rs_state3],
)
fsmthreestate = FiniteStateMachine(amap,rs_state2,rs_state1)

@test !finished(fsmthreestate)

@test can(fsmthreestate,rs_st2_2)
@test can(fsmthreestate,rs_st2_3)
@test can(fsmthreestate,rs_st123_3)
@test !can(fsmthreestate,rs_st13_2)
@test !can(fsmthreestate,rs_st1_1)
@test !can(fsmthreestate,rs_st3_3)

@test fire!(fsmthreestate,rs_st2_2)
@test fsmthreestate.current == rs_state2
@test fire!(fsmthreestate,rs_st123_3)
@test fsmthreestate.current == rs_state3
@test fire!(fsmthreestate,rs_st13_2)
@test fsmthreestate.current == rs_state2
@test !finished(fsmthreestate)
@test fire!(fsmthreestate,rs_st2_1)
@test fsmthreestate.current == rs_state1
@test finished(fsmthreestate)


text:
# "Dirichlet distribution"
# abstract type Dirichlet <: Dist end

# function dirichlet(ω::Ω, α)
#   gammas = [gammarv(ω[@uid][i], αi, 1.0) for (i, αi) in enumerate(α)]
#   Σ = sum(gammas)
#   [gamma/Σ for gamma in gammas]
# end
# # FIXME: Type
# dirichlet(α::MaybeRV{T}) where T = RandVar{T, Dirichlet}(dirichlet, (α,))

# _rand!(rng::AbstractRNG, d::Djl.MvNormal, x::VecOrMat) = add!(unwhiten!(d.Σ, randn!(rng, x)), d.μ)

# mvchol(x, μ, Σ::PDMat) = Distributions.unwhiten(Σ, x) .+ μ

# "Multivariate Normal Distribution with mean vector `μ` and covariance `Σ`"
# abstract type MvNormal <: Dist end

# mvnormal(ω::Ω, μ::Vector, Σ::PDMat) = mvchol(normal(ω, 0.0, 1.0, size(μ)), μ, Σ)
# # mvnormal(ω::Ω, μ, Σ) = rand(ω, MvNormal(μ, Σ))
# mvnormal(μ::MaybeRV{T1}, Σ::MaybeRV{T2}) where {T1, T2} =
#   RandVar{T1, MvNormal}(mvnormal, (μ, PDMat(Σ)))

# lift(:PDMat, 1)


text:
#  File: combined.jl
#  Author: www.erinzuun.com
#  Date: 21-Jun-2020

function barchart()
    #  import XL data
    [~,dates,returns] = getXLData('hfma_matlab_data.xlsx','CTA Index');

    n = size(returns,1); #  #  data points

    figure; #  create figure
    sDate = datenum(dates(1)); #  set start date for x-axis
    eDate = datenum(dates(end)); #  set end date for x-axis
    xData = linspace(sDate,eDate,n);

    bar(xData, returns); #  plot bar chart
    xlabel('Date'); #  add x label
    ylabel('RoR (#  )'); #  add y label
    box off;
    ytick = get(gca,'YTick'); #  format axes
    set(gca,'YTickLabel', sprintf('#  .2f|',ytick))
    datetick('x','yyyy','keeplimits');
    return bar
end

using StatsBase, StatsPlots, Distributions

function histogram()
    # import XL file
    [dates,returns] = getXLData('hfma_matlab_data.xlsx','CTA Index');

    mn = abs(min(returns)); #  get min percent absolute
    mx = abs(max(returns)); #  get max percent absolute

    figure; #  create figure
    y=returns; #  set returns to y

    x=-4.5:0.25:4.5;  #  bins -4.5 to 4.5 intervals 0.25
    freq = histc(y,x);  #  count frequency per bin
    freq = y.count(x);
    rfreq = (freq/sum(freq))*100; #  relative frequency

    bar(x,rfreq);  #  plot bar chart
    xlabel('Bin'); #  add x label
    ylabel('Rel. Freq. (#  )'); #  add y label
end

function barchart()
    #  import XL data
    [~,dates,returns] = getXLData('hfma_matlab_data.xlsx','CTA Index');

    n = size(returns,1); #  #  data points

    figure; #  create figure
    sDate = datenum(dates(1)); #  set start date for x-axis
    eDate = datenum(dates(end)); #  set end date for x-axis
    xData = linspace(sDate,eDate,n);

    bar(xData, returns); #  plot bar chart
    xlabel('Date'); #  add x label
    ylabel('RoR (#  )'); #  add y label
    box off;
    ytick = get(gca,'YTick'); #  format axes
    set(gca,'YTickLabel', sprintf('#  .2f|',ytick))
    datetick('x','yyyy','keeplimits');
end

#  EXCESS KURTOSIS
function fKurt(x,flag,f)
#  --------------------------------------------------------------------------
#  x:     returns
#  flag:  0 = sample, 1 = population
#  f:     reporting frequency
#  m4:    excess kurtosis (sample or population)
#  --------------------------------------------------------------------------

    m4 = (kurtosis(x,flag)-3)./f
    return m4
end

#  MEAN
function fMean(x,f)
#  -------------------------------------------------------------------
#  x:     returns
#  f:     reporting frequency
#  m1:    mean
#  -------------------------------------------------------------------

    m1 = mean(x).*f
    return m1
end

# STANDARD DEVIATION
function fStd(x,flag,f)
# --------------------------------------------------------------------------
# x:     returns
# flag:  0 = sample, 1 = population
# f:     reporting frequency
# m2:    standard deviation (sample or population)
# --------------------------------------------------------------------------

    m2 = std(x,flag).*sqrt(f)
    return m2
end

#  VALUE ADDED INDEX (VAI)
function fVAI(x)

#  --------------------------------------------------------------------------
#  x:       returns vector
#  vaidx:   Value Added Index (VAI)
#  --------------------------------------------------------------------------

    n = size(x,1); #  data points
    fVAI = ones(1,n); #  pre-allocate memory()
    fVAI(1) = 100; #  initialise to $100
    for i in 1:n
        fVAI(i+1) = fVAI(i) * (1+x[i]./100); #  calculate VAI
    end
    vaidx = fVAI;
    return vaidx
end

#  GET XL DATA
function getXLData(fname,sname)
    #  raw data
    rawdata = XLSX.readxlsx(fname)
    #  labels
    names = rawdata(1,:);
    names(1) = [];
    #  dates
    dates = rawdata(2:end,1);
    #  returns
    rwdata(:,1) = [];
    returns = rawdata(2:end,:);
    # #
    # returns = cell2mat(returns);
end

#  import XL data
[~,dates,returns] = getXLData('hfma_matlab_data.xlsx','CTA Index');

figure; # create figure
y=returns; # set returns to y

x=-4.5:0.25:4.5; # bins -4.5 to 4.5 intervals 0.25
freq = histc(y,x); # count frequency per bin
rfreq = (freq/sum(freq))*100; # relative frequency
nfreq=normpdf(x); # standard normal frequency per bin
rnfreq = (nfreq/sum(nfreq))*100; # relative normal frequency

bar(x,rfreq); # plot chart combination
hold on;
plot(x,rnfreq,'r');
hold off;
xlabel('Bin'); # add x label
ylabel('Rel. Freq. (# )'); # add y label


text:
function ellipk(x::DoubleFloat{T}) where {T<:IEEEFloat}
    !(-0.0 <= x <= 1.0) && throw(DomainError("$x not in 0..1"))
    x === one(T) && return DoubleFloat{T}(Inf)

    df_one    = one(DoubleFloat{T})
    df_halfpi = DoubleFloat{T}(pi)/2

    k = sqrt(abs(x - df_one))
    k = agm1(k)
    k = df_halfpi / k

    return k
end


text:
# This file is part of Kpax3. License is MIT.

function test_mcmc_split_init()
  # split cluster 1
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=2, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [13; 13; 13; 13; 13; 76]
  k = length(unique(R)) + 1

  ij = [1; 5]
  S = 3
  u = 4

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportsplit!(ij, k, data, priorC, settings, support)

  len = 4

  g = 0
  lp = zeros(Float64, 4, len, m)
  for b in 1:m, l in 1:state.k
    g = state.cl[l]

    lp[1, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[1, b], priorC.B[1, b])
    lp[2, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[2, b], priorC.B[2, b])
    lp[3, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[3, b], priorC.B[3, b])
    lp[4, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[4, b], priorC.B[4, b])
  end
  @test support.lp == lp

  wi = zeros(Float64, 4, m)
  for col in 1:m
    wi[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col])
    wi[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col])
    wi[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col])
    wi[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col])
  end

  ci = Float64[log(sum(exp.(support.wi.w[:, b]))) for b in 1:m]

  wj = zeros(Float64, 4, m)
  for col in 1:m
    wj[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col])
    wj[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col])
    wj[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col])
    wj[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col])
  end

  cj = Float64[log(sum(exp.(support.wj.w[:, b]))) for b in 1:m]

  @test support.vi == 1
  @test support.ni == float(data[:, ij[1]])
  @test support.ui == [ij[1]; 0; 0; 0; 0; 0]

  @test maximum(abs.(support.wi.w - wi)) <= ε
  @test maximum(abs.(support.wi.c - ci)) <= ε
  @test support.wi.z == zeros(Float64, 4, m)

  @test support.vj == 1
  @test support.nj == float(data[:, ij[2]])
  @test support.uj == [ij[2]; 0; 0; 0; 0; 0]

  @test maximum(abs.(support.wj.w - wj)) <= ε
  @test maximum(abs.(support.wj.c - cj)) <= ε
  @test support.wj.z == zeros(Float64, 4, m)

  nothing
end

test_mcmc_split_init()

function test_mcmc_split_updatei()
  # move the first unit (u) to cluster 1
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=2, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [13; 13; 13; 13; 13; 76]
  k = length(unique(R)) + 1

  ij = [1; 5]
  S = 3
  u = 4

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportsplit!(ij, k, data, priorC, settings, support)

  len = 4

  g = 0
  lp = zeros(Float64, 4, len, m)
  for b in 1:m, l in 1:state.k
    g = state.cl[l]

    lp[1, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[1, b], priorC.B[1, b])
    lp[2, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[2, b], priorC.B[2, b])
    lp[3, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[3, b], priorC.B[3, b])
    lp[4, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[4, b], priorC.B[4, b])
  end
  @test support.lp == lp

  lcp = zeros(Float64, 2)
  for b in 1:m
    lcp[1] += Kpax3.computeclusteriseqprobs!(data[b, u], b, priorC, support)
    lcp[2] += Kpax3.computeclusterjseqprobs!(data[b, u], b, priorC, support)
  end
  Kpax3.updateclusteri!(u, data, support)

  wi = zeros(Float64, 4, m)
  for col in 1:m
    wi[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col])
    wi[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col])
    wi[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col])
    wi[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col])
  end
  zi = copy(wi)

  ci = Float64[log(sum(exp.(support.wi.w[:, b]))) for b in 1:m]

  wj = zeros(Float64, 4, m)
  zj = zeros(Float64, 4, m)
  for col in 1:m
    wj[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col])
    wj[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col])
    wj[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col])
    wj[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col])

    zj[1, col] = wj[1, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col])
    zj[2, col] = wj[2, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col])
    zj[3, col] = wj[3, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col])
    zj[4, col] = wj[4, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col])
  end

  cj = Float64[log(sum(exp.(support.wj.w[:, b]))) for b in 1:m]

  @test support.vi == 2
  @test support.ni == float(data[:, ij[1]]) + float(data[:, u])
  @test support.ui == [ij[1]; u; 0; 0; 0; 0]

  @test maximum(abs.(support.wi.w - wi)) <= ε
  @test maximum(abs.(support.wi.c - ci)) <= ε
  @test maximum(abs.(support.wi.z - zi)) <= ε

  @test support.vj == 1
  @test support.nj == float(data[:, ij[2]])
  @test support.uj == [ij[2]; 0; 0; 0; 0; 0]

  @test maximum(abs.(support.wj.w - wj)) <= ε
  @test maximum(abs.(support.wj.c - cj)) <= ε
  @test maximum(abs.(support.wj.z - zj)) <= ε

  nothing
end

function test_mcmc_split_updatej()
  # move the first unit (u) to cluster 2
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=2, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [13; 13; 13; 13; 13; 76]
  k = length(unique(R)) + 1

  ij = [1; 5]
  S = 3
  u = 4

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportsplit!(ij, k, data, priorC, settings, support)

  len = 4

  g = 0
  lp = zeros(Float64, 4, len, m)
  for b in 1:m, l in 1:state.k
    g = state.cl[l]

    lp[1, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[1, b], priorC.B[1, b])
    lp[2, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[2, b], priorC.B[2, b])
    lp[3, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[3, b], priorC.B[3, b])
    lp[4, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[4, b], priorC.B[4, b])
  end
  @test support.lp == lp

  lcp = zeros(Float64, 2)
  for b in 1:m
    lcp[1] += Kpax3.computeclusteriseqprobs!(data[b, u], b, priorC, support)
    lcp[2] += Kpax3.computeclusterjseqprobs!(data[b, u], b, priorC, support)
  end
  Kpax3.updateclusterj!(u, data, support)

  wi = zeros(Float64, 4, m)
  zi = zeros(Float64, 4, m)
  for col in 1:m
    wi[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col])
    wi[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col])
    wi[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col])
    wi[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col])

    zi[1, col] = wi[1, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col])
    zi[2, col] = wi[2, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col])
    zi[3, col] = wi[3, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col])
    zi[4, col] = wi[4, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col])
  end

  ci = Float64[log(sum(exp.(support.wi.w[:, b]))) for b in 1:m]

  wj = zeros(Float64, 4, m)
  for col in 1:m
    wj[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col])
    wj[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col])
    wj[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col])
    wj[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col])
  end
  zj = copy(wj)

  cj = Float64[log(sum(exp.(support.wj.w[:, b]))) for b in 1:m]

  @test support.vi == 1
  @test support.ni == float(data[:, ij[1]])
  @test support.ui == [ij[1]; 0; 0; 0; 0; 0]

  @test maximum(abs.(support.wi.w - wi)) <= ε
  @test maximum(abs.(support.wi.c - ci)) <= ε
  @test maximum(abs.(support.wi.z - zi)) <= ε

  @test support.vj == 2
  @test support.nj == float(data[:, ij[2]]) + float(data[:, u])
  @test support.uj == [ij[2]; u; 0; 0; 0; 0]

  @test maximum(abs.(support.wj.w - wj)) <= ε
  @test maximum(abs.(support.wj.c - cj)) <= ε
  @test maximum(abs.(support.wj.z - zj)) <= ε

  nothing
end

test_mcmc_split_updatej()

function test_mcmc_updatelogmarglik()
  # move the first unit (u) to cluster 2 (test inverse split operator)
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=2, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [13; 13; 13; 13; 13; 76]
  k = length(unique(R)) + 1

  ij = [1; 4]
  S = 3
  u = 4

  ni = vec(sum(float(data[:, [1; 2; 3]]), dims=2))
  vi = 3

  nj = vec(sum(float(data[:, [4; 5]]), dims=2))
  vj = 2

  lpi = zeros(Float64, 4, m)
  lpj = zeros(Float64, 4, m)
  for b in 1:m
    lpi[1, b] = Kpax3.logmarglik(ni[b], vi, priorC.A[1, b], priorC.B[1, b])
    lpi[2, b] = Kpax3.logmarglik(ni[b], vi, priorC.A[2, b], priorC.B[2, b])
    lpi[3, b] = Kpax3.logmarglik(ni[b], vi, priorC.A[3, b], priorC.B[3, b])
    lpi[4, b] = Kpax3.logmarglik(ni[b], vi, priorC.A[4, b], priorC.B[4, b])

    lpj[1, b] = Kpax3.logmarglik(nj[b], vj, priorC.A[1, b], priorC.B[1, b])
    lpj[2, b] = Kpax3.logmarglik(nj[b], vj, priorC.A[2, b], priorC.B[2, b])
    lpj[3, b] = Kpax3.logmarglik(nj[b], vj, priorC.A[3, b], priorC.B[3, b])
    lpj[4, b] = Kpax3.logmarglik(nj[b], vj, priorC.A[4, b], priorC.B[4, b])
  end

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportsplit!(ij, k, data, priorC, settings, support)

  support.ni = copy(ni)
  support.vi = copy(vi)

  support.nj = copy(nj)
  support.vj = copy(vj)

  Kpax3.updatelogmargliki!(priorC, support)
  Kpax3.updatelogmarglikj!(priorC, support)

  @test support.lpi == lpi
  @test support.lpj == lpj

  fill!(support.lpi, 0.0)
  fill!(support.lpj, 0.0)

  Kpax3.updatelogmargliki!(ni, vi, priorC, support)
  Kpax3.updatelogmarglikj!(nj, vj, priorC, support)

  @test support.lpi == lpi
  @test support.lpj == lpj

  nothing
end

test_mcmc_updatelogmarglik()


text:
# use this to make small messages fast
function message(x)
  m = Message(sizeof(x))
  mb = Base.AbstractIOBuffer(m,true,true,true,false,length(m))
  write(mb,x)
  m
end

struct ZMQDataSink
  s::ZMQ.Socket
  channel_number::Int
end
function Base.write(zds::ZMQDataSink, dp)
  send_multipart(zds.s, [message("$(zds.channel_number)"), message(dp)])
end
function send_multipart(socket::Socket, parts::Vector{Message})
  for msg in parts[1:end-1]
   send(socket, msg, SNDMORE)
  end
  return send(socket, parts[end])
end
function write_header(zds::ZMQDataSink, r)
  send_multipart(zds.s,message.(["header$(zds.channel_number)","write_header called"]))
end
function write_header_end(zds::ZMQDataSink, ljh, analyzer)
  send_multipart(zds.s,message.(["header$(zds.channel_number)","write_header_end called"]))
end
function write_header_allchannel(zds::ZMQDataSink, x...)
  send_multipart(zds.s, message.(["header_all_channel", "write_header_allchannel called"]))
end
function Base.close(zds::ZMQDataSink)
  # for now don't actually close the socket in case, since there is no mechanism
  # to reopen it
  # closing the same ZMQ socket multiple times does not error or hang
end
Base.flush(zds::ZMQDataSink) = nothing

module ZMQDataSinkConsts
using ZMQ
initialized = false
context = nothing
socket = nothing
port = nothing
  "init_for_zmqdatasink(portin::Int;verbose=false, error_on_already_initialized=false)"
  function init_for_zmqdatasink(portin::Int;verbose=false, error_on_already_initialized=false)
    if initialized && error_on_already_initialized
      error("ZMQDataSinkConsts already initialized")
    elseif initialized
      return nothing
    end
    global port = portin
    global context = ZMQ.Context()
    global socket = ZMQ.Socket(context, ZMQ.PUB)
    ZMQ.bind(socket, "tcp://*:$port")
    global initialized = true
    verbose && println("Pope ZMQ publisher did bind on port $port")
    return nothing
  end
end
const init_for_zmqdatasink = ZMQDataSinkConsts.init_for_zmqdatasink

"make_zmqdatasink(channel::Int)
must call `init_for_zmqdatasink(port)` first
Returns a ZMQDataSink for `channel`"
function make_zmqdatasink(channel::Int)
  ZMQDataSinkConsts.initialized || error("must call init_for_zmqdatasink(port) first")
  ZMQDataSink(ZMQDataSinkConsts.socket,channel)
end

"make_buffered_hdf5_and_zmq_multi_sink(output_file, channel_number)
makes a BufferedHDF5Writer and a ZMQDataSink, for the same channel,
puts them both into a MultipleDataSink."
function make_buffered_hdf5_and_zmq_multisink(output_file, channel_number, analyzer)
  product_writer_a = Pope.make_buffered_hdf5_writer(output_file, channel_number, analyzer)
  product_writer_b = Pope.make_zmqdatasink(channel_number)
  Pope.MultipleDataSink(product_writer_a,product_writer_b)
end


text:
"""
Constant (forced) assimilation, given in ``μmol\\ m^{-2}\\ s^{-1}``.

# Examples

```julia
ConstantA(30.0)
```
"""
Base.@kwdef struct ConstantA{T} <: AbstractAModel
    A::T = 25.0
end

function inputs(::ConstantA)
    (:Cₛ,)
end

function outputs(::ConstantA)
    (:A, :Gₛ, :Cᵢ)
end

Base.eltype(x::ConstantA) = typeof(x).parameters[1]

"""
    photosynthesis!_(leaf::LeafModels{I,E,<:ConstantA,<:AbstractGsModel,S},constants = Constants())

Constant photosynthesis.

# Returns

Modify the leaf status in place for A, Gₛ and Cᵢ:

- A: carbon assimilation, set to leaf.photosynthesis.A (μmol[CO₂] m-2 s-1)
- Gₛ: stomatal conductance for CO₂ (mol[CO₂] m-2 s-1)
- Cᵢ: intercellular CO₂ concentration (ppm)

# Arguments

- `leaf::LeafModels{.,.,<:ConstantA,<:AbstractGsModel,.}`: A [`LeafModels`](@ref) struct holding the parameters for
the model with initialisations for:
    - `Cₛ` (mol m-2 s-1): surface CO₂ concentration.
    - `Dₗ` (mol m-2 s-1): vapour pressure difference between the surface and the air saturation
    vapour pressure in case you're using the stomatal conductance model of [`Medlyn`](@ref).
- `meteo`: meteorology structure, see [`Atmosphere`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

# Note

`Cₛ` (and `Dₗ` if you use [`Medlyn`](@ref)) must be initialised by providing them as keyword
arguments (see examples). If in doubt, it is simpler to compute the energy balance of the
leaf with the photosynthesis to get those variables. See [`energy_balance`](@ref) for more
details.

# Examples

```julia
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)
leaf = LeafModels(photosynthesis = ConstantA(),
            stomatal_conductance = Medlyn(0.03, 12.0),
            Cₛ = 400.0)

photosynthesis!_(leaf,meteo,Constants())

leaf.status.A
```
"""
function photosynthesis!_(leaf::LeafModels{I,E,<:ConstantA,<:AbstractGsModel,S}, meteo,
    constants = Constants()) where {I,E,S}

    # Net assimilation (μmol m-2 s-1)
    leaf.status.A = leaf.photosynthesis.A

    # Stomatal conductance (mol[CO₂] m-2 s-1)
    leaf.status.Gₛ = gs(leaf, meteo)

    # Intercellular CO₂ concentration (Cᵢ, μmol mol)
    leaf.status.Cᵢ = min(leaf.status.Cₛ, leaf.status.Cₛ - leaf.status.A / leaf.status.Gₛ)

    nothing
end


text:
function get_output(coefs,Df,sq,r2,f_terms,Gres,Tot,p,n)
n_term = length(coefs)
DF = RegCol("Df")
R2 = RegCol("R²")
sumsq = RegCol("SumOfSqs")
F = RegCol("F")
P = RegCol("P")

for i in 1:n_term
setcoef!(R2, coefs[i]=>r2[i])
setcoef!(DF, coefs[i]=>Df[i])
setcoef!(sumsq, coefs[i]=>sq[i])
setcoef!(F, coefs[i]=>f_terms[i])
setcoef!(P, coefs[i]=>p[i])
end
coefArray = hcat(Df,sq,r2,f_terms,p)
coefArray= vcat(coefArray,[n-1-sum(Df) Gres 1-sum(r2) NaN NaN ])
coefArray= vcat(coefArray,[n-1 Tot 1 NaN NaN ])

setcoef!(R2, "Residual"=>1-sum(r2))
setcoef!(DF, "Residual"=>n-1-sum(Df))
setcoef!(sumsq, "Residual"=>Gres)
setcoef!(R2, "Total"=>1)
setcoef!(DF, "Total"=>n-1)
setcoef!(sumsq, "Total"=>Tot)
return NamedArray(coefArray,(vcat(coefs...,"Residual","Total"),["Df","SumOfSqs","R2","F","P"])), hcat(DF,sumsq,R2,F,P)
end


text:
using Statistics
using Plots
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using JLD2
using Random
using Distributions
using FITSIO
using SparseArrays

push!(LOAD_PATH, pwd()*"/../main")
using DHC_2DUtils
push!(LOAD_PATH, pwd()*"/../scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs

function kldiv_from_samples(samps, qGaussian)
    #size(samps) = (Nsamp, Ndim)
    #qGaussian is a function which acts on each row of samps
    #MC estimate of KLP||Q: cant do this naively
    Nsamps = size(samps)[1]
    logqeval = zeros(Nsamps)
    for s=1:Nsamps
        qgeval = qGaussian(samps[s, 1:end])
        println(size(qgeval), " Val", qgeval)
        logqeval[s] = qgeval
    end
    println(logqeval)
    logpeval = fill(-log.(Nsamps), Nsamps)
    kld = sum(logpeval .- logqeval)./Nsamps
    return kld
end

function Gaussian_expr(meantarg, invcov)
    function Gwrapper(x)
        return -0.5.* (x- meantarg)' * invcov * (x- meantarg)
    end
    return Gwrapper
end

Nx=64
im = readsfd(Nx)
true_img = im[:, :, 1]

noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
#heatmap(init)
#heatmap(true_img)

filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(S1iso, Nf) = size(filter_hash["S1_iso_mat"])
(S2iso, Nfsq) = size(filter_hash["S2_iso_mat"])
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true, :iso=>false) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
if dhc_args[:iso]
    coeff_mask = falses(2+S1iso+S2iso)
    coeff_mask[S1iso+3:end] .= true
else #Not iso
    coeff_mask = falses(2+Nf+Nf^2)
    coeff_mask[Nf+3:end] .= true
end

white_noise_args = Dict(:loc=>0.0, :sig=>std(true_img), :Nsam=>1000, :norm=>false, :smooth=>false, :smoothval=>0.8) #Iso #only if you're using noise based covar
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])
recon_settings = Dict([("target_type", "sfd_dbn"), ("covar_type", "sfd_dbn"), ("log", true), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", 0.00456), ("white_noise_args", white_noise_args)]) #Add constraints
regs_true = DHC_compute_wrapper(log.(true_img), filter_hash; norm=false, dhc_args...)[coeff_mask]
s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings, safety=nothing)
qg = Gaussian_expr(s2mean, s2icov)
psamps = get_dbn_coeffs(log.(im), filter_hash, dhc_args, coeff_mask=coeff_mask)
#println("Check ", qg())
println("KLDiv = ", kldiv_from_samples(psamps, qg))

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
using Permanents
using Test
using LinearAlgebra

const ATOL = 1e6*eps()

function fourier_matrix(n; normalized = true)

    U = Array{ComplexF64}(undef, n, n)

    for i in 1:n
        for j in 1:n
            U[i,j] = exp((2im * pi / n)* (i-1) * (j-1))
        end
    end

    if normalized
        return 1/sqrt(n) * U
    else
        return U
    end
end

function rand_haar(n::Int)

    """generates a Haar distributed unitary matrix of size n*n"""
    # follows https://case.edu/artsci/math/esmeckes/Meckes_SAMSI_Lecture2.pdf

    qr(randn(ComplexF64, n,n)).Q
end

function rand_gram_matrix(n::Int)

	"""random gram matrix with full rank """

	function normalized_random_vector(n)
		v = rand(ComplexF64, n)
		1/norm(v) .* v
	end

	generating_vectors = 0.
	while det(generating_vectors) ≈ 0. #check it's a basis
	    generating_vectors = hcat([normalized_random_vector(n) for i = 1:n]...)

	end

	generating_vectors' * generating_vectors

end

@testset "Permanents.jl" begin
    theoretical_permanents_fourier_matrix = [1,0,-3,0,-5,0,-105,0,81,0,6765,0,175747,0,30375,0,25219857,0,142901109,0,4548104883,0]
    #from "on the permanent of schur matrix", graham

    for n = 1:6

        @test naive(ones(Int,n,n)) == factorial(n)
        @test naive_tensor(ones(n,n,n)) ≈ factorial(n)^2

        U = fourier_matrix(n, normalized = false)
        @test naive(U) ≈ theoretical_permanents_fourier_matrix[n] atol=ATOL
        @test ryser(U) ≈ theoretical_permanents_fourier_matrix[n] atol=ATOL
        @test fast_glynn_perm(U) ≈ theoretical_permanents_fourier_matrix[n] atol=ATOL
        @test multi_dim_ryser(U, ones(n,n)) ≈ abs(theoretical_permanents_fourier_matrix[n]).^2 atol=ATOL
    end

	### consistency of ryser_tensor(W) ≈  multi_dim_ryser(U, S) ###

    U = fourier_matrix(7,normalized = false)
    @test glynn_precision(U; rtol = 1e-3, maxiter = 10^5) ≈ theoretical_permanents_fourier_matrix[7] rtol = 1e-2

	n = 7
	U = rand_haar(n)
	S = rand_gram_matrix(n)

	W = Array{eltype(S)}(undef, (n,n,n))

	for ss in 1:n
	    for rr in 1:n
	        for j in 1:n
	            W[ss,rr,j] = U[ss, j] * conj(U[rr, j]) * S[rr, ss]
	        end
	    end
	end

	@test ryser_tensor(W) ≈  multi_dim_ryser(U, S)

	U = rand(4,4)
	@test positive_entry(U) ≈ naive(U) atol=1e-2

end