text:
"""
    loglikelihood(v::mIHTVariable)

Calculates the loglikelihood of observing `Y` given mean `μ` and precision matrix
`Γ` (inverse covariance matrix) under a multivariate Gaussian.

Caution: This function mutates storage variables `v.r_by_r1` and `v.r_by_r2`
"""
function loglikelihood(v::mIHTVariable)
    mul!(v.r_by_r1, v.resid, Transpose(v.resid)) # r_by_r = (Y - BX)(Y - BX)'
    mul!(v.r_by_r2, v.Γ, v.r_by_r1) # r_by_r2 = Γ(Y - BX)(Y - BX)'
    return nsamples(v) / 2 * logdet(v.Γ) - 0.5 * tr(v.r_by_r2) # logdet allocates!
end

"""
    update_xb!(v::mIHTVariable)

Using only support of B, update `v.BX` with the new proposed `B` (effect size
for genetic covariates) and `C` (non-genetic covariates).
"""
function update_xb!(v::mIHTVariable)
    copyto!(v.Xk, @view(v.X[v.idx, :]))
    mul!(v.BX, view(v.B, :, v.idx), v.Xk)
    mul!(v.CZ, v.C, v.Z)
end

"""
    update_μ!(v::mIHTVariable)

Update the mean `μ` with the linear predictors `BX` and `CZ`. Here `BX` is the
genetic contribution and `CZ` is the non-genetic contribution
"""
function update_μ!(v::mIHTVariable)
    @inbounds @simd for i in eachindex(v.μ)
        v.μ[i] = v.BX[i] + v.CZ[i]
    end
end

"""
    update_resid!(v)

Update residuals `resid = (Y - BX)` where `BX = μ` has already been computed.
"""
function update_resid!(v::mIHTVariable)
    y = v.Y # r × n
    μ = v.μ # r × n
    r = v.resid # r × n
    cv_wts = v.cv_wts # n × 1 cross validation weights
    @inbounds for j in 1:size(y, 2), i in 1:size(y, 1)
        r[i, j] = (y[i, j] - μ[i, j]) * cv_wts[j]
    end
end

"""
    score!(v::mIHTVariable)

Calculates the gradient `Γ(Y - XB)X'` for multivariate Gaussian model. Residuals
were updated in solve_Σ!.
"""
function score!(v::mIHTVariable)
    mul!(v.r_by_n1, v.Γ, v.resid) # r_by_n1 = Γ(Y - BX)
    update_df!(v) # v.df = Γ(Y - BX)X' (genetic gradient)
    mul!(v.df2, v.r_by_n1, Transpose(v.Z)) # v.df2 = Γ(Y - BX)Z' (nongenetic gradient)
end

"""
    update_df!(v::mIHTVariable)

Compute `v.df = Γ(Y - BX)X'` efficiently where `v.r_by_n1 = Γ(Y - BX)`. Here
`X` is `p × n` where `p` is number of SNPs and `n` is number of samples.

Note if `X` is a `SnpLinAlg`, currently only `X * V` with matrix `V` is
    efficiently implemented.
"""
function update_df!(v::mIHTVariable)
    T = eltype(v.X)
    if typeof(v.X) <: Union{Transpose{T, SnpLinAlg{T}}, Adjoint{T, SnpLinAlg{T}}}
        v.n_by_r .= Transpose(v.r_by_n1)
        SnpArrays.mul!(v.p_by_r, v.X, v.n_by_r) # note v.X is Transpose(SnpLinAlg)
        v.df .= Transpose(v.p_by_r)
    elseif typeof(v.X) <: Union{Transpose{T, AbstractMatrix{T}}, Adjoint{T, AbstractMatrix{T}}} # v.X is transposed numeric matrix
        mul!(v.df, v.r_by_n1, v.X.parent)
    else
        mul!(v.df, v.r_by_n1, Transpose(v.X)) # v.X is numeric matrix
    end
end

"""
    _iht_gradstep!(v::mIHTVariable, η::Float)

Computes the gradient step v.b = P_k(β + η∇f(β)) and updates idx and idc.
"""
function _iht_gradstep!(v::mIHTVariable, η::Float)
    # take gradient step: b = b + η ∇f
    BLAS.axpy!(η, v.df, v.B)
    BLAS.axpy!(η, v.df2, v.C)

    # project beta to sparsity
    project_k!(v)
end

function project_k!(v::mIHTVariable)
    # store complete beta [v.b v.c] in full_b, potentially replacing v.C with Inf
    vectorize!(v.full_b, v.B, v.C, v.zkeep)

    # project to sparsity: keeping k + number of nongentic covariates
    project_k!(v.full_b, v.k + v.zkeepn)

    # save model after projection
    unvectorize!(v.full_b, v.B, v.C, v.zkeep)

    # if more than k entries are selected per column, randomly choose k of them
    _choose!(v)

    #recombute support
    update_support!(v.idx, v.B)
    update_support!(v.idc, v.C)

    # make necessary resizing since grad step might include/exclude non-genetic covariates
    check_covariate_supp!(v)
end

"""
    vectorize!(a::AbstractVector, B::AbstractMatrix, C::AbstractMatrix, Ckeep::BitVector)

Without allocations, copies `vec(B)` into `a` and the copies `vec(C)` into
remaining parts of `a`.

`Ckeep` tracks entries of `C` that will not be projected. That is, entries of `a`
corresponding to entries in `C[:, i]` will be filled with `Inf` if `Ckeep[i] = true`.
"""
function vectorize!(a::AbstractVector, B::AbstractMatrix, C::AbstractMatrix, Ckeep::BitVector)
    size(C, 2) == length(Ckeep) ||
        error("size(C, 2) = $(size(C, 2)) and length(Ckeep) = $(length(Ckeep)) are different!")
    i = 1
    @inbounds @simd for j in eachindex(B)
        a[i] = B[j]
        i += 1
    end
    m = typemax(eltype(a))
    for j in 1:size(C, 2)
        if Ckeep[j]
            @inbounds @simd for l in 1:size(C, 1)
                a[i] = m
                i += 1
            end
        else
            @inbounds @simd for l in 1:size(C, 1)
                a[i] = C[l, j]
                i += 1
            end
        end
    end
    return nothing
end

"""
    unvectorize!(a::AbstractVector, B::AbstractMatrix, C::AbstractMatrix, Ckeep::BitVector)

Without allocations, copies the first `length(vec(B))` part of `a` into `B`
and the remaining parts of `a` into `C`.

`Ckeep` tracks entries of `C` that will not be projected. That is, entries in `C[i, :]`
will not be touched if `Ckeep[i] = true`.
"""
function unvectorize!(a::AbstractVector, B::AbstractMatrix, C::AbstractMatrix, Ckeep::BitVector)
    size(C, 2) == length(Ckeep) ||
        error("size(C, 2) = $(size(C, 2)) and length(Ckeep) = $(length(Ckeep)) are different!")
    i = 1
    @inbounds @simd for j in eachindex(B)
        B[j] = a[i]
        i += 1
    end
    @inbounds @simd for j in 1:size(C, 2)
        if !Ckeep[j]
            for l in 1:size(C, 1)
                C[l, j] = a[i]
                i += 1
            end
        end
    end
    return nothing
end

"""
    update_support!(idx::BitVector, b::AbstractMatrix)

Updates `idx` so that `idx[i] = true` if `i`th column of `v.B` contains non-0
values.
"""
function update_support!(idx::BitVector, b::AbstractMatrix{T}) where T
    r, p = size(b)
    fill!(idx, false)
    @inbounds for j in 1:p, i in 1:r
        if b[i, j] != zero(T)
            idx[j] = true
        end
    end
    return nothing
end

"""
    iht_stepsize(v)

Computes the best step size `η = ||∇f||^2_F / tr(X'∇f'Γ∇fX)` where
∇f = Γ(Y - BX)(Y - BX)'. Note the denominator can be rewritten as
`tr(X'∇f'LU∇fX) = ||v||^2_F` where `v = U∇fX`, `L = U'` is cholesky factor
of `Γ`.

Must be careful for cross validation weights in ∇f * X
"""
function iht_stepsize!(v::mIHTVariable{T, M}) where {T <: Float, M}
    # store part of X corresponding to non-zero component of B
    copyto!(v.dfidx, @view(v.df[:, v.idx]))

    # compute numerator of step size
    numer = zero(T)
    @inbounds for i in v.dfidx
        numer += abs2(i)
    end

    # denominator of step size
    denom = zero(T)
    mul!(v.r_by_n1, v.dfidx, v.Xk) # r_by_n1 = ∇f*X
    for j in 1:size(v.Y, 2), i in 1:ntraits(v)
        v.r_by_n1[i, j] *= v.cv_wts[j] # cross validation masking happens here
    end
    cholesky!(Symmetric(v.Γ, :U), Val(true)) # overwrite upper triangular of Γ with U, where LU = Γ, U = L'
    triu!(v.Γ) # set entries below diagonal to 0, so Γ = L'
    mul!(v.r_by_n2, v.Γ, v.r_by_n1) # r_by_n2 = L'*∇f*X
    @inbounds for i in eachindex(v.r_by_n2)
        denom += abs2(v.r_by_n2[i])
    end

    # for bad boundary cases (sometimes, k = 1 in cross validation generates weird η)
    η = numer / denom
    isinf(η) && (η = T(1e-8))
    isnan(η) && (η = T(1e-8))

    return η :: T
end

"""
    project_Σ!(A::AbstractMatrix)

Projects symmetric matrix `A` to the nearest positive semidefinite matrix.

# TODO: efficiency
"""
function project_Σ!(v::mIHTVariable, minλ = 0.01)
    λs, U = eigen(v.Γ)
    for (i, λ) in enumerate(λs)
        λ < minλ && (λs[i] = minλ)
    end
    v.Γ = U * Diagonal(λs) * U'
end

"""
    solve_Σ!(v::mIHTVariable)

Solve for `Σ = 1/n(Y-BX)(Y-BX)'` exactly rather than projecting
"""
function solve_Σ!(v::mIHTVariable)
    update_resid!(v) # v.resid = Y - BX
    mul!(v.r_by_r1, v.resid, Transpose(v.resid)) # r_by_r1 = (Y-BX)(Y-BX)'
    v.r_by_r1 ./= nsamples(v)
    LinearAlgebra.inv!(cholesky!(Symmetric(v.r_by_r1, :U))) # r_by_r1 = (1/n(Y-BX)(Y-BX)')^{-1}
    copyto!(v.Γ, v.r_by_r1)
end

"""
    check_covariate_supp!(v::mIHTVariable)

Possibly rescales `v.Xk` and `v.dfidx`, which needs to happen when non-genetic
covariates get included/excluded between different iterations
"""
function check_covariate_supp!(v::mIHTVariable{T, M}) where {T <: Float, M}
    n, r = size(v.X, 2), ntraits(v)
    nzidx = sum(v.idx)
    if nzidx != size(v.Xk, 1)
        v.Xk = zeros(T, nzidx, n)
        v.dfidx = zeros(T, r, nzidx) # TODO ElasticArrays.jl
    end
    # @inbounds for i in eachindex(v.zkeep)
    #     v.zkeep[i] && !v.idc[i] && error("A non-genetic covariate was accidentally set to 0! Shouldn't happen!")
    # end
end

"""
    _choose!(v::mIHTVariable)

When `B` has `≥k` non-zero entries in a row, randomly choose `k` among them.

Note: It is possible to have `≥k` non-zero entries after projection due to
numerical errors.
"""
function _choose!(v::mIHTVariable)
    sparsity = v.k + v.zkeepn
    B_nz = 0
    C_nz = 0
    B_nz_idx = Int[]
    C_nz_idx = Tuple{Int, Int}[]
    # find position of non-zero beta
    for i in eachindex(v.B)
        if v.B[i] != 0
            B_nz += 1
            push!(B_nz_idx, i)
        end
    end
    for j in 1:size(v.C, 2)
        v.zkeep[j] && continue
        for i in 1:size(v.C, 1)
            if v.C[i, j] != 0
                C_nz += 1
                push!(C_nz_idx, (i, j))
            end
        end
    end

    # if more non-zero than expected, randomly set entries (starting in B) to 0
    excess = B_nz + C_nz - sparsity
    if excess > 0
        shuffle!(B_nz_idx)
        shuffle!(C_nz_idx)
        for i in 1:excess
            if B_nz > 0
                v.B[B_nz_idx[i]] = 0
                B_nz -= 1
            else
                (ii, jj) = C_nz_idx[i]
                v.C[ii, jj] = 0
                C_nz -= 1
            end
        end
    end
    empty!(B_nz_idx)
    empty!(C_nz_idx)
end

"""
Function that saves variables that need to be updated each iteration
"""
function save_prev!(v::mIHTVariable, cur_logl, best_logl)
    copyto!(v.B0, v.B)     # B0 = B
    copyto!(v.idx0, v.idx) # idx0 = idx
    copyto!(v.idc0, v.idc) # idc0 = idc
    copyto!(v.C0, v.C)     # C0 = C
    copyto!(v.Γ0, v.Γ)     # Γ0 = Γ
    if cur_logl > best_logl
        copyto!(v.best_B, v.B) # best_B = B
        copyto!(v.best_C, v.C) # best_C = C
    end
    return max(cur_logl, best_logl)
end

checky(y::AbstractMatrix, d::MvNormal) = nothing

"""
When initializing the IHT algorithm, take `k` largest elements in magnitude of
the score as nonzero components of `v.B`. This function set v.idx = 1 for
those indices. `v.μ`, `v.df`, `v.df2`, and possibly `v.BX` are also initialized.
"""
function init_iht_indices!(v::mIHTVariable, init_beta::Bool, cv_idx::BitVector, verbose::Bool)
    fill!(v.B, 0)
    fill!(v.B0, 0)
    fill!(v.best_B, 0)
    fill!(v.BX, 0)
    fill!(v.Xk, 0)
    fill!(v.idx, false)
    fill!(v.idx0, false)
    copyto!(v.idc, v.zkeep)
    copyto!(v.idc0, v.zkeep)
    fill!(v.resid, 0)
    fill!(v.df, 0)
    fill!(v.df2, 0)
    fill!(v.dfidx, 0)
    fill!(v.C, 0)
    fill!(v.C0, 0)
    fill!(v.best_C, 0)
    fill!(v.CZ, 0)
    fill!(v.μ, 0)
    fill!(v.Γ, 0)
    fill!(v.Γ0, 0)
    for i in 1:ntraits(v)
        v.Γ[i, i] = 1
        v.Γ0[i, i] = 1
    end
    fill!(v.cv_wts, 0)
    v.cv_wts[cv_idx] .= 1
    fill!(v.full_b, 0)
    fill!(v.r_by_r1, 0)
    fill!(v.r_by_r2, 0)
    fill!(v.r_by_n1, 0)
    fill!(v.r_by_n2, 0)
    fill!(v.n_by_r, 0)
    fill!(v.p_by_r, 0)
    fill!(v.k_by_r, 0)
    fill!(v.k_by_k, 0)

    # initialize intercept to mean of each trait
    nz_samples = count(!iszero, v.cv_wts) # for cross validation masking
    for i in 1:ntraits(v)
        ybar = zero(eltype(v.C))
        @inbounds for j in 1:size(v.Y, 2)
            ybar += v.Y[i, j] * v.cv_wts[j]
        end
        v.C[i, 1] = ybar / nz_samples
    end
    mul!(v.CZ, v.C, v.Z)

    if init_beta
        initialize_beta!(v, cv_idx, verbose)
        project_k!(v)
        update_xb!(v)
    end

    # update mean and compute score (gradient)
    update_μ!(v)
    update_resid!(v)
    score!(v)

    if !init_beta
        # first `k` non-zero entries in each β are chosen based on largest gradient
        vectorize!(v.full_b, v.df, v.df2, v.zkeep)
        project_k!(v.full_b, v.k + v.zkeepn) # project k + number of nongentic covariates to keep
        unvectorize!(v.full_b, v.df, v.df2, v.zkeep)

        # compute support based on largest gradient
        update_support!(v.idx, v.df)
        update_support!(v.idc, v.df2) # should we initialize v.idc to all trues?
    end

    # make necessary resizing when necessary
    check_covariate_supp!(v)

    # store relevant components of x for first iteration
    copyto!(v.Xk, @view(v.X[v.idx, :]))
end

function check_convergence(v::mIHTVariable)
    the_norm = max(chebyshev(v.B, v.B0), chebyshev(v.C, v.C0)) #max(abs(x - y))
    scaled_norm = the_norm / (max(norm(v.B0, Inf), norm(v.C0, Inf)) + 1.0)
    return scaled_norm
end

function backtrack!(v::mIHTVariable, η::Float)
    # recompute gradient step with the new (halved) step size η
    copyto!(v.B, v.B0)
    copyto!(v.C, v.C0)
    copyto!(v.Γ, v.Γ0)
    _iht_gradstep!(v, η)

    # recompute BX and ZC, μ (includes genetic + nongenetic component), and Γ
    update_xb!(v)
    update_μ!(v)
    solve_Σ!(v)

    return loglikelihood(v)
end

"""
    is_multivariate(y::AbstractVecOrMat)

Returns true if response `y` can be modeled by a multivariate Gaussian
distributions. Currently simply checks whether there is >1 trait
"""
function is_multivariate(y::AbstractVecOrMat)
    size(y, 1) > 1 && size(y, 2) > 1
end

function save_best_model!(v::mIHTVariable)
    # compute η = BX with the best estimated model
    copyto!(v.B, v.best_B)
    copyto!(v.C, v.best_C)
    update_support!(v.idx, v.B)
    update_support!(v.idc, v.C)
    check_covariate_supp!(v)
    update_xb!(v)

    # update estimated mean μ with genotype predictors
    update_μ!(v)
end

function save_last_model!(v::mIHTVariable)
    # compute η = BX with the estimated model from last IHT iteration
    update_support!(v.idx, v.B)
    update_support!(v.idc, v.C)
    check_covariate_supp!(v)
    update_xb!(v)

    # update estimated mean μ with genotype predictors
    update_μ!(v)
end

"""
    initialize_beta!(v, cv_wts, verbose)

Initialze beta to univariate regression values. That is, `β[i, j]` is set to the estimated
beta with `y[j, :]` as response, and `x[:, i]` with an intercept term as covariate.

Progress will be displayed if `verbose=true`

Note: this function assumes quantitative (Gaussian) phenotypes.
"""
function initialize_beta!(
    v::mIHTVariable,
    cv_wts::BitVector,# cross validation weights; 1 = sample is present, 0 = not present
    verbose::Bool
    )
    y, x, z, B, C, T = v.Y, v.X, v.Z, v.B, v.C, eltype(v.B)
    xtx_store = [zeros(T, 2, 2) for _ in 1:Threads.nthreads()]
    xty_store = [zeros(T, 2) for _ in 1:Threads.nthreads()]
    xstore = [zeros(T, sum(cv_wts)) for _ in 1:Threads.nthreads()]
    ystore = zeros(T, sum(cv_wts))
    pmeter = verbose ? Progress(ntraits(v)*(nsnps(v) + ncovariates(v)), 5,
        "Initializing β to univariate regression values...") : nothing
    @inbounds for j in 1:ntraits(v) # loop over each y
        copyto!(ystore, @view(y[j, cv_wts]))
        # genetic covariates
        Threads.@threads for i in 1:nsnps(v)
            id = Threads.threadid()
            copyto!(xstore[id], @view(x[i, cv_wts]))
            linreg!(xstore[id], ystore, xtx_store[id], xty_store[id])
            B[j, i] = xty_store[id][2]
            verbose && next!(pmeter) # update progress
        end
        # non genetic covariates
        Threads.@threads for i in 1:ncovariates(v)
            id = Threads.threadid()
            copyto!(xstore[id], @view(z[i, cv_wts]))
            linreg!(xstore[id], ystore, xtx_store[id], xty_store[id])
            C[j, i] = xty_store[id][2]
            verbose && next!(pmeter) # update progress
        end
    end
    clamp!(C, -2, 2)
    clamp!(B, -2, 2)
    copyto!(v.B0, v.B)
    copyto!(v.C0, v.C)
end

"""
    debias!(v::mIHTVariable)

Solves the multivariate linear regression `Y = BX + E` by `B̂ = inv(X'X) X'Y` on the
support set of `B`. Since `B` is sparse, this is a low dimensional problem, and
the solution is unique.

Note: since `X` and `Y` are transposed in memory, we actually have `B̂ = inv(XX')XY'`
"""
function debias!(v::mIHTVariable{T, M}) where {T <: Float, M}
    # first rescale matrix dimension if needed
    supp_size = sum(v.idx)
    if supp_size != size(v.k_by_r, 1)
        v.k_by_r = Matrix{T}(undef, supp_size, size(v.k_by_r, 2))
        v.k_by_k = Matrix{T}(undef, supp_size, supp_size)
    end

    # try debiasing: B̂ = inv(XX')XY'. Do nothing if it fals
    mul!(v.k_by_r, v.Xk, Transpose(v.Y))
    mul!(v.k_by_k, v.Xk, Transpose(v.Xk))
    try
        ldiv!(cholesky!(Symmetric(v.k_by_k, :U)), v.k_by_r)
    catch
        return nothing
    end
    v.B[:, v.idx] .= v.k_by_r'

    # ensure B is k-sparse
    project_k!(v)
end


text:
#=
Problem: Divider
Max. score: 20

Shopee has N software engineers.
Shopee accommodates them by arranging N tables on a 1D plane.
However, many people raise concerns that the work environment is too noisy.
To mitigate this issue, Shopee decides to group the engineers into K groups.
A group is a non-overlapping segment of contiguous engineers.
Shopee will then put dividers between the groups.

Shopee wants to minimize the total noise value.
Please help Shopee to find the minimum total noise value possible.
=#

# Julia program to calculate sum of
# 5 integers obtained from console/user input

#Variables
line1 = ""
line2 = ""
allEng = []

#Functions
#Calculate noise in a group of engineers
function calculateNoise(l,r, arrEngineers)
    noise = sum(arrEngineers[l:r]) * (r-l+1)
    #println("Noise between index ",l," and ",r," : ",noise)
    return noise
end

#Calculate noise in all the groups
function calculateAllNoises(N,K,arrDivide,arrEngineers)
    start = 0
    last = N
    sum = 0
    for i in 1:K-1
        sum += calculateNoise(arrDivide[i],arrDivide[i+1],arrEngineers)
    end
    return sum
end

function divideArray(arr, n, k)

    # Dp to store the values
    dp = [[0 for i in range(500)]
             for i in range(500)]

    k -= 1
    # Fill up the dp table
    for i in -1:n-1
        for j in 0:k + 1

            # Intitilize maximum value
            dp[i][j] = 199999999999999999999999999

            # Max element and the summ
            min = 999999999999999999999999999999
            summ = 0

            # Run a loop from i to n
            for l in i:n

                # Find the maximum number
                # from i to l and the summ
                # from i to l
                min_ = calculateNoise(min_, arr[l],arr)
                summ += arr[l]

                # Find the summ of difference
                # of every element with the
                # maximum element
                diff = (l - i + 1) * min_ - summ

                # If the array can be divided
                if (j > 0)
                    dp[i][j]= min(dp[i][j], diff +
                                  dp[l + 1][j - 1])
                else
                    dp[i][j] = diff
                end
            end
        end
    end
    # Returns the minimum summ
    # in K parts
    return dp[0][k]
end

# Prompt to enter 1st line
#println("Enter N and K respectively on same line: ")

line1 = readline()
inputArr = split(line1, " ")
N = parse(Int64, inputArr[1])
K = parse(Int64, inputArr[2])
println("N:", N)
#rintln("K:", K)

# Prompt to enter 2nd line
#println("Enter ",N," numbers in a line")

line2 = readline()
inputArr = split(line2, " ")
# Taking Input from user
for i in 1:N
   push!(allEng,parse(Int64, inputArr[i]))
end
#println("Engineers: ",line2)

#Run noise check
noise = calculateNoise(2,3,allEng)
noise = divideArray(allEng, N,K)
println(noise)

ext: jl

max_stars_repo_path: 03 Programming Contest/divider.jl

max_issues_repo_path: 03 Programming Contest/divider.jl

max_forks_repo_path: 03 Programming Contest/divider.jl


text:
using Glove
using Test
using LinearAlgebra

# Anything <: AbstractFloat, just here for easy use
S = Float32

# Test from python Glove implementation.
# https://github.com/maciejkula/glove-python/tree/master/glove/tests
corpus1 = split(read("data/corpus_test1.txt", String))
table = LookupTable(corpus1)
@test table.word2id == Dict("a" => 1, "naive" => 2, "fox" => 3)
@test table.id2word == Dict(1 => "a", 2 => "naive",  3 => "fox")
@test length(table) == 3
@test table["fox"] == 3
@test table[1] == "a"

# create the co-occurence matrix and table in 1 pass on the corpus.
codict = weightedsums(S, table, corpus1)
@test codict[(1,1)] == 0.0
@test codict[(1,2)] == 1.0
@test codict[(1,3)] == 0.5

@test codict[(2,1)] == 1.0
@test codict[(2,2)] == 0.0
@test codict[(2,3)] == 1.0

@test codict[(3,1)] == 0.5
@test codict[(3,2)] == 1.0
@test codict[(3,3)] == 0.0

corpus2 = split(read("data/corpus_test2.txt", String))
table = LookupTable(corpus2)
codict = weightedsums(S, table, corpus2)

cokeys = collect(keys(codict))
covals = collect(values(codict))

n = length(codict)
covec = CooccurenceVector{S}(undef, n)

@inbounds for i=1:n
  covec[i] = (cokeys[i][1], cokeys[i][2], covals[i])
end

vocabsize = length(table)
vecsize = 10

model = Model(S, vocabsize, vecsize)

@test typeof(model) == Glove.Model{S}
@test eltype(model.W_main) == S
@test typeof(covec) == CooccurenceVector{S}
@test eltype(covec) == Tuple{Int, Int, S}

# The corpus is very, very small.
# So even with a large amount of iterations
# this test may fail.
adagrad!(S, model, covec, 500)

# make sure everything is still the same type
@test eltype(model.W_main) == S

# similar words
function similar_words(
  M::Matrix{T},
  table::LookupTable,
  word::S;
  n::Int=10) where {T, S<:Token}

    c_id = table[word]

    dists = vec(M[:, c_id]' * M) / norm(M[:, c_id]) / norm(M, 1)

    sorted_ids = sortperm(dists, rev=true)[1:n+1]
    sim_words = Token[]

    for id = sorted_ids
        if c_id == id
            continue
        end
        word = table[id]
        push!(sim_words, word)
    end
    sim_words
end

# paper recommends adding the main and context matrices.
M = model.W_main + model.W_ctx
top_words = similar_words(M, table, "trees", n=10)[1:3]
@test in("graph", top_words)
top_words = similar_words(M, table, "graph", n=10)[1:3]
@test in("trees", top_words)


text:
RA_COLOR = Dict(:cw => :green, :ccw => :blue, :none => :gray)
RA_SHAPE = Dict(:asc => :utriangle, :desc => :dtriangle, :none => :circle)
FORMAT = (grid=true, gridalpha=0.2, minorgrid=true, minorgridalpha=0.1, legend=:none)

function parse_connect(log, connect)
    if connect == :none
        []
    elseif connect == :all
        subsets(1:log["n"], 2)
    elseif connect isa Tuple{Z, Z}
        [connect]
    elseif connect isa Vector{Tuple{Z, Z}}
        connect
    else
        throw(ArgumentError("keyword `connect` does not support value $connect."))
    end
end

function add_horizontal(log::Dict{String, Any}; connect=:all)
	p = plot(; FORMAT..., title="Horizontal Position")
	ss = []
	for i in 1:log["n"]
		s = reduce(vcat, s' for s in log["ac_$i"])
		push!(ss, s)
		plot!(p, s[:,1], s[:,2]; color=i, aspect_ratio=:equal, FORMAT...)
		scatter!(p, [s[1,1]], [s[1,2]]; ms=3, mc=:black, shape=:star, FORMAT...)

		# advisories
		ras = log["cas_$i"]
		for (t, ra) in enumerate(ras)
			if any(ra .!= :none)
				color = RA_COLOR[ra[1]]
				shape = RA_SHAPE[ra[2]]
				scatter!(p, s[t,1:1], s[t,2:2]; ms=3, msw=0, mc=color, shape=shape, FORMAT...)
			end
		end
	end

    # add inter-aircraft connections
    for (i, j) in parse_connect(log, connect)
	    for t in 1:size(ss[1], 1)
            xs = (s -> s[t,1]).(ss[[i,j]])
            ys = (s -> s[t,2]).(ss[[i,j]])
		    plot!(p, xs, ys; lc=:black, ls=:dot, lw=0.5, ma=0.5, FORMAT...)
        end
	end

	return p
end

function add_vertical(log::Dict{String, Any})
	p = plot(title="Vertical Position")
    ts = log["t"]

	for i in 1:log["n"]
		s = reduce(vcat, s' for s in log["ac_$i"])
		plot!(p, ts, s[:,3]; color=i, FORMAT...)

		# advisories
		ras = log["cas_$i"]
		for (t, ra) in enumerate(ras)
			if any(ra .!= :none)
				color = RA_COLOR[ra[1]]
				shape = RA_SHAPE[ra[2]]
				scatter!(p, ts[t:t], s[t,3:3]; ms=3, msw=0, mc=color, shape=shape, FORMAT...)
			end
		end
	end
	return p
end

function add_separation(log::Dict{String, Any})
	return plot(log["t"], log["d"]; title="Separation", FORMAT...)
end

function add_velocity(log::Dict{String, Any})
	p = plot(title="Ground Velocity")
	for i in 1:log["n"]
		s = reduce(vcat, s' for s in log["ac_$i"])
		plot!(p, log["t"], sqrt.(s[:,4].^2 + s[:,5].^2); FORMAT...)
	end
	return p
end

function visualize(log::Dict{String, Any}; filename::String="", kwargs...)
    @assert !isempty(log)

	ps = [
        add_horizontal(log; kwargs...),
        add_vertical(log),
        add_separation(log),
        add_velocity(log)
    ]

	p = plot(ps..., layout=(2,2), size=(900,500))
    isempty(filename) || savefig(p, filename)
	return p
end

visualize(sim::Simulator; kwargs...) = visualize(sim.log; kwargs...)


text:
using Hecke

using Test

using LinearAlgebra

Hecke.assertions(true)

## Test the test_module on Windows

x = Hecke._adjust_path("GrpAb/Elem")
y = joinpath(Hecke.pkgdir, "test", "$x.jl")
@test isfile(y)

k, a = quadratic_field(5)
@test fmpz(1) - a == -(a - 1)
@test 1 - a == -(a - 1)

push!(Base.LOAD_PATH, "@v#.#")

using Random
using RandomExtensions

const rng = MersenneTwister()
const rand_seed = rand(UInt128)

# tests if rand(rng, args...) gives reproducible results
function reproducible(args...)
  Random.seed!(rng)
  x = rand(rng, args...)
  Random.seed!(rng, rng.seed)
  x == rand(rng, args...)
end

try
  using GAP
  @time include("FieldFactory.jl")
catch e
  if !(isa(e, ArgumentError))
    rethrow(e)
  else
    println("using GAP failed. Not running FieldFactory tests")
  end
end

@time include("NumField.jl")
@time include("AlgAss.jl")
@time include("AlgAssAbsOrd.jl")
@time include("AlgAssRelOrd.jl")
@time include("EllCrv.jl")
@time include("GrpAb.jl")
@time include("Grp.jl")
@time include("LinearAlgebra.jl")
@time include("Map.jl")
@time include("Misc.jl")
@time include("NfAbs.jl")
@time include("NfOrd.jl")
@time include("NfRel.jl")
@time include("RCF.jl")
@time include("Examples.jl")
@time include("Sparse.jl")
@time include("QuadForm.jl")
@time include("LocalField.jl")

#try
#  using Polymake
#  @time include("AlgAssRelOrd/Eichler.jl")
#catch e
#  if !(isa(e, ArgumentError))
#    rethrow(e)
#  else
#    println("using Polymake failed. Not running sophisticated norm equation tests")
#  end
#end


text:
function rsaencode(clearmsg::AbstractString, nmod::Integer, expub::Integer)
    bytes = parse(BigInt, "0x" * bytes2hex(collect(UInt8, clearmsg)))
    return powermod(bytes, expub, nmod)
end

function rsadecode(cryptmsg::Integer, nmod::Integer, dsecr::Integer)
    decoded = powermod(encoded, dsecr, nmod)
    return join(Char.(hex2bytes(hex(decoded))))
end

msg = "Rosetta Code."
nmod = big"9516311845790656153499716760847001433441357"
expub = 65537
dsecr = big"5617843187844953170308463622230283376298685"

encoded = rsaencode(msg, nmod, expub)
decoded = rsadecode(encoded, nmod, dsecr)
println("\n# $msg\n -> ENCODED: $encoded\n -> DECODED: $decoded")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using BrainWaveExperimental

fToProcess = ["BrainWaveExperimental.jl", "findABRPeaks.jl", "utils.jl"]
for fName in fToProcess
#fName = fToProcess[1]
    fIn = open(string("../src/", fName), "r")
    fOut = open(string("test_", fName), "w")
    lns = readlines(fIn)
    idxStart = (Int)[]
    idxStop = (Int)[]
    for i=1:length(lns)
        if lns[i] == "```julia"
            push!(idxStart, i+1)
        elseif lns[i] == "```"
            push!(idxStop, i-1)
        end
    end
    for i=1:length(idxStart)
        theseLines = lns[idxStart[i]:idxStop[i]]
        for currLine in theseLines
            write(fOut, string(currLine, "\n"))
        end
    end
    close(fIn); close(fOut)
end

for fName in fToProcess
    fNameTest = string("test_", fName)
    include(fNameTest)
end


text:
function main(window)
    push!(window.assets,"layout2")
    A = vbox("Hello, world!"|>
             pad(1em)|>
             fillcolor("#006699"),
             "Howdy!"|>
             pad(1em)|>
             fillcolor("#669933")
             )
    pages([A])
end


text:
import Compat: ∘                # FIXME remove in v0.6
using ArgCheck
using Parameters

export crra_u, CrraUtility, crra_u′

"""
CRRA/isoelastic utility, with risk aversion parameter `σ`.
"""
@inline function crra_u(c, σ)
    omσ = one(σ) - σ
    omσ == zero(omσ) ? log(c) : (c^omσ - 1)/omσ
end

"""
CRRA/isoelastic utility, with risk aversion parameter `σ`. Returns a
function, whcih can be evaluated at `c`. See also `crra_u`.
"""
immutable CrraUtility{T}
    σ::T
end

(u::CrraUtility)(c) = crra_u(c, u.σ)

# FIXME at σ=1, the continuous extension of the derivative in σ is
# -log(c)^2/2, need to implement this, with tests

"""
Return a transformation for the marginal CRRA/isoelastic utlity, with
risk aversion parameter `σ`. Use `Base.inv` to obtain inverse.
"""
function crra_u′(σ)
    @argcheck σ ≥ 0
    Power(-σ)
end


text:
# identity.jl
# test

using Test: @test, @testset
using LinearAlgebra: I

@testset "identity" begin
	X = rand(2,3,4) # AbstractArray
	@test I * X === X
	@test (5I) * X == 5*X
end


text:
using Compat
using GeoJSON
using GridOperations
using LibGEOS
using PeatUtils
using TiledIteration
using XYZ
using Glob
using Lumberjack

const overlap = 0.  # width of overlap [m]

# general settings
const high_res = 1.0  # cellsize for water grids & and high res dtm [m]
const nodata = -9999.0

function xyvec(cloud::XYZ.Cloud, sampleratio::Float64)
    n = length(cloud.positions)
    step = round(Int, cld(1, sampleratio))
    # have minimum of points needed
    if n < 10000 | fld(n, step) < 1000
        step = 1
    end
    idx_samples = 1:step:n
    nsample = length(idx_samples)
    A = Vector{Vector{Float64}}() # empy vector of vectors
    sizehint!(A, nsample)
    for i in idx_samples
        p = cloud.positions[i]
        push!(A, [p[1], p[2]])
    end
    A
end

# simplify raster from cloud
function cloud2arr(cloud::XYZ.Cloud, tile_bbox::XYZ.BoundingBox, cellsize::Float64; pointfilter=nothing, reducer=XYZ.reducer_medz, interpolate=true)
    index = XYZ.define_raster(cloud, tile_bbox, overlap, cellsize; snapgrid=25.0, epsg=epsg, pointfilter=pointfilter)
    raster = XYZ.rasterize(cloud, index; reducer=reducer, min_dens = 0)[:,:,1]
    if interpolate
        raster = GridOperations.interp_missing!(copy(raster), eltype(raster)(nodata); mode="kriging")
    end
    raster, index
end

# simplify tif creation
function cloud2tif(folder::String, fn::String, cloud::XYZ.Cloud, tile_bbox::XYZ.BoundingBox, cellsize::Float64, filter, reducer; interpolate=false)
    raster, index = cloud2arr(cloud, tile_bbox, cellsize, pointfilter=filter, reducer=reducer, interpolate=interpolate)
    XYZ.grid2tif(folder, fn, index, raster; nodata=nodata)
end

"""pipeline to classify cloud and create dtm"""
function lidar_pipeline(
        filename::String,
        epsg::Int,
        out_dir::String,
        gf_radius::Float64,
        gf_slope::Float64,
        gf_dh_max::Float64,
        gf_dh_min::Float64,
        )

    # file admin
    fileid = basename(splitext(filename)[1])
    in_dir = dirname(filename)
    filen = fileid

    out_dir = joinpath(out_dir, fileid)
    isdir(out_dir) || mkpath(out_dir)
    info("Ouput directory is $(out_dir)")

    # Set up Logging
    add_truck(LumberjackTruck(joinpath(out_dir, "pipeline.log")), "pipelinelogger")
    info("Processing $(filen)")
    info("with ground parameters r:$(gf_radius) \t s:$(gf_slope) \t min:$(gf_dh_min) \t max:$(gf_dh_max)")

    # read las
    info("--> read $(fileid)")
    info("--> write results to $(out_dir)")
    cloud, header = XYZ.read_pointcloud(filename)
    info("  $(cloud)")

    # reset classification
    XYZ.reset_class!(cloud)

    tile_bbox = XYZ.BoundingBox(xmin=header.x_min, ymin=header.y_min,
        xmax=header.x_max, ymax=header.y_max)

    ## CLASSIFY OUTLIERS
    info("--> 1.) classify outliers")
    # classify low and high outliers
    # evaluate points per in <cellsize> [m] cells
    # low outliers if a jump larger than <dz> [m] in lowest n <max_outliers> [-] points
    # high outliers if points more than <max_height> [m] above lowest non-outlier point
        XYZ.classify_outliers!(cloud;
        cellsize = 100.0,   # cellsize used for evaluation of points
        dz = 1.0,           # threshold in vertical distance between low points
        max_outliers = 15,  # max number of outliers in one cell = low  points to be evaluated
        max_height=100.0)   # threshold in vertical distance for high points

    ## CLASSIFY GROUND with Zhang algorithm. use only non-water, non-outlier, last return points
    info("--> 1.5) intermediate rasters")

    # make non-filled raster index for high resolution cells with non-outliers and last return points
    r = XYZ.define_raster(cloud, tile_bbox, overlap, high_res; pointfilter = XYZ.unclassified_lastreturn, epsg = epsg)
    if !(r.nrow >= 2 && r.ncol >= 2)
        info("Defined raster is too small to be classified, skipping")
        @show r
        return
    end

    boundarymask_r = trues(r.nrow, r.ncol)
    smallboundarymask_r = trues(r.nrow, r.ncol)

    zmin = XYZ.rasterize(cloud, r; reducer = XYZ.reducer_minz,  pointfilter=nothing)[:,:,1]
    XYZ.grid2tif(out_dir, "$(filen)_zmin.tif", r, zmin; nodata=nodata)

    # create las mask and dropouts & vegetation binary images
    haspoints = Array{Bool}(XYZ.rasterize(cloud, r; reducer = XYZ.reducer_count)[:,:,1] .> 0)
    las_mask = GridOperations.create_mask(haspoints)
    dropouts = Array{Bool}(las_mask .& broadcast(~, haspoints))  # nodata cells within las mask

    # zmin interpolated surface
    r_filled = XYZ.inpaint_missings_nn(r, cloud)[1]
    zmin_filled = XYZ.rasterize(cloud, r_filled; reducer = XYZ.reducer_minz,  pointfilter=nothing)[:,:,1]
    XYZ.grid2tif(out_dir, "$(filen)_zmin_nn.tif", r, zmin_filled; nodata=nodata)

    ## CLASSIFY GROUND with Zhang algorithm. use only non-water, non-outlier, last return points
    info("--> 2.) classify terrain (pmf)")

    # apply zhang to grid
    zmax_pmf, flags = GridOperations.pmf_filter(zmin_filled, boundarymask_r, gf_radius, gf_slope, gf_dh_max, gf_dh_min, r.cellsize)
    XYZ.grid2tif(out_dir, "$(filen)_zmax_pmf.tif", r, zmax_pmf; nodata=nodata)
    bin_vegetation = Array{Bool}(flags .> 0) .& broadcast(~, dropouts)

    # use maxz grid for classification
    XYZ.classify_below_surface!(cloud, r, zmax_pmf, 2)
    XYZ.classify_mask!(cloud, r, smallboundarymask_r, 28)  # unused class for edge effects

    # filter zmin grid
    zmin_filtered = copy(zmin)
    zmin_filtered[((bin_vegetation .| dropouts) .| .!smallboundarymask_r)] = nodata
    XYZ.grid2tif(out_dir, "$(filen)_zmin_filtered.tif", r, zmin_filtered; nodata=nodata)

    nothing
end


text:
"""
StochasticGene
"""
module StochasticGene

using Distributions
using StatsBase
using Statistics
using DelimitedFiles
using Dates
using Distributed
using LinearAlgebra
using PyPlot
using SparseArrays
using DifferentialEquations
using LSODA

export
    Model,
    GModel,
    GRModel,
    GRSModel,
    RNALiveCellSplice,
    ExperimentalData,
    HistogramData,
    RNALiveCell,
    runmodel,
    logprior,
    priordistribution,
    likelihoodtuple,
    crossentropy,
    likelihoodfn,
    readparams,
    readFISH,
    readLCPDF,
    runparallel,
    metropolis_hastings,
    sample,
    waicupdate!,
    mhStep!,
    proposal,
    initialproposal,
    sigmalogNormal,
    loglikelihood,
    MHOPtions,
    offtimeCDF,
    telegraphPDF

### Source files

# Type system
include("common.jl")

# commonly used functions
include("utilities.jl")

# Metropolis Hastings MCMC for computing posterior distributions of model parameters
include("metropolis_hastings.jl")

# Input output functions
include("io.jl")

# Chemical master equation solutions of stochastic models for likelihood functions in fitting algorithms
include("chemical_master.jl")

# Transition rate matrices of stochastic models defining master equations
include("transition_rate_matrices.jl")

# Probability distributions by direct simulation of stochastic models using Gillespie and Gibson-Bruck algorithms
include("ClassicTelegraph.jl")
include("simulator.jl")

# functions specific for Gene Trap experiments of Wan et al.
include("genetrap.jl")

# functions for scRNA and FISH experiments
include("rna.jl")

# include("/Users/carsonc/Dropbox/Larson/GeneTrap_analysis/code/GillespieSimulatorExInt.jl")
# include("scRNA.jl")
# include("tff1.jl")

"""
A Julia module for simulation and Bayesian inference of parameters of stochastic models of gene transcription.

API Overview:

"""

end #Module

max_stars_repo_name: nih-niddk-mbs/StochasticGene

max_issues_repo_name: nih-niddk-mbs/StochasticGene

max_forks_repo_name: nih-niddk-mbs/StochasticGene


text:
using Plexifications
using Base.Test: @test, @test_throws

@test_throws ErrorException begin
    inv(Prplex(1, 1))
end

@test begin
    iszerodivisor(Prplex(1,1))
end

@test begin
    !iszerodivisor(Prplex(0,1))
end

@test begin
    length(unreal(1)) == 0
end

@test begin
    a = Plexifications.randomBigFloat()
    isreal(Prplex(a))
end

@test begin
    a = 1
    !isreal(Prplex(a, a))
end

@test begin
    z = random(Prplex{Int})
    z == +(z)
end

@test begin
    a = Plexifications.randomBigFloat()
    real(Prplex(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, Prplex(1, 2))
    l = String(take!(io))
    r = "[1: 1, R: 2]"
    l == r
end

@test begin
    x = random(Prplex{BigFloat})
    y = random(Prplex{BigFloat})
    x != y
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(Prplex{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(Prplex{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    z = random(Prplex{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    z = random(Prplex{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    barR(barR(x)) == x
end

@test begin
    x = random(Prplex{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    l = barR(x * y)
    r = barR(y) * barR(x)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    l = barR(x + y)
    r = barR(x) + barR(y)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    l = barR(x - y)
    r = barR(x) - barR(y)
    l == r
end

@test begin
    x = random(Prplex{Rational{BigInt}})
    y = random(Prplex{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Prplex{Rational{BigInt}})
    y = random(Prplex{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Prplex{Rational{BigInt}})
    y = random(Prplex{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    z = random(Prplex{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(Prplex{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    z = random(Prplex{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    z = random(Prplex{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    z = random(Prplex{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    abs2(x) > big(0) || abs2(x) < big(0) || abs2(x) == big(0)
end

@test begin
    x = random(Prplex{BigInt})
    y = random(Prplex{BigInt})
    l = abs2(x * y)
    r = abs2(x) * abs2(y)
    l == r
end

@test_throws ErrorException begin
    inv(Prplex(0, 0.0))
end

@test_throws ErrorException begin
    random(Prplex{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Prplex{Int})
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(Prplex{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(Prplex{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Prplex{BigInt})
    z == z + zero(Prplex{BigInt})
end

@test begin
    z = random(Prplex{BigInt})
    z == z * one(Prplex{BigInt})
end

@test begin
    z = random(Prplex{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Prplex{BigInt})
    z == z * one(z)
end


text:
"""
    Configuration(configfile::String)

Parses configuration file `configfile`.
"""
struct Configuration
    equation_name::String
    scenario_name::String

    end_time::Float64
    timeintegrator_name::String
    courant::Float64

    plot_start::Float64
    plot_step::Float64
    grid_elements::Int64
    physicalsize::Array{Float64,1}
    cellsize::Array{Float64,1}

    function Configuration(configfile::String)
        config = YAML.load(open(configfile))
        equation_name = config["equation"]["equation"]
        scenario_name = config["equation"]["scenario"]

        end_time = config["simulation"]["end_time"]
        timeintegrator_name = config["solver"]["timeintegrator"]
        courant = config["solver"]["courant"]

        plot_start = config["output"]["start"]
        plot_step = config["output"]["step"]
        grid_elements = config["simulation"]["grid_elements"]
        physicalsize_1d = config["simulation"]["grid_size"]
        physicalsize = [physicalsize_1d, physicalsize_1d]
        cellsize = physicalsize ./ grid_elements

        new(
            equation_name,
            scenario_name,
            end_time,
            timeintegrator_name,
            courant,
            plot_start,
            plot_step,
            grid_elements,
            physicalsize,
            cellsize
        )
    end
end


text:
"""
		Implementation of Hyperspherical Variational Auto-Encoders

		Original paper: https://arxiv.org/abs/1804.00891

		SVAEtwocaps(q,g,zdim,hue,μzfromhidden,κzfromhidden)

		q --- encoder - in this case encoder only encodes from input to a hidden
						layer which is then transformed into parameters for the latent
						layer by `μzfromhidden` and `κzfromhidden` functions
		g --- decoder
		zdim --- dimension of the latent space
		hue --- Hyperspherical Uniform Entropy that is part of the KL divergence but depends only on dimensionality so can be computed in constructor
		μzfromhidden --- function that transforms the hidden layer to μ parameter of the latent layer by normalization
		κzfromhidden --- transforms hidden layer to κ parameter of the latent layer using softplus since κ is a positive scalar
"""
mutable struct SVAEtwocaps{V<:Val} <: SVAE
	q
	g
	zdim
	hue
	μzfromhidden
	κzfromhidden
	priorμ
	priorκ
	variant::V

	"""
	SVAEtwocaps(q, g, hdim, zdim, T) Constructor of the S-VAE where `zdim > 3` and T determines the floating point type (default Float32)
	"""
end
SVAEtwocaps(q, g, hdim::Int, zdim::Int, μ, v::Symbol = :unit, T = Float32) = SVAEtwocaps(q, g, zdim, convert(T, huentropy(zdim)), Adapt.adapt(T, Chain(Dense(hdim, zdim), x -> normalizecolumns(x))), Adapt.adapt(T, Dense(hdim, 1, x -> σ.(x) .* 100)), μ, Flux.param(Adapt.adapt(T, [1.])), Val(v))
SVAEtwocaps(inputDim, hiddenDim, latentDim, numLayers, nonlinearity, layerType, v::Symbol = :unit, T::DataType = Float32) = SVAEtwocaps(inputDim, hiddenDim, latentDim, numLayers, nonlinearity, layerType, Flux.param(Adapt.adapt(T, normalize(randn(latentDim)))), v, T)
function SVAEtwocaps(inputDim, hiddenDim, latentDim, numLayers, nonlinearity, layerType, μ::AbstractVector, v::Symbol = :unit, T = Float32)
	encoder = Adapt.adapt(T, FluxExtensions.layerbuilder(inputDim, hiddenDim, hiddenDim, numLayers, nonlinearity, "", layerType))
	decoder = nothing
	if v == :unit
		decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim, inputDim, numLayers + 1, nonlinearity, "linear", layerType))
	elseif v == :scalarsigma
		decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim, inputDim + 1, numLayers + 1, nonlinearity, "linear", layerType))
	end
	return SVAEtwocaps(encoder, decoder, hiddenDim, latentDim, μ, v, T)
end

Flux.@treelike(SVAEtwocaps)

export closestz, manifoldz, log_pz_from_z, pz_from_z

function closestz(m::SVAEtwocaps, x, steps = 100)
	z = param(Flux.data(zparams(m, x)[1]))
	ps = Flux.Tracker.Params([z])
	opt = ADAM()
	_lkl(model, x, z) = mean(log_pxexpectedz(model, x, z) .- log_pz_from_z(model, z))
	li = Flux.data(_lkl(m, x, z))
	Flux.train!((i) -> -_lkl(m, x, z), ps, 1:steps, opt)
	le = Flux.data(_lkl(m, x, z))
	println("initial = ",li, " final = ",le)
	Flux.data(z)
end

function manifoldz(m::SVAEtwocaps, x, steps = 100)
	z = param(Flux.data(zparams(m, x)[1]))
	ps = Flux.Tracker.Params([z])
	opt = ADAM()
	li = Flux.data(mean(log_pxexpectedz(m, x, z)))
	Flux.train!((i) -> -mean(log_pxexpectedz(m, x, z)), ps, 1:steps, opt)
	le = Flux.data(mean(log_pxexpectedz(m, x, z)))
	println("initial = ",li, " final = ",le)
	Flux.data(z)
end

jacobian_decoder(m::SVAEtwocaps{V}, z) where {V <: Val{:scalarsigma}} = Flux.Tracker.jacobian(a -> m.g(a)[1:end-1], z)

function log_pxexpectedz(m::SVAEtwocaps{V}, x) where {V <: Val{:scalarsigma}}
	xgivenz = m.g(zparams(m, x)[1])[1:end - 1, :]
	log_normal(x, xgivenz, collect(softplus.(xgivenz[end, :])'))
end

function log_pxexpectedz(m::SVAEtwocaps{V}, x, z) where {V <: Val{:scalarsigma}}
	xgivenz = m.g(z)[1:end - 1, :]
	log_normal(x, xgivenz, collect(softplus.(xgivenz[end, :])'))
end

function log_pz(m::SVAEtwocaps, x)
	μz, _ = zparams(m, x)
	return log_vmf_c(Flux.Tracker.data(μz), Flux.Tracker.data(m.priorμ), Flux.Tracker.data(m.priorκ[1]))
end

log_pz_from_z(m::SVAEtwocaps, z) = log_vmf_c(Flux.Tracker.data(z), Flux.Tracker.data(m.priorμ), Flux.Tracker.data(m.priorκ[1]))

pz(m::SVAEtwocaps, x) = exp.(log_pz(m, x))
pz_from_z(m::SVAEtwocaps, z) = exp.(log_pz_from_z(m, z))

function log_px(m::SVAEtwocaps{V}, x::Matrix, k::Int = 100) where {V <: Val{:scalarsigma}}
	x = [x[:, i] for i in 1:size(x, 2)]
	return map(a -> log_px(m, a, k), x)
end

function log_px(m::SVAEtwocaps{V}, x::Vector, k::Int = 100) where {V <: Val{:scalarsigma}}
	μz, κz = zparams(m, x)
	# println("μz: $μz")
	# println("κz: $κz")
	μz = repeat(Flux.Tracker.data(μz), 1, k)
	κz = repeat(Flux.Tracker.data(κz), 1, k)
	z = Flux.Tracker.data(samplez(m, μz, κz))
	xgivenz = Flux.Tracker.data(m.g(z))

	pxgivenz = log_normal(repeat(x, 1, k), xgivenz[1:end - 1, :], collect(softplus.(xgivenz[end, :])'))
	# println("pxgivenz: $pxgivenz")
	pz = log_vmf_wo_c(z, m.priorμ, m.priorκ[1])
	# println("pz $pz")
	qzgivenx = log_vmf_wo_c(z, μz[:, 1], κz[1])
	# println("qzgivenx: $qzgivenx")

	return log(sum(exp.(Flux.Tracker.data(pxgivenz .+ pz .- qzgivenx))))
end

function set_normal_μ(m::SVAEtwocaps, μ)
	κz = 1.
	T = eltype(m.hue)
	m.priorμ = Flux.param(Adapt.adapt(T, Flux.Tracker.data(μ)))
	m.priorκ = Flux.param(Adapt.adapt(T, [κz]))
end

function set_normal_μ_nonparam(m::SVAEtwocaps, μ)
	κz = 1.
	T = eltype(m.hue)
	m.priorμ = Adapt.adapt(T, Flux.Tracker.data(μ))
	m.priorκ = Adapt.adapt(T, [κz])
end

function set_normal_hypersphere(m::SVAEtwocaps, anomaly)
	μz, _ = zparams(m, anomaly)
	set_normal_μ(m, μz)
end

function wloss(m::SVAEtwocaps{V}, x, β, d) where {V <: Val{:unit}}
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	# zp = samplehsuniform(size(z))
	prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
	Ω = d(z, prior)
	xgivenz = m.g(z)
	return Flux.mse(x, xgivenz) + β * Ω
end

function printing_wloss(m::SVAEtwocaps{V}, x, β, d) where {V <: Val{:unit}}
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	# zp = samplehsuniform(size(z))
	prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
	Ω = d(z, prior)
	xgivenz = m.g(z)
	re = Flux.mse(x, xgivenz)
	println("loglkl: $re | Wass-dist: $β x $Ω")
	return re + β * Ω
end

function wloss(m::SVAEtwocaps{V}, x, d) where {V <: Val{:scalarsigma}}
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
	Ω = d(z, prior)
	xgivenz = m.g(z)
	return -mean(log_normal(x, xgivenz[1:end - 1, :], collect(softplus.(xgivenz[end, :])'))) + Ω
end

function printing_wloss(m::SVAEtwocaps{V}, x, d) where {V <: Val{:scalarsigma}}
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
	Ω = d(z, prior)
	xgivenz = m.g(z)
	re = -mean(log_normal(x, xgivenz[1:end - 1, :], collect(softplus.(xgivenz[end, :])')))
	println("loglkl: $re | Wass-dist: $Ω")
	return re + Ω
end

function wloss_semi_supervised(m::SVAEtwocaps{V}, x, y, β, d, α) where {V <: Val{:unit}}
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	xgivenz = m.g(z)

	if count(y .== 1) > 0
		anom_ids = findall(y .== 1)
		anom_ids = anom_ids[rand(1:length(anom_ids), length(y))]
		μzanom = μz[:, anom_ids]
		κzanom = κz[anom_ids]
		zanom = samplez(m, μzanom, collect(κzanom'))

		norm_ids = findall(y .== 0)
		norm_ids = norm_ids[rand(1:length(norm_ids), length(y))]
		μznorm = μz[:, norm_ids]
		κznorm = κz[norm_ids]
		znorm = samplez(m, μznorm, collect(κznorm'))

		anom_prior = samplez(m, ones(size(μz)) .* normalizecolumns(.-m.priorμ), ones(size(κz)) .* m.priorκ)
		norm_prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
		Ωnorm = d(znorm, norm_prior)
		Ωanom = d(zanom, anom_prior)
		return Flux.mse(x, xgivenz) + β * (α .* Ωnorm .+ (1 - α) .* Ωanom)
	else
		norm_prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
		Ωnorm = d(z, norm_prior)
		return Flux.mse(x, xgivenz) + β * Ωnorm
	end
end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
module ClampTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable, @ensureVariables, generateId, addNode!, associate!
using ForneyLab: SPClamp

# Integration helper
mutable struct MockNode <: FactorNode
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function MockNode(out; id=generateId(MockNode))
        @ensureVariables(out)
        self = new(id, Array{Interface}(undef, 1), Dict{Int,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)

        return self
    end
end

@testset "@ensureVariables" begin
    g = FactorGraph()
    nd_c = MockNode(constant(1.0))
    @test isa(nd_c.i[:out].partner.node, Clamp)

    nd_v = MockNode(1.0)
    @test isa(nd_v.i[:out].partner.node, Clamp)
end

@testset "Clamp" begin
    g = FactorGraph()
    nd = Clamp(Variable(), 1.0)

    @test isa(nd, Clamp)
    @test nd.value == 1.0
end

@testset "constant" begin
    g = FactorGraph()
    var = constant(1.0, id=:my_constant)
    nd = g.nodes[:my_constant]

    @test isa(var, Variable)
    @test isa(nd, Clamp)
    @test nd.value == 1.0
end

@testset "placeholder" begin
    g = FactorGraph()

    # Standard placeholder
    var = Variable(id=:y)
    placeholder(var, :y)
    nd = g.nodes[:placeholder_y]

    @test isa(nd, Clamp)
    @test g.placeholders[nd] == (:y, 0)

    # Placeholder without explicit variable initialization
    g = FactorGraph()

    placeholder(:y)
    nd = g.nodes[:placeholder_y]

    @test length(g.variables) == 1
    @test isa(nd, Clamp)
    @test g.placeholders[nd] == (:y, 0)

    # Indexed placeholder
    g = FactorGraph()

    var_i = Variable(id=:y)
    placeholder(var_i, :y, index=1)
    nd_i = g.nodes[:placeholder_y]

    @test isa(nd_i, Clamp)
    @test g.placeholders[nd_i] == (:y, 1)
end

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

@testset "SPClamp" begin
    @test SPClamp{Univariate} <: SumProductRule{Clamp{Univariate}}
    @test outboundType(SPClamp{Univariate}) == Message{PointMass, Univariate}
    @test isApplicable(SPClamp{Univariate}, DataType[])

    @test SPClamp{Multivariate} <: SumProductRule{Clamp{Multivariate}}
    @test outboundType(SPClamp{Multivariate}) == Message{PointMass, Multivariate}
    @test isApplicable(SPClamp{Multivariate}, DataType[])

    @test SPClamp{MatrixVariate} <: SumProductRule{Clamp{MatrixVariate}}
    @test outboundType(SPClamp{MatrixVariate}) == Message{PointMass, MatrixVariate}
    @test isApplicable(SPClamp{MatrixVariate}, DataType[])
end

end #module


text:
# gpu

import .CUDA
import .CUDA: CuArray, CuArrayStyle, @allowscalar
import Adapt: WrappedArray

const CuOneHotArray{K, N, var"N+1"} = OneHotArray{K, N, var"N+1", <: CuArray{OneHot{K}, N}}
const WrappedCuOneHotArray{K, N} = WrappedArray{Bool, N, CuOneHotArray{K, N}, CuOneHotArray{K, N}}

Base.print_array(io::IO, oa::Union{CuOneHotArray, WrappedCuOneHotArray}) = Base.print_array(io, adapt(Array, oa))

Base._show_nonempty(io::IO, oa::Union{CuOneHotArray, WrappedCuOneHotArray}, prefix::String) =
    Base._show_nonempty(io, adapt(Array, oa), prefix)

Base.convert(::Type{T}, oa::Union{CuOneHotArray, WrappedCuOneHotArray}) where {T<:Array} =
    Base.convert(T, adapt(Array, oa))

Base.Array{T, N}(oa::Union{CuOneHotArray, WrappedCuOneHotArray}) where {T, N} = Array{T, N}(adapt(Array, oa))
Base.collect(oa::Union{CuOneHotArray, WrappedCuOneHotArray}) = collect(adapt(Array, oa))

Base.BroadcastStyle(::Type{<:CuOneHotArray{K, N}}) where {K, N} = CuArrayStyle{N+1}()

# avoid array copy
function Base.copyto!(dst::CuArray, oa::OneHotArray)
    dst .= oa
    return dst
end

CUDA.CuArray{T, N, B}(oa::Union{CuOneHotArray, WrappedCuOneHotArray}) where {T, K, N, B, var"N+1", A <: CuArray} =
    copyto!(similar(oa, T), oa)

using Base.Cartesian

function Base.findmax(oa::CuOneHotArray; dims=:)
    if dims == Colon()
        i = findfirst(!iszero, parent(oa))
        isnothing(i) && return (false, CartesianIndex{ndims(oa)}())
        v = @allowscalar oa.onehots[i]
        return (true, CartesianIndex(Int(v), i))
    elseif isone(dims)
        onehots = parent(oa)
        vs = map(!iszero, onehots)
        function f1(x, cidx)
            iszero(x) ? CartesianIndex(1, cidx) : CartesianIndex(Int(x), cidx)
        end
        is = map(f1, onehots, CartesianIndices(onehots))
        return (reshape(vs, 1, size(vs)...), reshape(is, 1, size(is)...))
    else
        return findminmax(>, oa; init=false, dims=dims)
    end
end

function findminmax(binop, a::CuOneHotArray; init, dims)
    function f(t1, t2)
        (x, i), (y, j) = t1, t2

        iszero(i) && return t2
        iszero(j) && return t1
        binop(x, y) && return t1
        x == y && return (x, Base.min(i, j))
        return t2
    end

    indx = ndims(a) == 1 ? (eachindex(a), 1) :
                           (CartesianIndices(a), zero(CartesianIndex{ndims(a)}))

    if dims == Colon()
        mapreduce(tuple, f, a, indx[1]; init = (init, indx[2]))
    else
        res = CUDA.GPUArrays._mapreduce(tuple, f, a, indx[1];
                                        init = (init, indx[2]), dims=dims)
        vals = map(x->x[1], res)
        inds = map(x->x[2], res)
        CUDA.unsafe_free!(res)
        return (vals, inds)
    end
end


text:
using Revise
using Latexify
using DifferentialEquations
using DataFrames
using LaTeXStrings

revise()

ex = :(min(1,2))
latexify(ex)

"x/y" |> latexify
l = latexify("x/y")

join(1, ",")

f = @ode_def feedback begin
    dx = y/c_1 - x
    dy = x^c_2 - y
    end c_1 c_2

latexalign(f.syms, f.funcs)
latexify(f.syms, f.funcs)


text: include("soil.jl")


text:
"""
see https://arxiv.org/abs/1701.07035
"""
@with_kw struct Vumps{F} <: Algorithm
    tol_galerkin::Float64 = Defaults.tol
    tol_gauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
    orthmaxiter::Int = Defaults.maxiter
    finalize::F = Defaults._finalize
    verbose::Bool = Defaults.verbose
end

"
    find_groundstate(state,ham,alg,envs=environments(state,ham))

    find the groundstate for ham using algorithm alg
"

function find_groundstate(state::InfiniteMPS{A,B}, H::Hamiltonian,alg::Vumps,envs::P=environments(state,H))::Tuple{InfiniteMPS{A,B},P,Float64} where {A,B,P}
    galerkin::Float64  = 1+alg.tol_galerkin
    iter      = 1

    temp_ACs = similar(state.AC);
    temp_Cs = similar(state.CR);

    while true
        eigalg = Arnoldi(tol=galerkin/(4*sqrt(iter)))

        @sync for (loc,(ac,c)) in enumerate(zip(state.AC,state.CR))
            @Threads.spawn begin
                (acvals,acvecs) = @closure eigsolve(ac, 1, :SR, eigalg) do x
                    ac_prime(x, loc,state,envs)
                end
                temp_ACs[loc] = acvecs[1];
            end

            @Threads.spawn begin
                (crvals,crvecs) = @closure eigsolve(c, 1, :SR, eigalg) do x
                    c_prime(x, loc,state,envs)
                end

                temp_Cs[loc] = crvecs[1];
            end

        end

        for (i,(ac,c)) in enumerate(zip(temp_ACs,temp_Cs))
            QAc,_ = TensorKit.leftorth!(ac, alg=QRpos())
            Qc,_  = TensorKit.leftorth!(c, alg=QRpos())

            temp_ACs[i] = QAc*adjoint(Qc)
        end

        state = InfiniteMPS(temp_ACs,state.CR[end]; tol = alg.tol_gauge, maxiter = alg.orthmaxiter)
        recalculate!(envs,state);

        galerkin   = calc_galerkin(state, envs)
        alg.verbose && @info "vumps @iteration $(iter) galerkin = $(galerkin)"

        (state,envs) = alg.finalize(iter,state,H,envs) :: Tuple{InfiniteMPS{A,B},P};
        if galerkin <= alg.tol_galerkin || iter>=alg.maxiter
            iter>=alg.maxiter && @warn "vumps didn't converge $(galerkin)"
            return state, envs, galerkin
        end

        iter += 1
    end
end


text:
struct Shift{T}
    value :: T
end

(𝒮::Shift)(a::Sequence) = *(𝒮, a)
Base.:*(𝒮::Shift, a::Sequence) = shift(a, 𝒮.value)

function shift(a::Sequence, τ)
    𝒮 = Shift(τ)
    space_a = space(a)
    new_space = image(𝒮, space_a)
    CoefType = _coeftype(𝒮, space_a, eltype(a))
    c = Sequence(new_space, Vector{CoefType}(undef, dimension(new_space)))
    _apply!(c, 𝒮, a)
    return c
end

function shift!(c::Sequence, a::Sequence, τ)
    𝒮 = Shift(τ)
    space(c) == image(𝒮, space(a)) || return throw(DomainError)
    _apply!(c, 𝒮, a)
    return c
end

function project(𝒮::Shift, domain::VectorSpace, codomain::VectorSpace, ::Type{T}) where {T}
    _iscompatible(domain, codomain) || return throw(DimensionMismatch)
    ind_domain = _findposition_nzind_domain(𝒮, domain, codomain)
    ind_codomain = _findposition_nzind_codomain(𝒮, domain, codomain)
    C = LinearOperator(domain, codomain, sparse(ind_codomain, ind_domain, zeros(T, length(ind_domain)), dimension(codomain), dimension(domain)))
    _project!(C, 𝒮)
    return C
end

function project!(C::LinearOperator, 𝒮::Shift)
    _iscompatible(domain(C), codomain(C)) || return throw(DimensionMismatch)
    coefficients(C) .= zero(eltype(C))
    _project!(C, 𝒮)
    return C
end

_findposition_nzind_domain(𝒮::Shift, domain, codomain) =
    _findposition(_nzind_domain(𝒮, domain, codomain), domain)

_findposition_nzind_codomain(𝒮::Shift, domain, codomain) =
    _findposition(_nzind_codomain(𝒮, domain, codomain), codomain)

# Sequence spaces

image(𝒮::Shift{<:NTuple{N,Any}}, s::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    TensorSpace(map((τᵢ, sᵢ) -> image(Shift(τᵢ), sᵢ), 𝒮.value, spaces(s)))

_coeftype(𝒮::Shift, s::TensorSpace, ::Type{T}) where {T} =
    @inbounds promote_type(_coeftype(Shift(𝒮.value[1]), s[1], T), _coeftype(Shift(Base.tail(𝒮.value)), Base.tail(s), T))
_coeftype(𝒮::Shift, s::TensorSpace{<:Tuple{BaseSpace}}, ::Type{T}) where {T} =
    @inbounds _coeftype(Shift(𝒮.value[1]), s[1], T)

function _apply!(c::Sequence{<:TensorSpace}, 𝒮::Shift, a)
    space_a = space(a)
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    C = _no_alloc_reshape(coefficients(c), dimensions(space(c)))
    _apply!(C, 𝒮, space_a, A)
    return c
end

_apply!(C, 𝒮::Shift, space::TensorSpace{<:NTuple{N₁,BaseSpace}}, A::AbstractArray{T,N₂}) where {N₁,T,N₂} =
    @inbounds _apply!(C, Shift(𝒮.value[1]), space[1], Val(N₂-N₁+1), _apply!(C, Shift(Base.tail(𝒮.value)), Base.tail(space), A))

_apply!(C, 𝒮::Shift, space::TensorSpace{<:Tuple{BaseSpace}}, A::AbstractArray) =
    @inbounds _apply!(C, Shift(𝒮.value[1]), space[1], A)

for (_f, __f) ∈ ((:_nzind_domain, :__nzind_domain), (:_nzind_codomain, :__nzind_codomain))
    @eval begin
        $_f(𝒮::Shift{<:NTuple{N,Any}}, domain::TensorSpace{<:NTuple{N,BaseSpace}}, codomain::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
            Base.Iterators.ProductIterator($__f(𝒮, domain, codomain))
        $__f(𝒮::Shift, domain::TensorSpace, codomain) =
            @inbounds ($_f(Shift(𝒮.value[1]), domain[1], codomain[1]), $__f(Shift(Base.tail(𝒮.value)), Base.tail(domain), Base.tail(codomain))...)
        $__f(𝒮::Shift, domain::TensorSpace{<:Tuple{BaseSpace}}, codomain) =
            @inbounds ($_f(Shift(𝒮.value[1]), domain[1], codomain[1]),)
    end
end

function _project!(C::LinearOperator{<:SequenceSpace,<:SequenceSpace}, 𝒮::Shift)
    domain_C = domain(C)
    codomain_C = codomain(C)
    CoefType = eltype(C)
    @inbounds for (α, β) ∈ zip(_nzind_codomain(𝒮, domain_C, codomain_C), _nzind_domain(𝒮, domain_C, codomain_C))
        C[α,β] = _nzval(𝒮, domain_C, codomain_C, CoefType, α, β)
    end
    return C
end

@generated function _nzval(𝒮::Shift{<:NTuple{N,Any}}, domain::TensorSpace{<:NTuple{N,BaseSpace}}, codomain::TensorSpace{<:NTuple{N,BaseSpace}}, ::Type{T}, α, β) where {N,T}
    p = :(_nzval(Shift(𝒮.value[1]), domain[1], codomain[1], T, α[1], β[1]))
    for i ∈ 2:N
        p = :(_nzval(Shift(𝒮.value[$i]), domain[$i], codomain[$i], T, α[$i], β[$i]) * $p)
    end
    return p
end

# Fourier

image(::Shift, s::Fourier) = s

_coeftype(::Shift{T}, s::Fourier, ::Type{S}) where {T,S} =
    promote_type(typeof(cis(frequency(s)*zero(T))), S)

function _apply!(c::Sequence{<:Fourier}, 𝒮::Shift, a)
    τ = 𝒮.value
    if iszero(τ)
        coefficients(c) .= coefficients(a)
    else
        @inbounds c[0] = a[0]
        eiωτ = cis(frequency(a)*τ)
        eiωτj = one(eiωτ)
        @inbounds for j ∈ 1:order(a)
            eiωτj *= eiωτ
            c[j] = eiωτj * a[j]
            c[-j] = conj(eiωτj) * a[-j]
        end
    end
    return c
end

function _apply!(C, 𝒮::Shift, space::Fourier, ::Val{D}, A) where {D}
    τ = 𝒮.value
    if !iszero(τ)
        ord = order(space)
        eiωτ = cis(frequency(space)*τ)
        eiωτj = one(eiωτ)
        @inbounds for j ∈ 1:ord
            eiωτj *= eiωτ
            selectdim(C, D, ord+1+j) .*= eiωτj
            selectdim(C, D, ord+1-j) .*= conj(eiωτj)
        end
    end
    return C
end

function _apply!(C::AbstractArray{T,N}, 𝒮::Shift, space::Fourier, A) where {T,N}
    τ = 𝒮.value
    if iszero(τ)
        C .= A
    else
        ord = order(space)
        @inbounds selectdim(C, N, 1) .= selectdim(A, N, 1)
        eiωτ = cis(frequency(space)*τ)
        eiωτj = one(eiωτ)
        @inbounds for j ∈ 1:ord
            eiωτj *= eiωτ
            selectdim(C, N, ord+1+j) .= eiωτj .* selectdim(A, N, ord+1+j)
            selectdim(C, N, ord+1-j) .= conj(eiωτj) .* selectdim(A, N, ord+1-j)
        end
    end
    return C
end

function _nzind_domain(::Shift, domain::Fourier, codomain::Fourier)
    ord = min(order(domain), order(codomain))
    return -ord:ord
end
function _nzind_codomain(::Shift, domain::Fourier, codomain::Fourier)
    ord = min(order(domain), order(codomain))
    return -ord:ord
end
function _nzval(𝒮::Shift, domain::Fourier, ::Fourier, ::Type{T}, i, j) where {T}
    τ = 𝒮.value
    if iszero(τ)
        return one(T)
    else
        return convert(T, cis(frequency(domain)*τ*i))
    end
end

# Cartesian spaces

image(𝒮::Shift, s::CartesianPower) =
    CartesianPower(image(𝒮, space(s)), nb_cartesian_product(s))

image(𝒮::Shift, s::CartesianProduct) =
    CartesianProduct(map(sᵢ -> image(𝒮, sᵢ), spaces(s)))

_coeftype(𝒮::Shift, s::CartesianPower, ::Type{T}) where {T} =
    _coeftype(𝒮, space(s), T)

_coeftype(𝒮::Shift, s::CartesianProduct, ::Type{T}) where {T} =
    @inbounds promote_type(_coeftype(𝒮, s[1], T), _coeftype(𝒮, Base.tail(s), T))
_coeftype(𝒮::Shift, s::CartesianProduct{<:Tuple{VectorSpace}}, ::Type{T}) where {T} =
    @inbounds _coeftype(𝒮, s[1], T)

function _apply!(c::Sequence{<:CartesianPower}, 𝒮::Shift, a)
    @inbounds for i ∈ 1:nb_cartesian_product(space(c))
        _apply!(component(c, i), 𝒮, component(a, i))
    end
    return c
end
function _apply!(c::Sequence{CartesianProduct{T}}, 𝒮::Shift, a) where {N,T<:NTuple{N,VectorSpace}}
    @inbounds _apply!(component(c, 1), 𝒮, component(a, 1))
    @inbounds _apply!(component(c, 2:N), 𝒮, component(a, 2:N))
    return c
end
function _apply!(c::Sequence{CartesianProduct{T}}, 𝒮::Shift, a) where {T<:Tuple{VectorSpace}}
    @inbounds _apply!(component(c, 1), 𝒮, component(a, 1))
    return c
end

function _findposition_nzind_domain(𝒮::Shift, domain::CartesianSpace, codomain::CartesianSpace)
    u = map((dom, codom) -> _findposition_nzind_domain(𝒮, dom, codom), spaces(domain), spaces(codomain))
    len = sum(length, u)
    v = Vector{Int}(undef, len)
    δ = δδ = 0
    @inbounds for (i, uᵢ) in enumerate(u)
        δ_ = δ
        δ += length(uᵢ)
        view(v, 1+δ_:δ) .= δδ .+ uᵢ
        δδ += dimension(domain[i])
    end
    return v
end

function _findposition_nzind_codomain(𝒮::Shift, domain::CartesianSpace, codomain::CartesianSpace)
    u = map((dom, codom) -> _findposition_nzind_codomain(𝒮, dom, codom), spaces(domain), spaces(codomain))
    len = sum(length, u)
    v = Vector{Int}(undef, len)
    δ = δδ = 0
    @inbounds for (i, uᵢ) in enumerate(u)
        δ_ = δ
        δ += length(uᵢ)
        view(v, 1+δ_:δ) .= δδ .+ uᵢ
        δδ += dimension(codomain[i])
    end
    return v
end

function _project!(C::LinearOperator{<:CartesianSpace,<:CartesianSpace}, 𝒮::Shift)
    @inbounds for i ∈ 1:nb_cartesian_product(domain(C))
        _project!(component(C, i, i), 𝒮)
    end
    return C
end


text:
using Distributed
@everywhere using Mamba, LinearAlgebra

## Data
birats = Dict{Symbol, Any}(
  :N => 30, :T => 5,
  :x => [8.0, 15.0, 22.0, 29.0, 36.0],
  :Y => [151 199 246 283 320
         145 199 249 293 354
         147 214 263 312 328
         155 200 237 272 297
         135 188 230 280 323
         159 210 252 298 331
         141 189 231 275 305
         159 201 248 297 338
         177 236 285 350 376
         134 182 220 260 296
         160 208 261 313 352
         143 188 220 273 314
         154 200 244 289 325
         171 221 270 326 358
         163 216 242 281 312
         160 207 248 288 324
         142 187 234 280 316
         156 203 243 283 317
         157 212 259 307 336
         152 203 246 286 321
         154 205 253 298 334
         139 190 225 267 302
         146 191 229 272 302
         157 211 250 285 323
         132 185 237 286 331
         160 207 257 303 345
         169 216 261 295 333
         157 205 248 289 316
         137 180 219 258 291
         153 200 244 286 324],
  :mean => [0.0, 0.0],
  :var => [1.0e6 0.0
           0.0   1.0e6],
  :Omega => [200.0 0.0
               0.0 0.2]
)

## Model Specification
model = Model(

  Y = Stochastic(2,
    (beta, x, sigmaC, N, T) ->
      UnivariateDistribution[
        Normal(beta[i, 1] + beta[i, 2] * x[j], sigmaC)
        for i in 1:N, j in 1:T
      ],
    false
  ),

  beta = Stochastic(2,
    (mu_beta, Sigma, N) ->
      MultivariateDistribution[
        MvNormal(mu_beta, Sigma)
        for i in 1:N
      ],
    false
  ),

  mu_beta = Stochastic(1,
    (mean, var) -> MvNormal(mean, var)
  ),

  Sigma = Stochastic(2,
    Omega -> InverseWishart(2, Omega),
    false
  ),

  sigma2C = Stochastic(
    () -> InverseGamma(0.001, 0.001),
    false
  ),

  sigmaC = Logical(
    sigma2C -> sqrt(sigma2C)
  )

)

## Initial Values
inits = [
  Dict(:Y => birats[:Y], :beta => repeat([100 6], birats[:N], 1),
       :mu_beta => [0, 0], :Sigma => Matrix{Float64}(I, 2, 2), :sigma2C => 1.0),
  Dict(:Y => birats[:Y], :beta => repeat([50 3], birats[:N], 1),
       :mu_beta => [10, 10], :Sigma => 0.3 * Matrix{Float64}(I, 2, 2), :sigma2C => 10.0)
]

## Sampling Scheme
scheme = [AMWG([:beta, :mu_beta], repeat([10.0, 1.0], birats[:N] + 1)),
          AMWG(:Sigma, 1.0),
          Slice(:sigma2C, 10.0)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, birats, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
using NoiseReduction
using Test

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


text:
module Gates
"""
Minimal definition of circuit components provided for testing convenience.
QXZoo contains more comprehensive gate set and features
"""

h() = [[1., 1.] [1., -1.]]./sqrt(2)
x() = [[0., 1.] [1., 0.]]
z() = [[1., 0.] [0., -1.]]
I() = [[1., 0.] [0., 1.]]
cx() = [[1., 0., 0., 0.] [0., 1., 0., 0.] [0., 0., 0., 1.] [0., 0., 1., 0.]]
cz() = [[1., 0., 0., 0.] [0., 1., 0., 0.] [0., 0., 1., 0.] [0., 0., 0., -1.]]
swap() = [[1., 0., 0., 0.] [0., 0., 1., 0.] [0., 1., 0., 0.] [0., 0., 0., 1.]]

end


text:
module TwitterUploadMedia

using OAuth
using JSON
using Twitter: twitterauth,
        get_oauth,
        post_oauth,
        DataFrame,
        TwitterType

export  get_media_upload,
        post_media_upload,
        post_media_metadata_create,
        post_media_subtitles_delete,
        post_media_subtitles_create

include("mediamethods.jl")
include("types.jl")

end # module


text:
function get_stage()
    # Appveyor uses "True" for non-Ubuntu images.
    if get(ENV, "APPVEYOR", "") == "True" || get(ENV, "APPVEYOR", "") == "true"
        return "nonnumeric"
    end

    # Handle Travis and Github Actions specially.
    if get(ENV, "TRAVIS", "") == "true" || get(ENV, "GITHUB_ACTIONS", "") == "true"
        if "STAGE" in keys(ENV)
            return ENV["STAGE"]
        else
            return "all"
        end
    end

    return "all"
end

function do_test(stage_str)
    stg = get_stage()

    # If the tests are being run by Appveyor, don't run
    # any numerical tests.
    if stg == "nonnumeric"
        if stage_str == "numerical"
            return false
        else
            return true
        end
    end

    # Otherwise run the regular testing procedure.
    if stg == "all" || stg == stage_str
        return true
    end

    return false
end

macro stage_testset(stage_string::String, args...)
    if do_test(stage_string)
        return esc(:(@testset($(args...))))
    end
end

macro numerical_testset(args...)
    esc(:(@stage_testset "numerical" $(args...)))
end

macro turing_testset(args...)
    esc(:(@stage_testset "test" $(args...)))
end


text:
### Override for Graphs' version, copied from https://github.com/JuliaGraphs/Graphs.jl/blob/master/src/shortestpaths/dijkstra.jl with some adjustments
function dijkstra_shortest_paths(g::AbstractIndexedGraph,
    srcs::Vector{U},
    distvec::AbstractVector{T}=ones(Int, ne(g));
    allpaths=false,
    trackvertices=false
    ) where T <: Real where U <: Integer

    issubset(1:ne(g), eachindex(distvec)) || throw(ArgumentError("Bad indices in distvec"))
    nvg = nv(g)
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)

    pathcounts = zeros(nvg)
    preds = fill(Vector{U}(), nvg)
    H::TrackingHeap{Int64, T, 2, MinHeapOrder, NoTrainingWheels} = TrackingHeap(T; S=NoTrainingWheels)
    sizehint!(H, nvg)

    for src in srcs
        dists[src] = zero(T)
        pathcounts[src] = one(Float64)
        H[src] = zero(T)
    end

    closest_vertices = Vector{U}()  # Maintains vertices in order of distances from source

    trackvertices && sizehint!(closest_vertices, nvg)

    @inbounds while !isempty(H)
        u, d = pop!(H)

        if trackvertices
            push!(closest_vertices, u)
        end

        for e in outedges(g, u)
            alt = d + distvec[idx(e)]
            v = dst(e)

            visited = dists[v] != typemax(T)
            if !visited
                dists[v] = alt
                parents[v] = u

                pathcounts[v] += pathcounts[u]
                if allpaths
                    preds[v] = [u;]
                end
                H[v] = alt
            elseif alt < dists[v]
                dists[v] = alt
                parents[v] = u
                #615
                pathcounts[v] = pathcounts[u]
                if allpaths
                    resize!(preds[v], 1)
                    preds[v][1] = u
                end
                H[v] = alt
            elseif alt == dists[v]
                pathcounts[v] += pathcounts[u]
                if allpaths
                    push!(preds[v], u)
                end
            end
        end
    end

    if trackvertices
        for s in vertices(g)
            visited = dists[s] != typemax(T)
            if !visited
                push!(closest_vertices, s)
            end
        end
    end

    for src in srcs
        pathcounts[src] = one(Float64)
        parents[src] = 0
        empty!(preds[src])
    end

    return DijkstraState{T,U}(parents, dists, preds, pathcounts, closest_vertices)
end

function dijkstra_shortest_paths(g::AbstractIndexedGraph, src::Integer,
        distvec::AbstractVector=ones(Int, ne(g)); kw...)
    dijkstra_shortest_paths(g, [src;], distvec; kw...)
end

function dijkstra_shortest_paths(g::AbstractIndexedGraph, src,
        distmx::AbstractMatrix; kw...)
    Wt = convert(SparseMatrixCSC, transpose(distmx))
    distvec = nonzeros(Wt)
    dijkstra_shortest_paths(g, src, distvec; kw...)
end

### Homemade minimal version
function dijkstra(g::AbstractIndexedGraph{T}, s::Integer,
        w::AbstractVector{U}=ones(Int, ne(g))) where {T<:Integer, U<:Real}
    n = nv(g)
    inf = typemax(U)
    D = fill(inf, n); D[s] = 0
    parents = fill(0, n); parents[s] = 0
    Q::TrackingHeap{Int64, U, 2, MinHeapOrder, NoTrainingWheels} = TrackingHeaps.TrackingHeap(U; S=TrackingHeaps.NoTrainingWheels)
    sizehint!(Q, n)
    Q[s] = zero(U)
    while !isempty(Q)
        v, d = TrackingHeaps.pop!(Q)
        for e in outedges(g, v)
            i = dst(e)
            dd = d + w[idx(e)]
            if dd < D[i]
                Q[i] = D[i] = dd
                parents[i] = v
            end
        end
    end
    return D, parents
end


text:
module parent

using FromFile
@from "revise/subfile.jl" import child

g1() = child.f1()^2
g2() = child.f2()^2
g3() = child.f3()^2

end


text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

mutable struct SubDatabase{D<:AbstractDict} <: AbstractDatabase{D}
    parent::AbstractDatabase
    key::Symbol
    datadict::D
    function SubDatabase{D}(
        parent::AbstractDatabase,
        key::Union{Symbol,AbstractString},
        data::D = datadict(parent)[key],
    ) where {D<:AbstractDict}
        new(parent, Symbol(key), data)
    end
end

function SubDatabase{D}(
    parent::SubDatabase{D},
    datakeys::Union{Tuple,AbstractVector},
) where {D<:AbstractDict}
    for key ∈ Symbol.(datakeys)
        data = datadict(parent)
        parent = SubDatabase(parent, key, data[key])
    end
    parent
end

Base.parent(database::SubDatabase) = database.parent

sources(database::SubDatabase) = sources(parent(database))
sources!(database::SubDatabase, value) = sources!(parent(database), value)

default_kwargs(database::SubDatabase) = default_kwargs(parent(database))
default_kwargs!(database::SubDatabase, value) = default_kwargs!(parent(database), value)

computing_functions(database::SubDatabase) = computing_functions(parent(database))
computing_functions!(database::SubDatabase, value) =
    computing_functions!(parent(database), value)

modeldict(database::SubDatabase) = modeldict(parent(database))[database.key]
modeldict!(database::SubDatabase, value) = modeldict(parent(database))[database.key] = value

root(database::SubDatabase) = root(parent(database))

function _find_subdata(data::AbstractDict, key)
    result = get(data, key, missing)
    !ismissing(result) && return result

    results = skipmissing(_find_subdata(subdata, key) for subdata ∈ values(data))
    results = collect(results)

    length(results) == 1 && return first(results)
    length(results) > 1 &&
        throw(ErrorException("More than one value found for key $(key)."))
    throw(ErrorException("No value found for key $(key)."))
end

_find_subdata(data::AbstractDatabase, key) = _find_subdata(datadict(data), key)
_find_subdata(data, key) = missing

function Base.getproperty(
    database::AbstractDatabase{D},
    key::Symbol,
) where {D<:AbstractDict}
    if hasfield(typeof(database), key)
        getfield(database, key)
    elseif key == :model
        SubDatabase{D}(database, :model, modeldict(database))
    else
        data = _find_subdata(database, key)
        data = data isa AbstractDict ? SubDatabase{D}(database, key) : data
    end
end

function Base.propertynames(database::AbstractDatabase{D}) where {D<:AbstractDict}
    props = copy(collect(fieldnames(typeof(database))))
    data = datadict(database)
    append!(props, keys(data))
end

function Base.getindex(
    database::AbstractDatabase{D},
    key::Union{Symbol,AbstractString},
    keys::Union{Symbol,AbstractString}...,
) where {D<:AbstractDict}
    database = getproperty(database, Symbol(key))
    for key ∈ keys
        database = getproperty(database, Symbol(key))
    end
    database
end

function Base.setindex!(
    database::AbstractDatabase{D},
    value,
    key::Union{Symbol,AbstractString},
    keys::Union{Symbol,AbstractString}...,
) where {D<:AbstractDict}
    key = Symbol(key)
    if isempty(keys)
        datadict(database)[key] = value
    else
        data, key = get!(datadict(database), key, D()), Symbol(pop!(keys))
        while !isempty(keys)
            data, key = get!(data, key, D()), Symbol(pop!(keys))
        end
        data[key] = value
    end
end

function Base.get(database::AbstractDatabase, key::Union{Symbol,AbstractString}, default)
    get(datadict(database), Symbol(key), default)
end

function Base.get(
    f::Base.Callable,
    database::AbstractDatabase,
    key::Union{Symbol,AbstractString},
)
    get(f, datadict(database), Symbol(key))
end

function Base.get!(database::AbstractDatabase, key::Union{Symbol,AbstractString}, default)
    get!(datadict(database), Symbol(key), default)
end


text:
function hexcolor(rgb::UInt32)
    r = ((0xff0000 & rgb) >> 16) / 255
    g = ((0x00ff00 & rgb) >>  8) / 255
    b = ((0x0000ff & rgb)      ) / 255
    Plots.RGBA(r, g, b)
end

# Values taken from http://www.toutes-les-couleurs.com/code-couleur-rvb.php
troyes = hexcolor(0xfefdf0)
#troyes = Plots.RGBA(0.9961, 0.9922, 0.9412)
frambo = hexcolor(0xc72c48)
#frambo = Plots.RGBA(0.7804, 0.1725, 0.2824)
lichen = hexcolor(0x85c17e)
#lichen = Plots.RGBA(0.5216, 0.7569, 0.4941)
canard = hexcolor(0x048b9a)
#canard = Plots.RGBA(0.0157, 0.5451, 0.6039)
aurore = hexcolor(0xffcb60)
# colors taken from https://previews.123rf.com/images/capacitorphoto/capacitorphoto1410/capacitorphoto141000191/32438941-Graph-Icon-color-set-illustration--Stock-Vector.jpg
re = hexcolor(0xf59297)
gre = hexcolor(0xcbdf80)
blu = hexcolor(0x2394ce)
ora = hexcolor(0xfacb95)
yel = hexcolor(0xf2f08b)


text:
"""
    flatten(vecofvec) -> vec

Flatten a nested vector, returning a single vector containing the elements
of each nested vector.

This is used for **@piping.**

**Warning.**
The type of the returned vector should be checked before further use.

**Example**
```julia
flatten([[1, 2], [3, 4]])  # returns [1, 2, 3, 4]
```
"""
flatten(vecofvec) = vcat(vecofvec...)


text:
export
ruleSPNonlinearSOutNM,
ruleSPNonlinearSIn1MN,
ruleSPNonlinearSOutNGX,
ruleSPNonlinearSInGX,
ruleSPNonlinearSOutNMX,
ruleSPNonlinearSInMX,
ruleMNonlinearSInMGX,
prod!

const default_n_samples = 1000 # Default value for the number of samples

#----------------------
# Sampling Update Rules
#----------------------

function ruleSPNonlinearSOutNM(g::Function,
                               msg_out::Nothing,
                               msg_in1::Message; # Applies to any message except SampleList
                               dims::Any=nothing,
                               n_samples=default_n_samples)

    samples = g.(sample(msg_in1.dist, n_samples))
    weights = ones(n_samples)/n_samples

    return Message(variateType(dims), SampleList, s=samples, w=weights)
end

function ruleSPNonlinearSOutNM(g::Function,
                               msg_out::Nothing,
                               msg_in1::Message{SampleList}; # Special case for SampleList
                               dims::Any=nothing,
                               n_samples=default_n_samples)

    samples = g.(msg_in1.dist.params[:s])
    weights = msg_in1.dist.params[:w]

    return Message(variateType(dims), SampleList, s=samples, w=weights)
end

function ruleSPNonlinearSIn1MN(g::Function,
                               msg_out::Message,
                               msg_in1::Nothing;
                               dims::Any=nothing,
                               n_samples=default_n_samples)

    return Message(variateType(dims), Function, log_pdf = (z)->logPdf(msg_out.dist, g(z)))
end

function ruleSPNonlinearSOutNGX(g::Function,
                                msg_out::Nothing,
                                msgs_in::Vararg{Message{<:Gaussian}};
                                dims::Any=nothing,
                                n_samples=default_n_samples)

    samples_in = [sample(msg_in.dist, n_samples) for msg_in in msgs_in]

    samples = g.(samples_in...)
    weights = ones(n_samples)/n_samples

    return Message(variateType(dims), SampleList, s=samples, w=weights)
end

function ruleSPNonlinearSInGX(g::Function,
                              inx::Int64, # Index of inbound interface inx
                              msg_out::Message,
                              msgs_in::Vararg{Message{<:Gaussian}};
                              dims::Any=nothing,
                              n_samples=default_n_samples)

    # Extract joint statistics of inbound messages
    (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Return arrays with individual means and covariances
    (m_fw_in, V_fw_in, ds) = concatenateGaussianMV(ms_fw_in, Vs_fw_in) # Concatenate individual statistics into joint statistics
    W_fw_in = cholinv(V_fw_in) # Convert to canonical statistics

    # Construct joint log-pdf function and gradient
    (log_joint, d_log_joint) = logJointPdfs(m_fw_in, W_fw_in, msg_out.dist, g, ds)

    # Compute joint belief on in's by gradient ascent
    m_in = gradientOptimization(log_joint, d_log_joint, m_fw_in, 0.01)
    V_in = cholinv(-ForwardDiff.jacobian(d_log_joint, m_in))

    # Marginalize joint belief on in's
    (m_inx, V_inx) = marginalizeGaussianMV(m_in, V_in, ds, inx)
    W_inx = cholinv(V_inx) # Convert to canonical statistics
    xi_inx = W_inx*m_inx

    # Divide marginal on inx by forward message
    (xi_fw_inx, W_fw_inx) = unsafeWeightedMeanPrecision(msgs_in[inx].dist)
    xi_bw_inx = xi_inx - xi_fw_inx
    W_bw_inx = W_inx - W_fw_inx # Note: subtraction might lead to posdef violations

    return Message(variateType(dims), GaussianWeightedMeanPrecision, xi=xi_bw_inx, w=W_bw_inx)
end

function ruleSPNonlinearSOutNMX(g::Function,
                                msg_out::Nothing,
                                msgs_in::Vararg{Message};
                                dims::Any=nothing,
                                n_samples=default_n_samples)

    samples_in = [sample(msg_in.dist, n_samples) for msg_in in msgs_in]
    samples = g.(samples_in...)
    weights = ones(n_samples)/n_samples

    return Message(variateType(dims), SampleList, s=samples, w=weights)
end

function ruleSPNonlinearSInMX(g::Function,
                              inx::Int64, # Index of inbound interface inx
                              msg_out::Message,
                              msgs_in::Vararg{Message};
                              dims::Any=nothing,
                              n_samples=default_n_samples)
    arg_sample = (z) -> begin
        samples_in = []
        for i=1:length(msgs_in)
            if i==inx
                push!(samples_in, collect(Iterators.repeat([z], n_samples)))
            else
                push!(samples_in, sample(msgs_in[i].dist, n_samples))
            end
        end

        return samples_in
    end

    approximate_pdf(z) = sum(exp.(logPdf.([msg_out.dist],g.(arg_sample(z)...))))/n_samples

    return Message(variateType(dims), Function, log_pdf = (z)->log(approximate_pdf(z)))
end

# Special case for two inputs with one PointMass (no inx required)
function ruleSPNonlinearSInMX(g::Function,
                              msg_out::Message,
                              msg_in1::Message{PointMass},
                              msg_in2::Nothing;
                              dims::Any=nothing,
                              n_samples=default_n_samples)

    m_in1 = msg_in1.dist.params[:m]

    return Message(variateType(dims), Function, log_pdf = (z)->logPdf(msg_out.dist, g(m_in1, z)))
end

# Special case for two inputs with one PointMass (no inx required)
function ruleSPNonlinearSInMX(g::Function,
                              msg_out::Message,
                              msg_in1::Nothing,
                              msg_in2::Message{PointMass};
                              dims::Any=nothing,
                              n_samples=default_n_samples)

    m_in2 = msg_in2.dist.params[:m]

    return Message(variateType(dims), Function, log_pdf = (z)->logPdf(msg_out.dist, g(z, m_in2)))
end

function ruleMNonlinearSInMGX(g::Function,
                              msg_out::Message,
                              msgs_in::Vararg{Message{<:Gaussian}})

    # Extract joint statistics of inbound messages
    (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Return arrays with individual means and covariances
    (m_fw_in, V_fw_in, ds) = concatenateGaussianMV(ms_fw_in, Vs_fw_in) # Concatenate individual statistics into joint statistics
    W_fw_in = cholinv(V_fw_in) # Convert to canonical statistics

    # Construct log-pdf function and gradient
    (log_joint, d_log_joint) = logJointPdfs(m_fw_in, W_fw_in, msg_out.dist, g, ds)

    # Compute joint marginal belief on in's by gradient ascent
    m_in = gradientOptimization(log_joint, d_log_joint, m_fw_in, 0.01)
    W_in = -ForwardDiff.jacobian(d_log_joint, m_in)

    return ProbabilityDistribution(Multivariate, GaussianMeanPrecision, m=m_in, w=W_in)
end

#---------------------------
# Custom inbounds collectors
#---------------------------

# Sampling approximation
function collectSumProductNodeInbounds(node::Nonlinear{Sampling}, entry::ScheduleEntry)
    inbounds = Any[]

    # Push function to calling signature
    # This function needs to be defined in the scope of the user
    push!(inbounds, Dict{Symbol, Any}(:g => node.g,
                                      :keyword => false))

    multi_in = isMultiIn(node) # Boolean to indicate a nonlinear node with multiple stochastic inbounds
    inx = findfirst(isequal(entry.interface), node.interfaces) - 1 # Find number of inbound interface; 0 for outbound

    if (inx > 0) && multi_in # Multi-inbound backward rule
        push!(inbounds, Dict{Symbol, Any}(:inx => inx, # Push inbound identifier
                                          :keyword => false))
    end

    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry
    for node_interface in node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if (node_interface == entry.interface != node.interfaces[1]) && multi_in
            # Collect the breaker message for a backward rule with multiple inbounds
            haskey(interface_to_schedule_entry, inbound_interface) || error("The nonlinear node's backward rule uses the incoming message on the input edge to determine the approximation point. Try altering the variable order in the scheduler to first perform a forward pass.")
            push!(inbounds, interface_to_schedule_entry[inbound_interface])
        elseif node_interface == entry.interface
            # Ignore inbound message on outbound interface
            push!(inbounds, nothing)
        elseif isClamped(inbound_interface)
            # Hard-code outbound message of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, Message))
        else
            # Collect message from previous result
            push!(inbounds, interface_to_schedule_entry[inbound_interface])
        end
    end

    # Push custom arguments if defined
    if (node.dims !== nothing)
        push!(inbounds, Dict{Symbol, Any}(:dims => node.dims[inx + 1],
                                          :keyword => true))
    end
    if (node.n_samples !== nothing)
        push!(inbounds, Dict{Symbol, Any}(:n_samples => node.n_samples,
                                          :keyword   => true))
    end
    return inbounds
end

#---------------------------
# Custom product definitions
#---------------------------

function prod!(
    x::ProbabilityDistribution{V, Function},
    y::ProbabilityDistribution{V, Function}) where V<:VariateType # log-pdf for z cannot be predefined, because it cannot be overwritten

    return ProbabilityDistribution(V, Function, log_pdf=(s)->x.params[:log_pdf](s) + y.params[:log_pdf](s))
end

@symmetrical function prod!(
    x::ProbabilityDistribution{Univariate}, # Includes function distributions
    y::ProbabilityDistribution{Univariate, <:Gaussian},
    z::ProbabilityDistribution{Univariate, GaussianMeanPrecision}=ProbabilityDistribution(Univariate, GaussianMeanPrecision, m=0.0, w=1.0))

    # Optimize with gradient ascent
    log_joint(s) = logPdf(y,s) + logPdf(x,s)
    d_log_joint(s) = ForwardDiff.derivative(log_joint, s)
    m_initial = unsafeMean(y)

    mean = gradientOptimization(log_joint, d_log_joint, m_initial, 0.01)
    prec = -ForwardDiff.derivative(d_log_joint, mean)

    z.params[:m] = mean
    z.params[:w] = prec

    return z
end

@symmetrical function prod!(
    x::ProbabilityDistribution{Multivariate}, # Includes function distributions
    y::ProbabilityDistribution{Multivariate, <:Gaussian},
    z::ProbabilityDistribution{Multivariate, GaussianMeanPrecision}=ProbabilityDistribution(Multivariate, GaussianMeanPrecision, m=[0.0], w=mat(1.0)))

    # Optimize with gradient ascent
    log_joint(s) = logPdf(y,s) + logPdf(x,s)
    d_log_joint(s) = ForwardDiff.gradient(log_joint, s)
    m_initial = unsafeMean(y)

    mean = gradientOptimization(log_joint, d_log_joint, m_initial, 0.01)
    prec = -ForwardDiff.jacobian(d_log_joint, mean)

    z.params[:m] = mean
    z.params[:w] = prec

    return z
end

#---------------------------------
# Gradient optimization subroutine
#---------------------------------

function gradientOptimization(log_joint::Function, d_log_joint::Function, m_initial, step_size)
    dim_tot = length(m_initial)
    m_total = zeros(dim_tot)
    m_average = zeros(dim_tot)
    m_new = zeros(dim_tot)
    m_old = m_initial
    satisfied = false
    step_count = 0
    m_latests = if (dim_tot == 1) Queue{Float64}() else Queue{Vector}() end

    while !satisfied
        m_new = m_old .+ step_size.*d_log_joint(m_old)
        if log_joint(m_new) > log_joint(m_old)
            proposal_step_size = 10*step_size
            m_proposal = m_old .+ proposal_step_size.*d_log_joint(m_old)
            if log_joint(m_proposal) > log_joint(m_new)
                m_new = m_proposal
                step_size = proposal_step_size
            end
        else
            step_size = 0.1*step_size
            m_new = m_old .+ step_size.*d_log_joint(m_old)
        end
        step_count += 1
        enqueue!(m_latests, m_old)
        if step_count > 10
            m_average = sum(x for x in m_latests)./10
            if sum(sqrt.(((m_new.-m_average)./m_average).^2)) < dim_tot*0.1
                satisfied = true
            end
            dequeue!(m_latests);
        end
        if step_count > dim_tot*250
            satisfied = true
        end
        m_old = m_new
    end

    return m_new
end

#--------
# Helpers
#--------

function logJointPdfs(m_fw_in::Vector, W_fw_in::AbstractMatrix, dist_out::ProbabilityDistribution, g::Function, ds::Vector)
    log_joint(x) = -0.5*sum(intdim.(ds))*log(2pi) + 0.5*logdet(W_fw_in) - 0.5*(x - m_fw_in)'*W_fw_in*(x - m_fw_in) + logPdf(dist_out, g(split(x, ds)...))
    d_log_joint(x) = ForwardDiff.gradient(log_joint, x)

    return (log_joint, d_log_joint)
end


text:
using Test
using WignerSymbols
using LinearAlgebra
using Random

Random.seed!(1234)
smalljlist = 0:1//2:10
largejlist = 0:1//2:1000

cache = WignerCache()

# test recurrence relations: Phys Rev E 57, 7274 (1998)
@testset "wigner3j: test recurrence relations" begin
    for k = 1:10
        j2 = convert(BigFloat, rand(0:1//2:1000))
        j3 = convert(BigFloat, rand(0:1//2:1000))
        m2 = convert(BigFloat, rand(-j2:j2))
        m3 = convert(BigFloat, rand(-j3:j3))

        for j in max(abs(j2-j3),abs(m2+m3))+1:(j2+j3)-1
            X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(m2+m3)^2))
            Y = (2*j+1)*((m2+m3)*(j2*(j2+1)-j3*(j3+1)) - (m2-m3)*j*(j+1))
            Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(m2+m3)^2))
            tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat)
            @test (X*wigner3j(cache, BigFloat,j+1,j2,j3,-m2-m3,m2,m3) +
                        Z*wigner3j(cache, BigFloat,j-1,j2,j3,-m2-m3,m2,m3)) ≈
                            (-Y*wigner3j(cache, BigFloat,j,j2,j3,-m2-m3,m2,m3)) atol=tol
        end
    end
end

@testset "wigner3j: test orthogonality relations" begin
    # equivalent to Clebsch-Gordan orthogonality, now test using Float32
    for j1 in smalljlist, j2 in smalljlist
        d1::Int = 2*j1+1
        d2::Int = 2*j2+1
        M = zeros(Float32, (d1*d2, d1*d2))
        ind2 = 1
        for m1 in -j1:j1, m2 in -j2:j2
            ind1 = 1
            @inbounds for j3 in abs(j1-j2):(j1+j2), m3 in -j3:j3
                d3::Int = 2*j3+1
                M[ind1,ind2] += sqrt(d3) * wigner3j(cache, Float32, j1, j2, j3, m1, m2, m3)
                ind1 += 1
            end
            ind2 += 1
        end
        @test M'*M ≈ one(M) # orthogonality relation type 1
        @test M*M' ≈ one(M) # orthogonality relation type 2
    end
end