text:
"""
   vbmf_sparse_parameters

Compound type for vbmf_sparse computation. Contains the following fields:\n

    L::Int - number of rows of the original matrix
    M::Int - number of columns of the original matrix
    H::Int - internal product dimension
    MH::Int - dimension of vec(A')
    H1::Int - number of columns of B that should belong to a contaminated class
    labels::Array{Int64, 1} - which columns of Y are labeled as non-contaminated

    AHat::Array{Float64, 2} - mean value of A, size (M, H)
    ATVecHat::Array{Float64, 1} - mean value of vec(A^T), size (MH,1)
    SigmaATVec::Array{Float64, 2} - covariance of vec(A^T), size (MH, MH)
    diagSigmaATVec::Array{Float64,1} - diagonal of covariance of vec(A^T), size (MH)
    invSigmaATVec::Array{Float64, 2} - inverse of covariance of vec(A^T), size (MH, MH)
    SigmaA::Array{Float64,2} - covariance of A matrix, size (H, H)

    BHat::Array{Float64, 2} - mean value of B, size (L, H)
    SigmaB::Array{Float64, 2} - covariance of B matrix, size (H, H)

    CA::Array{Float64, 1} - diagonal of the prior covariance of vec(A^T), size (MH)
    alpha0::Float64 - shape of CA gamma prior
    beta0::Float64 - scale of CA gamma prior
    alpha::Float64 - shape of CA gamma posterior
    beta::Array{Float64,1} - scale of CA gamma posterior, size (MH, 1)

    CB::Array{Float64, 1} - diagonal of prior covariance of B, size (H)
    gamma0::Float64 - shape of CB gamma prior
    delta0::Float64 - scale of CB gamma prior
    gamma::Float64 - shape of CB posterior
    delta::Array{Float64,1} - scale of CB posterior, size (H)

    sigmaHat::Float64 - variance of data - homoscedastic case
    eta0::Float64 - shape of sigma gamma prior
    zeta0::Float64 - scale of sigma gamma prior
    eta::Float64 - shape of homoscedastic sigma posterior
    zeta::Float64 - scale of homoscedastic sigma posterior
    sigmaVecHat::Array{Float64,1} - variance of data (rows) - heteroscedastic case, size (L)
    etaVec::Array{Float64,1} - shape of heteroscedastic sigma posterior, size (L)
    zetaVec::Array{Float64,1} - scale of heteroscedastic sigma posterior, size (L)

    YHat::Array{Float64, 2} - estimate of Y, size (L, M)
    trYTY::Float64 - trace(Y^T*Y), saved so that it does not have to be recomputed
"""
type vbmf_sparse_parameters
    L::Int
    M::Int
    H::Int
    MH::Int
    H1::Int
    labels::Array{Int64,1}

    AHat::Array{Float64, 2}
    ATVecHat::Array{Float64, 1}
    SigmaATVec::Array{Float64, 2}
    diagSigmaATVec::Array{Float64,1}
    invSigmaATVec::Array{Float64, 2}
    SigmaA::Array{Float64,2}

    BHat::Array{Float64, 2}
    SigmaB::Array{Float64, 2}

    CA::Array{Float64, 1}
    alpha0::Float64
    beta0::Float64
    alpha::Float64
    beta::Array{Float64,1}

    CB::Array{Float64, 1}
    gamma0::Float64
    delta0::Float64
    gamma::Float64
    delta::Array{Float64,1}

    sigmaHat::Float64
    eta0::Float64
    zeta0::Float64
    eta::Float64
    zeta::Float64
    sigmaVecHat::Array{Float64,1}
    etaVec::Array{Float64,1}
    zetaVec::Array{Float64,1}

    YHat::Array{Float64, 2}
    trYTY::Float64

    vbmf_sparse_parameters() = new()
end

"""
    vbmf_sparse_init(Y::Array{Float64,2}, H::Int; ca::Float64 = 1.0,
    alpha0::Float64 = 1e-10, beta0::Float64 = 1e-10, cb::Float64 = 1.0,
    gamma0::Float64 = 1e-10, delta0::Float64 = 1e-10,
    sigma::Float64 = 1.0, eta0::Float64 = 1e-10, zeta0::Float64 = 1e-10,
    H1::Int = 0, labels::Array{Int64,1} = Array{Int64,1}())

Returns an initialized structure of type vbmf_sparse_parameters.
"""
function vbmf_sparse_init(Y::Array{Float64,2}, H::Int; ca::Float64 = 1.0,
    alpha0::Float64 = 1e-10, beta0::Float64 = 1e-10, cb::Float64 = 1.0,
    gamma0::Float64 = 1e-10, delta0::Float64 = 1e-10,
    sigma::Float64 = 1.0, eta0::Float64 = 1e-10, zeta0::Float64 = 1e-10,
    H1::Int = 0, labels::Array{Int64,1} = Array{Int64,1}())
    params = vbmf_sparse_parameters()
    L, M = size(Y)

    params.L, params.M = L, M
    params.H = H
    params.MH = M*H
    params.H1 = H1
    params.labels = labels

    params.AHat = randn(M, H)
    # now set zeroes to where the columns of Y labeled as non-infected are
    # in A, these are rows
    params.AHat[labels, end-H1+1:end] = 0.0
    params.ATVecHat = reshape(params.AHat', M*H)
    params.SigmaATVec = eye(M*H, M*H)
    params.diagSigmaATVec = ones(M*H)
    params.invSigmaATVec = eye(M*H, M*H)
    params.SigmaA = zeros(H, H)

    params.BHat = randn(L, H)
    params.SigmaB = zeros(H, H)

    params.CA = ca*ones(M*H)
    params.alpha0 = alpha0
    params.beta0 = beta0
    params.alpha = params.alpha0 + 0.5
    params.beta = beta0*ones(M*H)

    params.CB = cb*ones(H)
    params.gamma0 = gamma0
    params.delta0 = delta0
    params.gamma = gamma0 + L/2
    params.delta = delta0*ones(H)

    params.sigmaHat = sigma
    params.eta0 = eta0
    params.zeta0 = zeta0
    params.eta = eta0 + L*M/2
    params.zeta = zeta0
    params.sigmaVecHat = sigma*ones(L)
    params.etaVec = (eta0 + M/2)*ones(L)
    params.zetaVec = zeta0*ones(L)

    params.YHat = params.BHat*params.AHat'
    params.trYTY = traceXTY(Y, Y)

    return params
end

"""
   copy(params_in::vbmf_sparse_parameters)

Copy function for vbmfa_sparse_parameters.
"""
function copy(params_in::vbmf_sparse_parameters)
    params = vbmf_sparse_parameters()

    for field in fieldnames(params_in)
        setfield!(params, field, copy(getfield(params_in, field)))
    end

    return params
end

"""
    updateA!(Y::Array{Float64,2}, params::vbmf_sparse_parameters; full_cov::Bool, diag_var::Bool = false)

Updates mean and covariance of vec(A^T) and also of the A matrix. If full_cov is true,
then inverse of full covariance matrix is computed, otherwise just the diagonal is estimated.
"""
function updateA!(Y::Array{Float64,2}, params::vbmf_sparse_parameters; full_cov::Bool = false, diag_var::Bool = false)
    # either just the diagonal or the full covariance
    if full_cov
        # compute the inverse of covariance
        if diag_var
            params.invSigmaATVec = kron(eye(params.M),
                (params.BHat'*diagm(params.sigmaVecHat)*params.BHat + params.L*mean(params.sigmaVecHat)*params.SigmaB)) + diagm(params.CA)
        else
            params.invSigmaATVec = params.sigmaHat*kron(eye(params.M),
                (params.BHat'*params.BHat + params.L*params.SigmaB)) + diagm(params.CA)
        end

        params.SigmaATVec = inv(params.invSigmaATVec)
        params.diagSigmaATVec = diag(params.SigmaATVec)

        # now the mean
        if diag_var
            params.ATVecHat = params.SigmaATVec*reshape(params.BHat'*diagm(params.sigmaVecHat)*Y, params.H*params.M)
        else
            params.ATVecHat = params.sigmaHat*params.SigmaATVec*reshape(params.BHat'*Y, params.H*params.M)
        end

        # now compute the covariance of AHat
        params.SigmaA = zeros(params.H, params.H)
        for m in 1:params.M
            params.SigmaA += view(params.SigmaATVec, (m-1)*params.H+1:m*params.H, (m-1)*params.H+1:m*params.H)
        end

    else # this updates just the diagonal of the covariance of A
        # compute the inverse of covariance
        # only use for large problems, for small problems the full covariance is more precise
        if diag_var
            for h in 1:params.H
                # instead of the full matrix multiplication as above, we compute just
                # the diagonal
                params.diagSigmaATVec[h] = norm2(params.BHat[:,h].*params.sigmaVecHat) + params.L*mean(params.sigmaVecHat)*params.SigmaB[h,h]
            end
        else         # this is computed in case the variance of data is homoscedastic
            for h in 1:params.H
                # instead of the full matrix multiplication as above, we compute just
                # the diagonal
                params.diagSigmaATVec[h] = params.sigmaHat*norm2(params.BHat[:,h]) + params.L*params.SigmaB[h,h]
            end
        end
        # we filled just the first H elements, so we just copy it to fill the rest
        params.diagSigmaATVec[params.H+1:end] = repeat(params.diagSigmaATVec[1:params.H], inner = params.M-1)
        # also add the CA vector
        params.diagSigmaATVec += params.CA

        # finally, invert it
        params.diagSigmaATVec = 1./params.diagSigmaATVec

        # now the mean
        if diag_var
            params.ATVecHat = params.diagSigmaATVec.*reshape(params.BHat'*diagm(params.sigmaVecHat)*Y, params.H*params.M)
        else
            params.ATVecHat = params.sigmaHat*params.diagSigmaATVec.*reshape(params.BHat'*Y, params.H*params.M)
        end

        # now compute the covariance of AHat
        params.SigmaA = zeros(params.H, params.H)
        for m in 1:params.M
            params.SigmaA += diagm(view(params.diagSigmaATVec, (m-1)*params.H+1:m*params.H))
        end
    end

    # now set zeroes to where the columns of Y labeled as non-infected are
    # in A, these are rows
    params.AHat = reshape(params.ATVecHat, params.H, params.M)'
    params.AHat[params.labels, end-params.H1+1:end] = 0.0
    params.ATVecHat = reshape(params.AHat', params.M*params.H)
end

"""
    updateB!(Y::Array{Float64,2}, params::vbmf_sparse_parameters, diag_var::Bool = false)

Updates mean and covariance of the B matrix.
"""
function updateB!(Y::Array{Float64,2}, params::vbmf_sparse_parameters;  diag_var::Bool = false)
    # first, compute the inverse of the covariance
    if diag_var
        params.SigmaB = diagm(params.CB) + mean(params.sigmaVecHat)*(params.AHat'*params.AHat + params.SigmaA)
        #invert it
        params.SigmaB = inv(params.SigmaB)
        # now compute the mean
        params.BHat = diagm(params.sigmaVecHat)*Y*params.AHat*params.SigmaB
    else
        params.SigmaB = diagm(params.CB) + params.sigmaHat*(params.AHat'*params.AHat + params.SigmaA)
        #invert it
        params.SigmaB = inv(params.SigmaB)
        params.BHat = params.sigmaHat*Y*params.AHat*params.SigmaB
    end
end

"""
    updateYHat!(params::vbmf_sparse_parameters)

Updates estimate of Y.
"""
function updateYHat!(params::vbmf_sparse_parameters)
    params.YHat = params.BHat*params.AHat'
end

"""
    updateCA!(params::vbmf_sparse_parameters)

Updates the estimate of CA.
"""
function updateCA!(params::vbmf_sparse_parameters)
    params.beta = params.beta0*ones(params.M*params.H) +
    1/2*(params.ATVecHat.*params.ATVecHat + params.diagSigmaATVec)
    params.CA = params.alpha*ones(params.M*params.H)./params.beta
end

"""
    updateCB!(params::vbmf_sparse_parameters)

Updates the estimate of CB.
"""
function updateCB!(params::vbmf_sparse_parameters)
    for h in 1:params.H
        params.delta[h] = params.delta0 + 1/2*(params.BHat[:,h]'*params.BHat[:,h])[1] + 1/2*params.SigmaB[h,h]
        params.CB[h] = params.gamma/params.delta[h]
    end
end

"""
    updateSigma!(Y::Array{Float64,2}, params::vbmf_sparse_parameters, diag_var::Bool = false)

Updates estimate of the measurement variance.
"""
function updateSigma!(Y::Array{Float64,2}, params::vbmf_sparse_parameters; diag_var::Bool = false)
    if diag_var
        for l in 1:params.L
            params.zetaVec[l] = params.zeta0 + 1/2*norm2(Y[l,:]) - sum(Y[l,:].*(params.AHat*params.BHat[l,:])) +
                1/2*traceXTY(params.AHat'*params.AHat + params.SigmaA,
                    params.BHat[l,:]*params.BHat[l,:]' + params.SigmaB)

            params.sigmaVecHat[l] = params.etaVec[l]/params.zetaVec[l]
        end
    else
        params.zeta = params.zeta0 + 1/2*params.trYTY - traceXTY(params.BHat, Y*params.AHat) +
            1/2*traceXTY(params.AHat'*params.AHat + params.SigmaA,
                params.BHat'*params.BHat + params.L*params.SigmaB)

        params.sigmaHat = params.eta/params.zeta
    end
end

"""
    vbmf_sparse!(Y::Array{Float64, 2}, params::vbmf_sparse_parameters, niter::Int, eps::Float64 = 1e-6,
    est_var = false, full_cov::Bool = false, logdir = "", desc = "")

Computes variational bayes matrix factorization of Y = AB' + E. Independence of A and B is assumed.
Estimate of variance sigma can be either a scalar or a vector - for row variance estimation.
ARD property is imposed upon columns of B and also
upon vec(A) through gamma prior CB and CA and estimation of covariance.
The prior model is following:

    p(Y|A,B) = N(Y|BA^T, 1/sigma*I) or p(Y|A,B) = N(Y|BA^T, inv(diag(sigma)))
    p(vec(A^T)) = N(vec(A^T)|0, invCA), CA = diag(c_a)
    p(B_h) = MN(B|0, I, invCB), CB = diag(c1, ..., cH)
    p(CA) = G(C_A| alpha0, beta0)
    p(CB) = G(C_B| gamma0, delta0)

The params argument with initialized data is modified and contains the resulting estimates after the
algorithm stops.
"""
function vbmf_sparse!(Y::Array{Float64, 2}, params::vbmf_sparse_parameters, niter::Int; eps::Float64 = 1e-6,
    diag_var::Bool = false, full_cov::Bool = false, logdir = "", desc = "", verb = false, est_cb::Bool = true)
    priors = Dict()

    # create the log dictionary
    loging = false
    if logdir !=""
        loging = true
        logVar = create_log(params)
    end

    # choice of convergence control variable
    convergence_var = :BHat
    #convergence_var = :dY

    if convergence_var == :dY
        old = Y - params.YHat
    else
        old = copy(getfield(params, convergence_var))
    end
    d = eps + 1.0 # delta
    i = 1

    # run the loop for a given number of iterations
    while ((i <= niter) && (d > eps))
        updateA!(Y, params, full_cov = full_cov, diag_var = diag_var)
        updateB!(Y, params, diag_var = diag_var)

        updateCA!(params)
        if est_cb
            updateCB!(params)
        end

        updateSigma!(Y, params, diag_var = diag_var)

        if loging
            update_log!(logVar, params)
        end

        # check convergence
        if convergence_var == :dY
            updateYHat!(params)
            d = delta(Y - params.YHat, old)
            old = Y - params.YHat
        else
            d = delta(getfield(params, convergence_var), old)
            old = copy(getfield(params, convergence_var))
        end
        i += 1
    end

    # finally, compute the estimate of Y
    updateYHat!(params)

    # convergence info
    if verb
        print("Factorization finished after ", i-1, " iterations, eps = ", d, "\n")
    end

    # save inputs and outputs
    if loging
        println("Saving outputs and inputs under ", logdir, "/")
        save_log(logVar, Y, priors, logdir, desc = desc)
    end

    return d
end

"""
    vbmf_sparse(Y::Array{Float64, 2}, params_in::vbmf_sparse_parameters, niter::Int, eps::Float64 = 1e-6,
    est_var = false, full_cov::Bool = false, logdir = "", desc = "")

Calls vbmf_sparse!() but copies the params_in argument so that it is not modified and can be reused.
"""
function vbmf_sparse(Y::Array{Float64, 2}, params_in::vbmf_sparse_parameters, niter::Int; eps::Float64 = 1e-6,
    diag_var::Bool = false, full_cov::Bool = false, logdir = "", desc = "", verb = false, est_cb::Bool = true)
    # make a copy of input params
    params = copy(params_in)

    # run the algorithm
    d = vbmf_sparse!(Y, params, niter, eps = eps, diag_var = diag_var, full_cov = full_cov,
        logdir = logdir, desc = desc, verb = verb, est_cb = est_cb)

    return params, d
end

"""
    lowerBound(Y::Array{Float64,2}, params::vbmf_sparse_parameters)

Compute the lower bound for logarithm of data distribution.
"""
function lowerBound(Y::Array{Float64,2}, params::vbmf_sparse_parameters)
    L = 0.0
    # E[lnp(Y|params)]
    L += - params.L*params.M/2*ln2pi + params.L*params.M/2*gammaELn(params.eta, params.zeta)
    L += - params.sigmaHat/2*(params.trYTY - 2*traceXTY(params.BHat, Y*params.AHat)
         + traceXTY(params.AHat'*params.AHat + params.SigmaA, params.BHat'*params.BHat + params.L*params.SigmaB))
    # E[lnp(vec(A'))]
    L += - params.MH/2*ln2pi + 1/2*sum(map(gammaELn, params.alpha*ones(size(params.beta)), params.beta))
    L += - (1/2*params.CA'*(params.ATVecHat.^2 + params.diagSigmaATVec))[1]
    # E[lnp(B)]
    L += - params.L*params.H/2*ln2pi
    L += params.L/2*sum(map(gammaELn, params.gamma*ones(size(params.delta)), params.delta))
    L += - 1/2*traceXTY(diagm(params.CB), (params.BHat'*params.BHat + params.L*params.SigmaB))
    # E[lnp(sigma)]
    L += params.eta0*log(params.zeta0) - lgamma(params.eta0)
    L += (params.eta0 - 1)*gammaELn(params.eta, params.zeta) - params.zeta0*params.sigmaHat
    # E[lnp(CA)]
    L += params.MH*(params.alpha0*log(params.beta0) - lgamma(params.alpha0))
    L += (params.alpha0 - 1)*sum(map(gammaELn, params.alpha*ones(size(params.beta)), params.beta))
    L += - params.beta0*sum(params.CA)
    # E[lnp(CB)]
    L += params.H*(params.gamma0*log(params.delta0) - lgamma(params.gamma0))
    L += (params.gamma0 - 1)*sum(map(gammaELn, params.gamma*ones(size(params.delta)), params.delta))
    L += - params.gamma0*sum(params.CB)

    # H(vec(A'))
    L += normalEntropy(params.diagSigmaATVec)
    # H(B)
    L += normalEntropy(kron(params.SigmaB,eye(params.L)))
    # H(sigma)
    L += gammaEntropy(params.eta, params.zeta)
    # H(CA)
    L += sum(map(gammaEntropy, params.alpha*ones(size(params.beta)), params.beta))
    # H(CB)
    L += sum(map(gammaEntropy, params.gamma*ones(size(params.delta)), params.delta))
    return L
end

"""
    lowerBoundTrimmed(Y::Array{Float64,2}, params_in::vbmf_sparse_parameters, trim = 1e-1)

AHat with abs(AHat[i,j]) < trim and the corresponding elements of CA are not involved in the computation.
"""
function lowerBoundTrimmed(Y::Array{Float64,2}, params_in::vbmf_sparse_parameters, trim = 1e-1)
    params = copy(params_in)

    trim_inds = (abs(params.ATVecHat) .> trim)
    params.ATVecHat = params.ATVecHat[trim_inds]
    params.MH = size(params.ATVecHat)[1]
    params.beta = params.beta[trim_inds]
    params.CA = params.CA[trim_inds]
    params.diagSigmaATVec = params.diagSigmaATVec[trim_inds]

    return lowerBound(Y, params)
end


text:
# by Lazaro Alonso
# using GLMakie # HIDE
using GLMakie
GLMakie.activate!() # HIDE
let
    x = 1:10
    y = 1:10
    z = 1:10
    f(x,y,z) = x^2 + y^2 + z^2
    vol = [f(ix,iy,iz) for ix in x, iy in y, iz in z]
    fig, ax, _ = volume(x, y, z, vol, colormap = :plasma, colorrange = (minimum(vol), maximum(vol)),
        figure = (; resolution = (800,800)),
        axis=(; type=Axis3, perspectiveness = 0.5,  azimuth = 7.19, elevation = 0.57,
            aspect = (1,1,1)))

    fig
    save(joinpath(@__DIR__, "output", "volume.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie"]) # HIDE

ext: jl

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
## using_cairomakie
using CairoMakie

## cairomakie_lineplot
f, ax, l = lines(1:10, 1:10)
mktempdir() do path
    save(joinpath(path, "test_cairomakie.png"), f)
end


text: # Ask how to manage the Julia development environments to [JinraeKim](https://github.com/JinraeKim)


text:
using WaterLily
using LinearAlgebra: norm2

function circle(radius=8;Re=250,n=12,m=6)
    center, ν = radius*m/2, radius/Re
    body = AutoBody((x,t)->norm2(x .- center) - radius)
    Simulation((n*radius+2,m*radius+2), [1.,0.], radius; ν, body)
end

using UnicodePlots,SparseArrays
function flood!(canvas,bitarray,color)
    x,y,_ = findnz(sparse(bitarray))
    points!(canvas,x,y,color)
end
function ascii_vort(sim,lines=10,body=false)
    a = sim.flow.σ; width,height = size(a)
    canvas = AsciiCanvas((2*lines*width)÷height,lines,width=width,height=height)
    @inside a[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
    flood!(canvas,a.<-1,:red)
    flood!(canvas,a.>1,:blue)
    !body && return canvas
    @inside a[I] = sum(WaterLily.ϕ(i,CartesianIndex(I,i),sim.flow.μ₀) for i ∈ 1:2)/2
    flood!(canvas,a.<0.25,:white)
    return canvas
end

ansi_moveup(n::Int) = string("\e[", n, "A")
ansi_enablecursor = "\e[?25h"
ansi_disablecursor = "\e[?25l"

function sim_ascii!(sim;duration=10,step=0.25,remeasure=false,lines=10,body=false)
    t₀ = round(sim_time(sim))
    print(ansi_disablecursor,repeat("\n",lines+1))
    for tᵢ in range(t₀,t₀+duration;step)
        sim_step!(sim,tᵢ;remeasure)
        println(ansi_moveup(lines),ascii_vort(sim,lines,body))
    end
    print(ansi_enablecursor)
end


text:
UNARY = [:+, :-, :~, :!, :abs, :sign, :sqrt, :cbrt,
          :log, :log2, :log10, :log1p,
          :exp, :exp2, :exp10, :expm1,
          :cos, :sin, :tan, :cosd, :sind, :tand,
          :acos, :asin, :atan, :acosd, :asind, :atand,
          :isnan, :isinf
        ]
MATH_DOTONLY    = [:.+, :.-, :.*, :./, :.%, :.^]
MATH_ALL        = [MATH_DOTONLY; [:+, :-, :*, :/, :%]]
COMPARE_DOTONLY = [:.>, :.<, :.==, :.>=, :.<=, :.!=]
BOOLEAN_OPS     = [:&; :|; :$; COMPARE_DOTONLY]

for op in [UNARY; MATH_ALL; BOOLEAN_OPS]
    eval(Expr(:import, :Base, op))
end # for

import Base.diff

###### Unary operators and functions ################

# TimeArray
for op in UNARY
    @eval begin
        function ($op){T,N}(ta::TimeArray{T,N})
            cnames  = [string($op) * name for name in ta.colnames]
            vals = ($op)(ta.values)
            TimeArray(ta.timestamp, vals, cnames, ta.meta)
        end # function
    end # eval
end # loop

###### Numerical operations and comparisons #########

# TimeArray <--> Scalar
for op in [MATH_ALL; COMPARE_DOTONLY]
    @eval begin
        function ($op){T<:Number,N}(ta::TimeArray{T,N}, var::Number)
            cnames  = [name * string($op) * string(var) for name in ta.colnames]
            vals = ($op)(ta.values, var)
            TimeArray(ta.timestamp, vals, cnames, ta.meta)
        end # function
    end # eval
end # loop

# Scalar <--> TimeArray
for op in [MATH_ALL; COMPARE_DOTONLY]
    @eval begin
        function ($op){T<:Number,N}(var::Number, ta::TimeArray{T,N})
            cnames  = [string(var) * string($op) * name for name in ta.colnames]
            vals = ($op)(var, ta.values)
            TimeArray(ta.timestamp, vals, cnames, ta.meta)
        end # function
    end # eval
end # loop

# ND TimeArray <--> MD TimeArray
for op in [MATH_DOTONLY; COMPARE_DOTONLY]
    @eval begin
        function ($op){S<:Number,T<:Number,N,M}(ta1::TimeArray{S,N}, ta2::TimeArray{T,M})

            # first test metadata matches
            meta = ta1.meta == ta2.meta ? ta1.meta : Void

            # determine array widths and name cols accordingly
            w1, w2  = length(ta1.colnames), length(ta2.colnames)
            if w1 == w2
              cnames = [ta1.colnames[i]*string($op)*ta2.colnames[i] for i=1:w1]
            elseif w1==1
              cnames = [ta1.colnames[1]*string($op)*ta2.colnames[i] for i=1:w2]
            elseif w2==1
              cnames = [ta1.colnames[i]*string($op)*ta2.colnames[1] for i=1:w1]
            else
              error("arrays must have the same number of columns, or one must be a single column")
            end

            # obtain shared timestamp
            idx1, idx2 = overlaps(ta1.timestamp, ta2.timestamp)
            tstamp = ta1[idx1].timestamp

            # retrieve values that match the Int array matching dates
            vals1, vals2 = ta1[idx1].values, ta2[idx2].values

            # compute output values
            vals = ($op)(vals1, vals2)
            TimeArray(tstamp, vals, cnames, meta)

        end # function
    end # eval
end # loop

###### Boolean operations and comparisons ###############

# TimeArray <--> Bool
for op in BOOLEAN_OPS
    @eval begin
        function ($op){N}(ta::TimeArray{Bool,N}, var::Bool)
            cnames = [name * string($op) * string(var) for name in ta.colnames]
            vals   = ($op)(ta.values, var)
            TimeArray(ta.timestamp, vals, cnames, ta.meta)
        end # function
    end # eval
end # loop

# Bool <--> TimeArray
for op in BOOLEAN_OPS
    @eval begin
        function ($op){N}(var::Bool, ta::TimeArray{Bool,N})
            cnames = [string(var) * string($op) * name for name in ta.colnames]
            vals   = ($op)(var, ta.values)
            TimeArray(ta.timestamp, vals, cnames, ta.meta)
        end # function
    end # eval
end # loop

# Boolean ND TimeArray <--> Boolean ND TimeArray
# Doesn't support broadcasting since it isn't supported by Base either
for op in BOOLEAN_OPS
    @eval begin
        function ($op){N}(ta1::TimeArray{Bool,N}, ta2::TimeArray{Bool,N})
            # test column count matches
            length(ta1.colnames) == length(ta2.colnames) ?
              ncols=length(ta1.colnames) :
              error("arrays must have the same number of columns")
            # test metadata matches
            ta1.meta == ta2.meta ? meta = ta1.meta : error("metadata doesn't match")
            cnames = [ta1.colnames[i]*string($op)*ta2.colnames[i] for i=1:ncols]
            idx1, idx2 = overlaps(ta1.timestamp, ta2.timestamp)
            # obtain shared timestamp
            tstamp = ta1[idx1].timestamp
            # retrieve values that match the Int array matching dates
            vals1  = ta1[idx1].values
            vals2  = ta2[idx2].values
            vals = ($op)(vals1, vals2)
            TimeArray(tstamp, vals, cnames, meta)
        end # function
    end # eval
end # loop

###### lag, lead ################

function lag{T,N}(ta::TimeArray{T,N}, n::Int=1; padding::Bool=false, period::Int=0)

    if period != 0
      warn("the period kwarg is deprecated, use lag(ta::TimeArray, period::Int) instead")
      n = period
    end #if

    if padding
        paddedvals = [NaN*ones(n, length(ta.colnames)); ta.values[1:end-n, :]]
        ta = TimeArray(ta.timestamp, paddedvals, ta.colnames, ta.meta)
    else
        ta = TimeArray(ta.timestamp[1+n:end], ta.values[1:end-n, :], ta.colnames, ta.meta)
    end #if

    N == 1 && (ta = ta[ta.colnames[1]])
    return ta

end #lag

function lead{T,N}(ta::TimeArray{T,N}, n::Int=1; padding::Bool=false, period::Int=0)

    if period != 0
      warn("the period kwarg is deprecated, use lead(ta::TimeArray, period::Int) instead")
      n = period
    end #if

    if padding
        paddedvals = [ta.values[1+n:end, :]; NaN*ones(n, length(ta.colnames))]
        ta = TimeArray(ta.timestamp, paddedvals, ta.colnames, ta.meta)
    else
        ta = TimeArray(ta.timestamp[1:end-n], ta.values[1+n:end, :], ta.colnames, ta.meta)
    end #if

    N == 1 && (ta = ta[ta.colnames[1]])
    return ta

end #lead

###### diff #####################

# TODO: Support higher-order differencing?
function diff(ta::TimeArray; padding::Bool=false)
    cols = ta.colnames
    ta = ta .- lag(ta, padding=padding)
    ta.colnames[:] = cols
    return ta
end #diff

###### percentchange ############

function percentchange(ta::TimeArray, returns::Symbol=:simple; padding::Bool=false, method::AbstractString="")

    if method != ""
        warn("the method kwarg is deprecated, use percentchange(ta, :methodname) instead")
        returns = symbol(method)
    end #if

    cols = ta.colnames
    ta =  returns == :log ? diff(log(ta), padding=padding) :
          returns == :simple ? expm1(percentchange(ta, :log, padding=padding)) :
          error("returns must be either :simple or :log")
    ta.colnames[:] = cols

   return ta

end #percentchange

###### moving ###################

function moving{T}(ta::TimeArray{T,1}, f::Function, window::Int; padding::Bool=false)
    tstamps = padding ? ta.timestamp : ta.timestamp[window:end]
    vals    = zeros(ta.values[window:end])
    for i=1:length(vals)
        vals[i] = f(ta.values[i:i+(window-1)])
    end
    padding && (vals = [NaN*ones(window-1); vals])
    TimeArray(tstamps, vals, ta.colnames, ta.meta)
end

function moving{T}(ta::TimeArray{T,2}, f::Function, window::Int; padding::Bool=false)
    tstamps = padding ? ta.timestamp : ta.timestamp[window:end]
    vals    = zeros(ta.values[window:end, :])
    for i=1:size(vals,1), j=1:size(vals, 2)
        vals[i, j] = f(ta.values[i:i+(window-1), j])
    end
    padding && (vals = [NaN*ones(ta.values[1:(window-1), :]); vals])
    TimeArray(tstamps, vals, ta.colnames, ta.meta)
end

###### upto #####################

function upto{T}(ta::TimeArray{T,1}, f::Function)
    vals = zeros(ta.values)
    for i=1:length(vals)
        vals[i] = f(ta.values[1:i])
    end
    TimeArray(ta.timestamp, vals, ta.colnames, ta.meta)
end

function upto{T}(ta::TimeArray{T,2}, f::Function)
    vals = zeros(ta.values)
    for i=1:size(vals, 1), j=1:size(vals, 2)
        vals[i, j] = f(ta.values[1:i, j])
    end
    TimeArray(ta.timestamp, vals, ta.colnames, ta.meta)
end

###### basecall #################

basecall{T,N}(ta::TimeArray{T,N}, f::Function; cnames=ta.colnames) =  TimeArray(ta.timestamp, f(ta.values), cnames, ta.meta)

###### uniform observations #####

function uniformspaced(ta::TimeArray)
    gap1 = ta.timestamp[2] - ta.timestamp[1]
    i, n, is_uniform = 2, length(ta), true
    while is_uniform & (i < n)
        is_uniform = gap1 == (ta.timestamp[i+1] - ta.timestamp[i])
        i += 1
    end #while
    return is_uniform
end #uniformlyspaced

function uniformspace{T,N}(ta::TimeArray{T,N})
    min_gap = minimum(ta.timestamp[2:end] - ta.timestamp[1:end-1])
    newtimestamp = ta.timestamp[1]:min_gap:ta.timestamp[end]
    emptyta = TimeArray(collect(newtimestamp), zeros(length(newtimestamp), 0), String[], ta.meta)
    ta = merge(emptyta, ta, :left)
    N == 1 && (ta = ta[ta.colnames[1]])
    return ta
end #uniformlyspace

###### dropnan ####################

dropnan(ta::TimeArray, method::Symbol=:all) =
    method == :all ? ta[find(any(!isnan(ta.values), 2))] :
    method == :any ? ta[find(all(!isnan(ta.values), 2))] :
    error("dropnan method must be :all or :any")


text:
using RCall, MixedModels, Test
using StatsBase: zscore
using Tables: columntable
using GLM
import GLM: Link
const LMM = LinearMixedModel
const GLMM = GeneralizedLinearMixedModel

logistic(x)  = 1 / (1 + exp(-x))

@testset "glmerMod" begin
    reval("""
    if(!("lme4" %in% installed.packages())){
        # use tmp for tests if lme4 isn't available
        .libPaths("/tmp")
        lib <- .libPaths()[1L]
        install.packages("lme4",repos="https://cloud.r-project.org", libs=lib)
    }
    if(!("afex" %in% installed.packages())){
        .libPaths("/tmp")
        lib <- .libPaths()[1L]
        install.packages("afex",repos="https://cloud.r-project.org", libs=lib)
    }
    """)

    ### from R ###
    @testset "get glmerMod" begin
        @testset "Binomial" begin

            reval(raw"""
            attach(cbpp)
            cbpp$rate <- with(cbpp, incidence/size)
            rlmm <- glmer(rate ~ period + (1 | herd), weights = size,
                          family = binomial, data = cbpp)
            """)
            rlmm = rcopy(R"rlmm")
            @test rcopy(R"fitted(rlmm)") ≈ fitted(rlmm) atol=0.001
            @test_broken rcopy(R"-2 * logLik(rlmm)") ≈ deviance(rlmm) atol=0.001
            # this is where the weights have their biggest impact
            @test stderror(rlmm) ≈ rcopy(R"""coef(summary(rlmm))[,"Std. Error"]""") atol=0.1
            @test rcopy(R"logLik(rlmm)") ≈ loglikelihood(rlmm) atol=0.001

            @testset "alternative response specifications" begin
                # R"""
                # # from lme4 ?cbpp
                # ## response as a matrix
                # m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                #               family = binomial, data = cbpp, nAGQ=0)
                #  ## response as a vector of probabilities and usage of argument "weights"
                #  m1p <- glmer(incidence / size ~ period + (1 | herd), weights = size,
                #               family = binomial, data = cbpp)
                # """

                @testset "cbind" begin
                end

                @testset "proportion computed in line" begin
                end
            end

            @testset "contrasts" begin
                reval(raw"""
                contrasts(cbpp$period) <- contr.helmert(levels(cbpp$period))
                """)
                rlmm = rcopy(R"fm1 <- glmer(rate ~ period + (1 | herd), weights = size,family = binomial, data = cbpp)");
                @test fixef(rlmm) ≈  rcopy(R"fixef(fm1)");
            end

            @testset "Bernoulli" begin
                @testset "nAGQ and scalar RE" begin
                    reval("""
                    rlmm <- glmer(r2 ~ Anger + Gender + btype + situ + (1|id),
                                        family = binomial, data=VerbAgg, nAGQ=4)
                    """)

                    rlmm = rcopy(R"rlmm")
                    @test rcopy(R"fitted(rlmm)") ≈ fitted(rlmm) atol=0.001
                    @test rcopy(R"-2 * logLik(rlmm)") ≈ deviance(rlmm) atol=0.001
                    # this is where the weights have their biggest impact
                    @test stderror(rlmm) ≈ rcopy(R"""coef(summary(rlmm))[,"Std. Error"]""") atol=0.05
                    @test_broken rcopy(R"logLik(rlmm)") ≈ loglikelihood(rlmm) atol=0.001
                end

                @testset "Laplace and probit link" begin
                    # TODO add tests for each link
                    reval("""
                    rlmm <- glmer(r2 ~ Anger + Gender + btype + situ + (1|id),
                                        family = binomial(link="probit"), data=VerbAgg)
                    """)

                    rlmm = rcopy(R"rlmm")
                    @test Link(rlmm.resp) isa GLM.ProbitLink
                    @test rcopy(R"fitted(rlmm)") ≈ fitted(rlmm) atol=0.001
                    @test_broken rcopy(R"-2 * logLik(rlmm)") ≈ deviance(rlmm) atol=0.001
                    # this is where the weights have their biggest impact
                    @test stderror(rlmm) ≈ rcopy(R"""coef(summary(rlmm))[,"Std. Error"]""") atol=0.05
                    @test rcopy(R"logLik(rlmm)") ≈ loglikelihood(rlmm) atol=0.001
                end
            end

        end

        @testset "Poisson and fast fit" begin
            reval("""
            # from the lme4 docs
            form <- TICKS ~ YEAR + HEIGHT + (1|BROOD) + (1|INDEX) + (1|LOCATION)
            rlmm  <- glmer(form, family="poisson",data=grouseticks, nAGQ=0)
            """)
            rlmm = rcopy(R"rlmm")
            @test rcopy(R"fitted(rlmm)") ≈ fitted(rlmm) atol=0.001
            @test_broken rcopy(R"-2 * logLik(rlmm)") ≈ deviance(rlmm) atol=0.001
            # 7.5% difference isn't great....
            @test stderror(rlmm) ≈ rcopy(R"""coef(summary(rlmm))[,"Std. Error"]""") rtol=0.075
            @test rcopy(R"logLik(rlmm)") ≈ loglikelihood(rlmm) atol=0.05
        end

    end

    ### from Julia ###
    @testset "put glmerMod" begin
        @testset "Bernoulli" begin
            # NOTE: For Bernoulli, some of the tests are doubled in order to
            #       check nAQG and fast options
            # TODO: remove upon next MixedModels.jl release
            dat = dataset(:verbagg)

            @testset "unfitted model" begin
                # we max out at 1 scalar RE for the nAGQ tests
                jlmm = GLMM(@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)),
                            dat, Bernoulli());
                jm = Tuple([jlmm, dat]);
                # unfitted model
                @test_throws ArgumentError @rput jm
            end

            @testset "nAGQ" begin
                # we max out at 1 scalar RE for the nAGQ tests
                jlmm = GLMM(@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)),
                            dat, Bernoulli());
                fit!(jlmm; fast=false, nAGQ=9)
                jm = (jlmm, dat);
                @rput jm;
                # @test_warn Regex(".*categorical.*") @rput jm;

                @test  rcopy(R"""jm@devcomp$dims["nAGQ"]""") == jlmm.optsum.nAGQ
                # note the really high tolerances
                @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) atol=0.1
                @test rcopy(R"-2 * logLik(jm)") ≈ deviance(jlmm) atol=0.1
            end

            @testset "Laplace" begin
                # we max out at 1 scalar RE for the nAGQ tests
                jlmm = GLMM(@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)),
                            dat, Bernoulli());
                fit!(jlmm);
                jm = Tuple([jlmm, dat])
                @rput jm;

                # note the really high tolerances
                @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) atol=0.1
                @test rcopy(R"-2 * logLik(jm)") ≈ deviance(jlmm) atol=0.5;
            end

            @testset "columntable" begin
                jlmm = GLMM(@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)),
                            dat, Bernoulli());
                fit!(jlmm);
                jm = Tuple([jlmm, columntable(dat)]);
                @rput jm;
            end
        end

        @testset "Binomial" begin
            # TODO: remove upon next MixedModels.jl release
            dat  = dataset(:cbpp)
            dat.rate = dat.incid ./ dat.hsz;
            jlmm = fit(MixedModel, @formula(rate ~ 1 + period + (1|herd)),
                      dat, Binomial(); wts=float(dat.hsz), fast=true)

            jm = (jlmm, dat)
            @rput jm;
            # @test_warn Regex(".*categorical.*") @rput jm;
            @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) atol=0.001
            @test_broken rcopy(R"-2 * logLik(jm)") ≈ deviance(jlmm) atol=0.001
            # this is where the weights have their biggest impact
            @test stderror(jlmm) ≈ rcopy(R"""coef(summary(jm))[,"Std. Error"]""") atol=0.1
            @test rcopy(R"logLik(jm)") ≈ loglikelihood(jlmm) atol=0.001
        end

        @testset "Poisson and fast fit" begin
            # TODO: remove upon next MixedModels.jl release
            center(v::AbstractVector) = v .- (sum(v) / length(v))
            dat = dataset(:grouseticks);
            dat.ch = center(dat.height);
            jlmm = GLMM(@formula(ticks ~ 1+year+ch+ (1|index) + (1|brood) + (1|location)),
                        dat, Poisson());

            jm = (jlmm, dat);
            fit!(jlmm, fast=true) # problems with this one in fast=false
            @rput jm;
            # @test_warn Regex(".*categorical.*") @rput jm;
            @test rcopy(R"""jm@devcomp$dims["nAGQ"]""") == 0;
            # note the really high tolerances
            @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) atol=0.001
            @test_broken rcopy(R"-2 * logLik(jm)") ≈ deviance(jlmm) atol=0.5
            @test rcopy(R"logLik(jm)") ≈ loglikelihood(jlmm) atol=0.001
        end

        @testset "Gaussian" begin
            # TODO: remove upon next MixedModels.jl release
            # TODO: dispersion check
            # dat = dataset(:verbagg);
            # jlmm = GLMM(@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item)),
            #             dat, Bernoulli());
            #
            # jm = (jlmm, dat);
            # fit!(jlmm, fast=true)
            # @rput jm;
            # # @test_warn Regex(".*categorical.*") @rput jm;
            # # note the really high tolerances
            # @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) atol=0.1
            # @test rcopy(R"-2 * logLik(jm)") ≈ deviance(jlmm) atol=0.5
        end

        @testset "InverseGaussian" begin
        end

        @testset "Gamma" begin
        end

        @testset "contrasts" begin
            dat = dataset(:verbagg);

            jlmm = fit(MixedModel, @formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item)),
                        dat, Bernoulli(), contrasts=Dict(:gender => EffectsCoding(),
                                                         :btypes => EffectsCoding()));
            jm = (jlmm, dat);
            @rput jm;
            @test fixef(jlmm) ≈ rcopy(R"fixef(jm)");
        end
    end
end


text:
module MCVI

using JSON
using POMDPs
using ParticleFilters
using Random
using Printf
using Distributed
using LinearAlgebra
import Statistics
import POMDPs: solve, action, rand, simulate, updater, initialize_belief, update

const Reward = Float64

# implementation warning handled by POMDPrequire
function init_lower_action end

include("subspace.jl")
include("belief.jl")
include("policy.jl")
include("simulate.jl")
include("alphaedge.jl")
include("updater.jl")
include("solver.jl")

export MCVISolver,
       MCVIPolicy,

       solve,
       action,
       MCVIUpdater,
       MCVISimulator,
       updater,
       dump_json,
       upper_bound,
       lower_bound

end # module


text:
import Term.Style: apply_style

@testset "\e[34mMACROS" begin
    @test (@green "string") == apply_style("{green}string{/green}")
    @test (@blue "string") == apply_style("{blue}string{/blue}")
    @test (@red "string") == apply_style("{red}string{/red}")

    @test (@bold "string") == apply_style("{bold}string{/bold}")
    @test (@italic "string") == apply_style("{italic}string{/italic}")
    @test (@underline "string") == apply_style("{underline}string{/underline}")

    @test (@style "string" bold underline) ==
          apply_style("{bold underline}string{/bold underline}")
    @test (@style "string" bold greeen) == apply_style("{bold greeen}string{/bold greeen}")
    @test (@style "string" red on_blue underline) ==
          apply_style("{red on_blue underline}string{/red on_blue underline}")
    @test (@style "string" bold italic) == apply_style("{bold italic}string{/bold italic}")

    # test with interpolations
    x = "test"
    @test (@red "aa $x aa") == "\e[31maa test aa\e[39m"
    @test (@bold "aa $x aa") == "\e[1maa test aa\e[22m"
    @test (@style "aa $x aa" green bold) == "\e[1m\e[32maa test aa\e[22m\e[39m"
end


text:
using BinDeps
import CondaBinDeps

@BinDeps.setup
libnetcdf = library_dependency(
    "libnetcdf",
    aliases = ["libnetcdf4", "libnetcdf-7", "netcdf"],
)

#CondaBinDeps.Conda.add_channel("conda-forge")
provides(CondaBinDeps.Manager, "libnetcdf", libnetcdf)
provides(AptGet, "libnetcdf-dev", libnetcdf, os = :Linux)
provides(Yum, "netcdf-devel", libnetcdf, os = :Linux)

@BinDeps.install Dict(:libnetcdf => :libnetcdf)


text:
if !@isdefined plotyes
    plotyes = true
end
if plotyes
    @eval using PyPlot
end
using WhereTheWaterFlows
if !(@isdefined WWF)
    const WWF = WhereTheWaterFlows
end

"A artificial DEM"
function ele(x, y; withpit=false, randfac=0.0)
    out = - (x^2 - 1)^2 - (x^2*y - x - 1)^2 + 6 + 0.1*x + 3*y
    if withpit
        out -= 2*exp(-(x^2 + y^2)*50)
    end
    out += randfac*randn()
    return out<0 ? 0.0 : out
end
dx = 0.01
xs = -1.5:dx:1
ys = -0.5:dx:3.0
dem = ele.(xs, ys', randfac=0.1, withpit=true);
plotyes && WWF.heatmap(xs,ys,dem)
area, slen, dir, nout, nin, pits, c, bnds  = WWF.waterflows(dem, drain_pits=true);
@assert size(dem)==(length(xs), length(ys))
plotyes && WWF.plotit(xs,ys,dem)
plotyes && WWF.plotarea(xs, ys, area, pits)

plotyes && WWF.heatmap(xs, ys, c)

demf = WWF.fill_dem(dem, pits, dir) #, small=1e-6)
plotyes && WWF.heatmap(xs,ys, demf.-dem)


text:
export batchsize,lr,δ,momentum,λ,▽max,order
const batchsize=50
const lr=1e-2 #learning rate
const δ=0.  #decrease constant lr=lr/(1+δ*t)
const momentum=0.
const λ=0. #regularization coefficient
const ▽max=0.1 #largest gradient
global order
global sharedParams=[]


text:
#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")

import GZip
import JSON

function iterlines(fun, filename; maxlines=typemax(Int))
    if endswith(filename, ".gz")
        f = GZip.open(filename)
        i = 0
        while !eof(f) && i < maxlines
            line = readline(f)
            i += 1
            if length(line) == 0
                continue
            end
            fun(line)
        end
        close(f)
    else
        open(filename) do f
            i = 0
            while !eof(f) && i < maxlines
                line = readline(f)
                i += 1
                fun(line)
            end
        end
    end
end

function parsetweet(line)
    if line[1] == '{'
        tweet = JSON.parse(line)
    else
        key, value = split(line, '\t', limit=2)
        tweet = JSON.parse(value)
        tweet["key"] = key
    end

    tweet
end

function itertweets(fun, filename::String; maxlines=typemax(Int))
    iterlines(filename, maxlines=maxlines) do line
        tweet = parsetweet(line)
        fun(tweet)
    end
end

function itertweets(fun, file; maxlines=typemax(Int))
    i = 0
    while !eof(file) && i < maxlines
        line = readline(file)
        i += 1
        try
            tweet = parsetweet(line)
            fun(tweet)
        catch
            continue
        end
    end
end

function getkeypath(dict, key)
    if contains(key, ".")
        v = dict
        for k in split(key, '.')
            v = v[k]
        end

        return v
    else
        return dict[key]
    end
end

function savepath{T}(filename::String, obj::T)
    d = dirname(filename)
    !isdir(d) && mkpath(d)

    open(filename, "w") do f
        save(f, obj)
    end
end

function loadpath{T}(filename::String, ::Type{T})
    return open(filename) do f
        load(f, T)
    end
end

function loaddatabase(filename)
    db = DocumentType[]
    if filename == "-"
        while !eof(f)
            push!(db, load(STDIN, DocumentType))
        end
    else
        open(filename) do f
            while !eof(f)
                push!(db, load(f, DocumentType))
            end
        end
    end

    db
end

function loadindex(args)
    local index
    env = args["env"]
    indexpath = joinpath(env, "index")
    if isfile(indexpath)
        open(indexpath) do f
            index = load(f, LocalSearchIndex, DocumentType, CosineDissimilarity(), loaditems=false)
        end

        open(joinpath(env, "db")) do f
            load(f, index)
        end
        return index
    else
        index = LocalSearchIndex(DocumentType, CosineDissimilarity(), recall=args["recall"], k=args["knn"])
    end

    index.neighborhood_algo = LogSatNeighborhood(1.1)
    index
end

function saveindex(args, index, fastlinks)
    env = args["env"]
    open(joinpath(env, "index"), "w") do f
        save(f, index, saveitems=false)
    end

    open(joinpath(env, "fastlinks"), "w") do f
        for fastlink in fastlinks
            save(f, fastlink)
        end
    end
end

max_stars_repo_name: INGEOTEC/ranix

max_issues_repo_name: INGEOTEC/ranix

max_forks_repo_name: INGEOTEC/ranix


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

export cooketriplet, doubleconvexlensonly

#Hardcode these glass types so this code will work even if it is not possible to download the glasses at build time.
const Examples_N_BK7 = GlassCat.Glass("$(@__MODULE__).Examples_N_BK7", 2, 1.03961212, 0.00600069867, 0.231792344, 0.0200179144, 1.01046945, 103.560653, 0.0, 0.0, NaN, NaN, 0.3, 2.5, 1.86e-6, 1.31e-8, -1.37e-11, 4.34e-7, 6.27e-10, 0.17, 20.0, -0.0009, 2.3, 1.0, 7.1, 1.0, 1, 1.0, [(0.3, 0.05, 25.0), (0.31, 0.25, 25.0), (0.32, 0.52, 25.0), (0.334, 0.78, 25.0), (0.35, 0.92, 25.0), (0.365, 0.971, 25.0), (0.37, 0.977, 25.0), (0.38, 0.983, 25.0), (0.39, 0.989, 25.0), (0.4, 0.992, 25.0), (0.405, 0.993, 25.0), (0.42, 0.993, 25.0), (0.436, 0.992, 25.0), (0.46, 0.993, 25.0), (0.5, 0.994, 25.0), (0.546, 0.996, 25.0), (0.58, 0.995, 25.0), (0.62, 0.994, 25.0), (0.66, 0.994, 25.0), (0.7, 0.996, 25.0), (1.06, 0.997, 25.0), (1.53, 0.98, 25.0), (1.97, 0.84, 25.0), (2.325, 0.56, 25.0), (2.5, 0.36, 25.0)], 1.5168, 2.3, 0.0, 0, 64.17, 0, 2.51, 0.0)

const Examples_N_SK16 =GlassCat.Glass("$(@__MODULE__).Examples_N_SK16", 2, 1.34317774, 0.00704687339, 0.241144399, 0.0229005, 0.994317969, 92.7508526, 0.0, 0.0, NaN, NaN, 0.31, 2.5, -2.37e-8, 1.32e-8, -1.29e-11, 4.09e-7, 5.17e-10, 0.17, 20.0, -0.0011, 3.2, 1.4, 6.3, 4.0, 1, 53.3, [(0.31, 0.02, 25.0), (0.32, 0.11, 25.0), (0.334, 0.4, 25.0), (0.35, 0.7, 25.0), (0.365, 0.86, 25.0), (0.37, 0.89, 25.0), (0.38, 0.93, 25.0), (0.39, 0.956, 25.0), (0.4, 0.97, 25.0), (0.405, 0.974, 25.0), (0.42, 0.979, 25.0), (0.436, 0.981, 25.0), (0.46, 0.984, 25.0), (0.5, 0.991, 25.0), (0.546, 0.994, 25.0), (0.58, 0.994, 25.0), (0.62, 0.993, 25.0), (0.66, 0.994, 25.0), (0.7, 0.996, 25.0), (1.06, 0.995, 25.0), (1.53, 0.973, 25.0), (1.97, 0.88, 25.0), (2.325, 0.54, 25.0), (2.5, 0.26, 25.0)], 1.62041, 3.3, 4.0, 0, 60.32, 0, 3.58, 0.0)

const Examples_N_SF2 = GlassCat.Glass("$(@__MODULE__).Examples_N_SF2", 2, 1.47343127, 0.0109019098, 0.163681849, 0.0585683687, 1.36920899, 127.404933, 0.0, 0.0, NaN, NaN, 0.365, 2.5, 3.1e-6, 1.75e-8, 6.62e-11, 7.51e-7, 8.99e-10, 0.277, 20.0, 0.0081, 1.0, 1.4, 6.68, 1.0, 1, 1.0, [(0.365, 0.007, 25.0), (0.37, 0.06, 25.0), (0.38, 0.4, 25.0), (0.39, 0.68, 25.0), (0.4, 0.83, 25.0), (0.405, 0.865, 25.0), (0.42, 0.926, 25.0), (0.436, 0.949, 25.0), (0.46, 0.961, 25.0), (0.5, 0.975, 25.0), (0.546, 0.986, 25.0), (0.58, 0.987, 25.0), (0.62, 0.984, 25.0), (0.66, 0.984, 25.0), (0.7, 0.987, 25.0), (1.06, 0.997, 25.0), (1.53, 0.984, 25.0), (1.97, 0.93, 25.0), (2.325, 0.76, 25.0), (2.5, 0.67, 25.0)], 1.64769, 1.2, 0.0, 0, 33.82, 0, 2.718, 0.0)

const Examples_N_SF14 = GlassCat.Glass("$(@__MODULE__).Examples_N_SF14", 2, 1.69022361, 0.0130512113, 0.288870052, 0.061369188, 1.7045187, 149.517689, 0.0, 0.0, NaN, NaN, 0.365, 2.5, -5.56e-6, 7.09e-9, -1.09e-11, 9.85e-7, 1.39e-9, 0.287, 20.0, 0.013, 1.0, 2.2, 9.41, 1.0, 1, 1.0, [(0.365, 0.004, 25.0), (0.37, 0.04, 25.0), (0.38, 0.33, 25.0), (0.39, 0.61, 25.0), (0.4, 0.75, 25.0), (0.405, 0.79, 25.0), (0.42, 0.87, 25.0), (0.436, 0.91, 25.0), (0.46, 0.938, 25.0), (0.5, 0.964, 25.0), (0.546, 0.983, 25.0), (0.58, 0.987, 25.0), (0.62, 0.987, 25.0), (0.66, 0.987, 25.0), (0.7, 0.985, 25.0), (1.06, 0.998, 25.0), (1.53, 0.98, 25.0), (1.97, 0.88, 25.0), (2.325, 0.64, 25.0), (2.5, 0.57, 25.0)], 1.76182, 1.0, 0.0, 0, 26.53, 0, 3.118, 0.0)

"""
    hemisphere()

Create a geometric hemisphere
"""
function hemisphere()::CSGTree
    sph = Sphere(10.0)
    pln = Plane(0.0, 0.0, -1.0, 0.0, 0.0, 0.0)
    # CSV operations create a csggenerator which instantiates the csg tree after applying a rigid body transformation.
    # This allows you to make as many instances of the object as you want with different transformations. We just want
    # the CSGTree object rather than a generator.
    return (sph ∩ pln)()
end

"""
    opticalhemisphere()

Create an optical hemisphere that has optical material properties so it will reflect and refract light. In the previous
example the hemisphere object had optical properties of Air, which is the default optical interface, so it won't refract
or reflect light.
"""
function opticalhemisphere()::CSGOpticalSystem
    sph = Sphere(10.0, interface = FresnelInterface{Float64}(Examples_N_BK7, Air))
    pln = Plane(0.0, 0.0, -1.0, 0.0, 0.0, 0.0, interface = FresnelInterface{Float64}(Examples_N_BK7, Air))
    assy = LensAssembly{Float64}((sph ∩ pln)())
    return CSGOpticalSystem(assy, Rectangle(1.0, 1.0, SVector{3,Float64}(0.0, 0.0, 1.0), SVector{3,Float64}(0.0, 0.0, -11.0)))
end

function cooketriplet(::Type{T} = Float64, detpix::Int = 1000) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Standard", "Stop", "Standard", "Standard", "Image"],
            Radius = [Inf, 26.777, 66.604, -35.571, 35.571, 35.571, -26.777, Inf],
            # OptimizeRadius = [false, true, true, true, true, true, true, false],
            Thickness = [Inf, 4.0, 2.0, 4.0, 2.0, 4.0, 44.748, missing],
            # OptimizeThickness = [false, true, true, true, true, true, true, false],
            Material = [Air, Examples_N_SK16, Air, Examples_N_SF2, Air, Examples_N_SK16, Air, missing],
            SemiDiameter = [Inf, 8.580, 7.513, 7.054, 6.033, 7.003, 7.506, 15.0]
        ),
        detpix,
        detpix
    )
end

function cooketripletfirstelement(::Type{T} = Float64) where {T<:Real}
    AxisymmetricOpticalSystem(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf, -35.571, 35.571, Inf],
            Thickness = [Inf, 4.0, 44.748, missing],
            Material = [Air, Examples_N_SK16, Air, missing],
            SemiDiameter = [Inf, 7.054, 6.033, 15.0]
        )
    )
end

function convexplano(::Type{T} = Float64) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf, 60.0, Inf, Inf],
            Thickness = [Inf, 10.0, 57.8, missing],
            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf, 9.0, 9.0, 15.0]
        )
    )
end

function doubleconvex(frontradius::T, rearradius::T) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [convert(T, Inf64), frontradius, rearradius, convert(T, Inf64)],
            # OptimizeRadius = [false, true, true, false],
            Thickness = [convert(T, Inf64), convert(T, 10.0), convert(T, 57.8), missing],
            # OptimizeThickness = [false, false, false, false],

            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [convert(T, Inf64), convert(T, 9.0), convert(T, 9.0), convert(T, 15.0)]
        )
    )
end

function doubleconvexconic(::Type{T} = Float64) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, 60, -60, Inf64],
            # OptimizeRadius = [false, true, true, false],
            Thickness = [Inf64, 10.0, 57.8, missing],
            # OptimizeThickness = [false, false, false, false],
            Conic = [missing, 0.01, 0.01, missing],
            # OptimizeConic = [false, true, true, false],
            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        )
    )
end

function doubleconvexlensonly(frontradius::T,rearradius::T) where{T<:Real}
    AxisymmetricLens{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [convert(T, Inf64), frontradius, rearradius, convert(T, Inf64)],
            # OptimizeRadius = [false, true, true, false],
            Thickness = [convert(T, Inf64), convert(T, 10.0), convert(T, 57.8), missing],
            # OptimizeThickness = [false, false, false, false],
            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [convert(T, Inf64), convert(T, 9.0), convert(T, 9.0), convert(T, 15.0)]
        )
    )
end

function doubleconvex(
    ::Type{T} = Float64;
    temperature::Unitful.Temperature = GlassCat.TEMP_REF_UNITFUL,
    pressure::T = convert(T, PRESSURE_REF)
) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, 60, -60, Inf64],
            # OptimizeRadius = [false, true, true, false],
            Thickness = [Inf64, 10.0, 57.8, missing],
            # OptimizeThickness = [false, true, true, false],
            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        );
        temperature,
        pressure
    )
end

function doubleconcave(::Type{T} = Float64) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, -41.0, 41.0, Inf64],
            Thickness = [Inf64, 10.0, 57.8, missing],
            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        )
    )
end

function planoconcaverefl(::Type{T} = Float64) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, Inf64, -41.0, Inf64],
            Thickness = [Inf64, 10.0, -57.8, missing],
            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 25.0],
            Reflectance = [missing, missing, 1.0, missing]
        )
    )
end

function concaveplano(::Type{T} = Float64) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", "Standard", "Standard", "Image"],
            Radius = [Inf64, -41.0, Inf64, Inf64],
            Thickness = [Inf64, 10.0, 57.8, missing],
            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        )
    )
end

function planoplano(::Type{T} = Float64) where {T<:Real}
    AxisymmetricOpticalSystem{T}(
        DataFrame(
            SurfaceType = ["Object", 1, 2, "Image"],
            Radius = [Inf64, Inf64, Inf64, Inf64],
            Thickness = [Inf64, 10.0, 57.8, missing],
            Material = [Air, Examples_N_BK7, Air, missing],
            SemiDiameter = [Inf64, 9.0, 9.0, 15.0]
        )
    )
end

function autodrawrays(lens::AxisymmetricOpticalSystem = cooketriplet(), angle = 10; kwargs...)
    f1 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourcenum = 1)
    Vis.drawtracerays(lens, raygenerator = f1, test = true, trackallrays = true, colorbysourcenum = true; kwargs...)
    f2 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourceangle = angle / 180 * π, sourcenum = 2)
    Vis.drawtracerays!(lens, raygenerator = f2, test = true, trackallrays = true, colorbysourcenum = true; kwargs...)
end

function autospotdiag(lens::AxisymmetricOpticalSystem = cooketriplet(); kwargs...)
    f1 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourcenum = 1)
    f2 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourceangle = 5 / 180 * π, sourcenum = 2)
    f3 = HexapolarField(lens, collimated = true, wavelength = 0.45, sourceangle = 10 / 180 * π, sourcenum = 3)
    Vis.spotdiaggrid(lens, [f1, f2, f3]; kwargs...)
end

# Display the spot diagram of a simple cooketriplet lens
function hexapolarspotdiagramexample(lens = cooketriplet(), numrings::Int = 5, angle = 0.0)
    Vis.spotdiag(lens, samples = numrings, sourceangle = angle)
end

function cartesiangridspotdiagramexample(lens = cooketriplet(), numsamples::Int = 5, angle = 0.0)
    Vis.spotdiag(lens, hexapolar = false, samples = numsamples, sourceangle = angle)
end

function prism_refraction()
    # build the triangular prism
    int = FresnelInterface{Float64}(Examples_N_SF14, Air)
    s = 2.0
    prism = (
        Plane(
            SVector(0.0, -1.0, 0.0),
            SVector(0.0, -s, 0.0),
            interface = int,
            vishalfsizeu = 2 * s,
            vishalfsizev = 2 * s
        ) ∩
        Plane(
            SVector(0.0, sind(30), cosd(30)),
            SVector(0.0, s * sind(30), s * cosd(30)),
            interface = int,
            vishalfsizeu = 2 * s,
            vishalfsizev = 2 * s
        ) ∩
        Plane(
            SVector(0.0, sind(30), -cosd(30)),
            SVector(0.0, s * sind(30), -s * cosd(30)),
            interface = int,
            vishalfsizeu = 2 * s,
            vishalfsizev = 2 * s
        )
    )
    sys = CSGOpticalSystem(LensAssembly(prism()), Rectangle(15.0, 15.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -20.0), interface = opaqueinterface()))
    # create some 'white' light
    rays = Vector{OpticalRay{Float64,3}}(undef, 0)
    for i in 0:7
        λ = ((i / 7) * 200 + 450) / 1000
        r = OpticalRay(SVector(0.0, -3.0, 10.0), SVector(0.0, 0.5, -1.0), 1.0, λ)
        push!(rays, r)
    end
    raygen = RayListSource(rays)
    # draw the result
    Vis.drawtracerays(sys, raygenerator = raygen, test = true, trackallrays = true)
end

function fresnel(convex = true; kwargs...)
    lens = FresnelLens(Examples_N_BK7, 0.0, convex ? 15.0 : -15.0, 1.0, 8.0, 0.8, conic = 0.1)
    sys = CSGOpticalSystem(LensAssembly(lens()), Rectangle(15.0, 15.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, -25.0), interface = opaqueinterface()))
    Vis.drawtracerays(sys; test = true, trackallrays = true, numdivisions = 30, kwargs...)
end

function eyetrackHOE(nrays = 5000, det = false, showhead = true, zeroorder = false; kwargs...)
    # TODO update for new specs from Chris
    hoehalfwidth = 50.0 #25.0
    hoehalfheight = 50.0 #22.5
    hoecenter = SVector(-8.0 - 25.0, 0.0, -10.0 - 25.0)
    rect = Rectangle(hoehalfheight, hoehalfwidth, SVector(0.0, 1.0, 0.0), hoecenter)
    er = 15.0
    cornea_rad = 7.85
    corneavertex = SVector(0.0, er, 0.0)
    sourceloc = SVector(-33.0, er, 0.0)
    camloc = SVector(20.0, 3.0, -11.0)
    camdir = corneavertex - camloc
    camdir_norm = normalize(camdir)

    interfaces = []

    # offset = SVector(-5.0, 10.0, -10.0)
    # for θ in 0:(π / 6):(2π)
    #     ledloc = SVector(20 * cos(θ) + offset[1], 0 + offset[2], 15 * sin(θ) + offset[3])
    #     int = HologramInterface(ledloc, ConvergingBeam, sourceloc, DivergingBeam, 0.78, 100.0, Air, Examples_N_BK7, Air, Air, Air, 0.05, zeroorder)
    #     push!(interfaces, int)
    # end

    dirs = [SVector(0.7713, 0.6350, -0.0437), SVector(0.5667, 0.8111, -0.1445), SVector(0.3400, 0.9349, -0.1017), SVector(0.1492, 0.9878, 0.0445), SVector(0.0249, 0.9686, 0.2474), SVector(-0.0184, 0.8855, 0.4643), SVector(0.0254, 0.7537, 0.6567), SVector(0.1548, 0.5964, 0.7876), SVector(0.3570, 0.4462, 0.8207), SVector(0.5959, 0.3470, 0.7242), SVector(0.7976, 0.3449, 0.4948), SVector(0.8680, 0.4555, 0.1978)]

    for d in dirs
        int = HologramInterface(normalize(d), CollimatedBeam, sourceloc, DivergingBeam, 0.78, 100.0, Air, Examples_N_BK7, Air, Air, Air, 0.05, zeroorder)
        # int = HologramInterface(corneavertex - 10 * d, ConvergingBeam, sourceloc, DivergingBeam, 0.78, 100.0, Air, Examples_N_BK7, Air, Air, Air, 0.05, zeroorder)
        push!(interfaces, int)
    end

    mint = MultiHologramInterface(interfaces...)
    obj = MultiHologramSurface(rect, mint)
    cornea = leaf(Sphere(cornea_rad, interface = FresnelInterface{Float64}(EYE.CORNEA, Air, reflectance = 1.0, transmission = 0.0)), translation(0.0, er + cornea_rad, 0.0))()

    # cam settings
    fnum = 2.0
    fov = 80
    sensorrad = 1.0
    barrellength = sensorrad / tand(fov / 2)
    aprad = barrellength / fnum / 2
    camrad = max(sensorrad, aprad)

    camap = Annulus(aprad, camrad, camdir_norm, camloc)
    distfromcamtoeye = norm(camdir)
    focallength = 1 / (1 / distfromcamtoeye + 1 / barrellength)
    camlens = ParaxialLensEllipse(focallength, aprad, aprad, -camdir_norm, camloc)
    barrelloc = camloc - barrellength / 2 * camdir_norm
    barreltop = Plane(camdir_norm, camloc)
    barrelbot = Plane(-camdir_norm, camloc - 3 * barrellength * camdir_norm)
    barrelrot = OpticSim.rotmatbetween(SVector(0.0, 0.0, 1.0), camdir_norm)
    cambarrel = (
        barrelbot ∩
        barreltop ∩
        leaf(Cylinder(camrad, barrellength, interface = opaqueinterface(Float64)), Transform(barrelrot, barrelloc))
    )()
    camdet = Circle(sensorrad, camdir_norm, camloc - barrellength * camdir_norm, interface = opaqueinterface(Float64))

    # sourceleft = hoecenter[1] + hoehalfwidth - sourceloc[1]
    # sourceright = hoecenter[1] - hoehalfwidth - sourceloc[1]
    # sourceleftθ = atan(sourceleft, sourceloc[2])
    # sourcerightθ = atan(sourceright, sourceloc[2])
    # midθ = (sourceleftθ + sourcerightθ) / 2
    # sourcedir = normalize(SVector(er * tan(midθ), -er, 0.0))
    # sourceextentθ = abs(midθ - sourcerightθ)
    # source = CosineOpticalSource(RandomSource(OriginPoint{Float64}(1, position = sourceloc, direction = sourcedir), nrays, sourceextentθ), 1.0, 0.78)

    rays = Vector{OpticalRay{Float64,3}}(undef, nrays)
    @simd for i in 1:nrays
        p = point(rect, rand() * 2 - 1, rand() * 2 - 1)
        rays[i] = OpticalRay(sourceloc, p - sourceloc, 1.0, 0.78)
    end
    source = RayListSource(rays)

    sys = CSGOpticalSystem(LensAssembly(obj, cornea, camlens, cambarrel, camap), camdet, 800, 800)
    if det
        Vis.show(OpticSim.traceMT(sys, source))
    else
        Vis.drawtracerays(sys; raygenerator = source, trackallrays = true, kwargs...)
        # for θ in 0:(π / 6):(2π)
        #     ledloc = SVector(20 * cos(θ) + offset[1], 0 + offset[2], 15 * sin(θ) + offset[3])
        #     Vis.draw!(leaf(Sphere(1.0), translation(ledloc...)), color = :red)
        # end
        for d in dirs
            # Vis.draw!(leaf(Sphere(1.0), translation((corneavertex - 10 * d)...)), color = :red)
            Vis.draw!((corneavertex - 50 * d, corneavertex), color = :red)
        end
        if showhead
            Vis.draw!(joinpath(@__DIR__, "../../OBJ/glasses.obj"), scale = 100.0, transform = Transform(OpticSim.rotmatd(90, 0, 0), [27.0, 45.0, -8.0]), color = :black)
            Vis.draw!(joinpath(@__DIR__, "../../OBJ/femalehead.obj"), scale = 13.0, transform = Transform(OpticSim.rotmatd(0, 0, 180), [27.0, 105.0, -148.0]), color = :white)
        end
        Vis.display()
    end
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, "convert_ames"))
Pkg.instantiate()

using DataFrames, CSV, MLJBase, CategoricalArrays

df = CSV.read(joinpath(@__DIR__, "reduced_ames.csv"))

schema(df)

price  = df.target
quality = df.OverallQual
area1 = map(df.GrLivArea) do a round(Int, a) end
area2 = map(df.x1stFlrSF) do a round(Int, a) end
area3 = map(df.TotalBsmtSF) do a round(Int, a) end
area4 = map(df.BsmtFinSF1) do a round(Int, a) end
area5 = map(df.GarageArea) do a round(Int, a) end
lot_area = map(df.LotArea) do a round(Int, a) end
garage_cars = map(df.GarageCars) do a round(Int, a) end
suburb = df.Neighborhood
council_code = map(df.MSSubClass) do a parse(Int, a[2:end]) end
year_built = map(df.YearBuilt) do a round(Int, a) end
year_upgraded =  map(df.YearRemodAdd) do a round(Int, a) end
zone = df.MSSubClass

df2 = DataFrame(price=price,
                area1=area1,
                area2=area2,
                area3=area3,
                area4=area4,
                area5=area5,
                lot_area=lot_area,
                year_built=year_built,
                year_upgraded=year_upgraded,
                quality=quality,
                garage_cars=garage_cars,
                suburb=suburb,
                council_code=council_code,
                zone=zone)

CSV.write(joinpath(@__DIR__, "ames.csv"), df)


text:
# function for parameters

"""
    check_parameters(par::PTParameters)

Check errors in `par`.
Return `nothing` if there is no error.
"""
function check_parameters(par::PTParameters, warn::Bool=false)
   if par.var_p<0
      error("negative phenotypic variance")
   end
   if par.var_g_sim<0
      error("negative QMSim genetic variance")
   end
   if par.h2_poly<0 || par.h2_poly>1 || par.h2_qtl<0 || par.h2_qtl>1 || (par.h2_poly+par.h2_qtl)<0 || (par.h2_poly+par.h2_qtl)>1
      error("heritability out of range")
   end
   if par.rep>0 && (par.rep < par.h2_poly || par.rep < par.h2_qtl || par.rep < (par.h2_poly + par.h2_qtl))
      error("repeatability>0 but smaller than heritability")
   end
   if par.rep < 0 || par.rep > 1
      error("repeatability out of range")
   end

   if warn && get_var_pe(par)<=0.0
      @warn "PE variance is set to 0 (rep==0 or rep==h2). Be aware that the repeatability model would not be appropriate if you have repeated records."
   end
end

function get_var_poly(par::PTParameters)
   return par.h2_poly * par.var_p
end

function get_var_qtl(par::PTParameters)
   return par.h2_qtl * par.var_p
end

function get_var_g_sim(par::PTParameters)
   return par.var_g_sim
end

function get_var_pe(par::PTParameters)
   var_u = get_var_poly(par) + get_var_qtl(par)
   return max((par.rep * par.var_p)-var_u, 0.0)
end

function get_var_error(par::PTParameters)
   return par.var_p - (get_var_poly(par) + get_var_qtl(par) + get_var_pe(par))
end


text:
# Kernel Functions =========================================================================

abstract type Kernel{T<:AbstractFloat} end

function string(κ::Kernel{T}) where {T}
    args = [string(getfield(κ,θ)) for θ in fieldnames(typeof(κ))]
    kernelname = typeof(κ).name.name
    string(kernelname, "{", string(T), "}(", join(args, ","), ")")
end

function show(io::IO, κ::Kernel)
    print(io, string(κ))
end

@inline eltype(::Type{<:Kernel{E}}) where {E} = E
@inline eltype(κ::Kernel) = eltype(typeof(κ))

"""
    ismercer(κ::Kernel)

Returns `true` if kernel `κ` is a Mercer kernel; `false` otherwise.
"""
ismercer(::Kernel) = false

"""
    isnegdef(κ::Kernel)

Returns `true` if the kernel `κ` is a negative definite kernel; `false` otherwise.
"""
isnegdef(::Kernel) = false

"""
    isstationary(κ::Kernel)

Returns `true` if the kernel `κ` is a stationary kernel; `false` otherwise.
"""
isstationary(κ::Kernel) = isstationary(basefunction(κ))

"""
    isisotropic(κ::Kernel)

Returns `true` if the kernel `κ` is an isotropic kernel; `false` otherwise.
"""
isisotropic(κ::Kernel)  = isisotropic(basefunction(κ))

# Mercer Kernels ===========================================================================

abstract type MercerKernel{T<:AbstractFloat} <: Kernel{T} end

@inline ismercer(::MercerKernel) = true

const mercer_kernels = [
    "exponential",
    "exponentiated",
    "rationalquadratic",
    "matern",
    "polynomial"
]

for kname in mercer_kernels
    include(joinpath("kernelfunctions", "mercer", "$(kname).jl"))
end

# Negative Definite Kernels ================================================================

abstract type NegativeDefiniteKernel{T<:AbstractFloat} <: Kernel{T} end

@inline isnegdef(::NegativeDefiniteKernel) = true

const negdef_kernels = [
    "log",
    "power"
]

for kname in negdef_kernels
    include(joinpath("kernelfunctions", "negativedefinite", "$(kname).jl"))
end

# Other Kernels ============================================================================

const other_kernels = [
    "sigmoid"
]

for kname in other_kernels
    include(joinpath("kernelfunctions", "$(kname).jl"))
end


text:
function valuableyesno(d, tg, av, yn)
    @unpack K = d
    K[av["requestid"]]["valuableyesno"] = yn
    if yn == 1
        # sendInvoice(chat_id = av["id"],
        #             title = "Talka",
        #             description = "Prisidėk prie platformos kūrimo!",
        #             payload = "payloadas",
        #             provider_token = ENV["TALKON_PROVIDER_TOKEN"],
        #             currency = "EUR",
        #             prices = [Dict("label"=>"kaina","amount"=>100)],
        #     suggested_tip_amounts = [100,200,300], max_tip_amount=500)
        msg = "Ми раді! | Džiaugiamės!"
        sendMessage(tg, chat_id = av["id"], text = msg)
    end
    tbegin(tg, av)
end


text:
import Pkg
Pkg.activate(".")
using ElasticSearch
using Test, JSON

@testset "ElasticSearch.jl" begin
      es = ES(user = "elastic", password = "elastic")
      @test es.url == "http://elastic:elastic@localhost:9200"
      @test [["a", "b"] |> SourceBlock] |> Query |> JSON.json ==
            """{"_source":["a","b"]}"""
      @test [[SortBody("_score")] |> SortBlock] |> Query |> JSON.json ==
            """{"sort":[{"_score":{"order":"asc"}}]}"""

      @test (TermsBody("name", "foo") |> json) ==
            (TermsBody("name", ["foo"]) |> json) ==
            """{"terms":{"name":["foo"]}}"""

      @test TermsBody("name", "foo") |> MustQueryBlock|>json ==
      """{"bool":{"must":[{"terms":{"name":["foo"]}}]}}"""

      @test [
            SizeBlock(10),
            FromBlock(0),
            RangeBody("correlation", Dict("gte"=>0.5)) |> MustQueryBlock,
            SourceBlock(["id"]),
            SortBlock(["_score"])
      ] |> Query |>json  === "{\"sort\":[\"_score\"],\"size\":10,\"from\":0,\"query\":{\"bool\":{\"must\":[{\"range\":{\"correlation\":{\"gte\":0.5}}}]}},\"_source\":[\"id\"]}"
end


text:
module ComprModeType
@enum ComprModeTypeEnum::UInt8 none=0 zlib=1 zstd=2
const default = zstd
end

mutable struct SEALHeader
  magic::UInt16
  header_size::UInt8
  version_major::UInt8
  version_minor::UInt8
  compr_mode::UInt8
  reserved::UInt16
  size::UInt64
end

SEALHeader() = SEALHeader(0, 0, 0, 0, 0, 0, 0)

function load_header!(header::SEALHeader, buffer::DenseVector{UInt8})
  io = IOBuffer(buffer)
  header.magic = read(io, UInt16)
  header.header_size = read(io, UInt8)
  header.version_major = read(io, UInt8)
  header.version_minor = read(io, UInt8)
  header.compr_mode = read(io, UInt8)
  header.reserved = read(io, UInt16)
  header.size = read(io, UInt64)
  return header
end


text:
# Script executed by Biowulf to run code

@everywhere include("/home/carsonc/StochasticGene/src/StochasticGene.jl")

include("/home/carsonc/StochasticGene/runfiles/scriptfunctions.jl")

@time fit_rna(parse(Int,ARGS[13]),ARGS[1],ARGS[14],ARGS[2],parse(Int,ARGS[3]),parse(Float64,ARGS[4]),ARGS[5],ARGS[6],ARGS[7],ARGS[8],ARGS[9],parse(Int,ARGS[10]),parse(Bool,ARGS[11]),parse(Bool,ARGS[12]))

# Arguments
# 1: gene
# 2: cond
# 3: G
# 4: maxtime
# 5: infolder
# 6: resultfolder
# 7: datafolder
# 8: inlabel
# 9: label
# 10: nsets (number of rate sets)
# 11: runcycle (bool)
# 12: transient (bool)
# 13: nchains
# 14: fittedparams (string), e.g. "1-2-4-7"

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:
##############################
# Probabilistic model

## Model parameters to fit
####################

"""
$(TYPEDEF)
$(FIELDS)

Type to hold information about the fitting parameters. Should be
flexible to easily change what parameters to fit.

A sample of the parameters is just a vector but can be interpreted
together with this type.  This is what is done in `make_fields1d`,
`make_loglikelihood` and `make_fwdm_fn`.

Contains:
- th0: initial parameter values
  - for the 1D fields btilde, fsl and temp, th0 corresponds to the
    deviation from the prior expectation, which is in gb.fsl, etc.
  - for other parameters it is just the straight value.
- names: name of parameters with their indices into th0 (to allow
  groups of parameters with the same name)
- th_ele: elevation of parameters. If just one parameter, then it is
  constant, otherwise linearly interpolated to get values at all
  elevations.
- gb::Bands: reference to the glacier Bands which this corresponds to

Convenience constructors (all slow)
- `Theta0(names, counts, gb, pp, pm, pn[, th])` constructs a theta0 with given names.
   Defaults th to zero-vectors.
- `Theta0(dict, gb, pp, pm, pn)` dict maps names to para-values
"""
struct Theta0{T}
    th0::Vector{T}  # the initial parameters
    names::OrderedDict{Symbol,UnitRange{Int64}} # name and index into th
    th_ele::OrderedDict{Symbol,LinSpace{Float64}} # elevation of parameters
                                           # if they are spatially
                                           # distributed (in 1D).

    ## Work arrays: these are updated
    # 1d fields defined on elevation bands gb.bands.  These are fed to
    # the forward model.
    gb::Bands
    pp::Phys
    pm::MPara
    pn::Num
    function Theta0{T}(th0, names, th_ele, gb, pp, pm, pn) where T
        @assert length(unique(keys(names)))==length(names)
        @assert length(th0)==mapreduce(length, +, 0, values(names))
        new{T}(th0, names, th_ele, gb, pp, pm, pn)
    end
end
# this constructor is slow, but the most convenient constructor.
function Theta0(names, counts, gb::Bands, pp::Phys, pm::MPara, pn::Num, th0::AbstractVector{T}=zeros(F,sum(count))) where T
    d = OrderedDict{Symbol,UnitRange{Int64}}()
    ii = 1
    for (name,count) in zip(names, counts)
        d[Symbol(name)] = ii:ii+count-1
        ii = ii+count
    end
    # make ranges of 1D fields thetas:
    th_ele = OrderedDict{Symbol,LinSpace{Float64}}()
    for name in names
        inds = d[name]
        li = length(inds)
        if li>1
            th_ele[name] = linspace(first(gb.ele), last(gb.ele), li)
        elseif li==1
            mid = (first(gb.ele) + last(gb.ele))/2
            th_ele[name] = linspace(mid,mid,1)
        else
            error("This should not happen...")
        end
    end

    return Theta0{T}(th0, d, th_ele, gb, pp, pm, pn)
end
function Theta0(d::Associative, gb::Bands, pp, pm, pn)
    names = collect(keys(d))
    counts = [length(v) for v in values(d)]
    th0 = F[]
    for vv in values(d)
        append!(th0, vv)
    end
    return Theta0(names, counts, gb, pp, pm, pn, th0)
end
Base.show(io::IO, theta0::Theta0) = println(io, "thetas with n=$(length(theta0.th0)) and keys:\n $((keys(theta0.names)...,))")
function get_varnames(theta0::Theta0; pad=true, addnumbers=true)
    varnames = String[]
    for name in keys(theta0.names)
        if length(theta0[name])==1
            push!(varnames, "$(string(name))")
        else
            for i=1:length(theta0[name])
                push!(varnames, addnumbers ? "$(string(name))$i" : "$(string(name))")
            end
        end
    end
    # pad with spaces:
    if pad
        maxlen = mapreduce(length, max, 0, varnames)
        return map(s->rpad(s, maxlen, " "), varnames)
    else
        return varnames
    end
end
"""
Get elevation for index i
"""
function get_ele(t0::Theta0, i::Int)
    for (n,inds) in t0.names
        if i in inds
            ele = t0.th_ele[n][findfirst(i.==inds)]
            return ele
        end
    end
    error("Index $i not in $t0")
end
# Array interface to directly index into thetas:
Base.IndexStyle(::Type{Theta0}) = IndexLinear()
Base.size(a::Theta0) = size(a.th0)
# Base.getindex(a::Theta0, i::Int) = a.th0[i]
# Base.setindex!(a::Theta0, v, i::Int) = a.th0[i]=v
Base.getindex(a::Theta0, s::Symbol) = @view a.th0[a.names[s]]

"""
Enumeration for the different types of model fitting

- h : fit to thickness only
- iv : fit to IV only
- h_iv : fit to h & IV
- length : fit to glacier length only (also included in h, iv and h_iv)
           Note: using this only necessitates the 1D model evaluation which is about
           100x faster than the 2d model.
- prior : fit to prior only
"""
baremodule FitTarget
using Base
@enum _FitTarget h iv h_iv length prior
end

#####
# Forward model wrapping
#####
"""
    make_fwdm_fn(theta0::Theta0, pm, pp, pn)

Create a function which runs the forward model for a theta as
argument.  One for the 2D model and one for the 1D, the latter only
runs the 1D part and does not produce a blob (i.e. ice thickness and
velocity maps).
"""
function make_fwdm_fn(theta0::Theta0, pp, pm, pn)
    gb = theta0.gb
    gl = gb.gl
    pn = Num(pn, verbose=false)
    pm = deepcopy(pm) # this gets updated in-place

    # check that we know what to set
    for (name,_) in theta0.names
        if !(name in [:btilde, :fsl, :temp]) && !isdefined(pm,name)
            error("Not implemented setting non-MPara parameters. Offending name: $name")
        end
    end

    # pre-calculate an update structure for pm
    args_to_MPara = []
    inds = Tuple{Int,Int}[]
    for (ii,name) in enumerate(fieldnames(pm))
        if haskey(theta0.names,name)
            its = theta0.names[name]
            push!(inds, (ii, its[1]))
            push!(args_to_MPara, nothing)
        else
            push!(args_to_MPara, getfield(pm,name))
        end
    end

    # https://github.com/JuliaLang/julia/issues/15276
    fwd = let inds = inds, args_to_MPara = args_to_MPara, pm = pm,
        theta0 = theta0, gb = gb, pp = pp, pn = pn
        function (theta::Vector)
            # modify pm to allow fitting of its parameters
            for (ias, its) in inds
                args_to_MPara[ias] = theta[its]
            end
            pm_ = MPara(args_to_MPara...)::MPara

            # calculate 1d fields (TODO: maybe use some work-arrays here)
            btilde, temp, fsl = make_fields1d(theta, theta0)

            # TODO: update other parameters?

            # run forward model:
            fwdsol = fwdm(gb, pp, pm_, pn, btilde, fsl, temp)
            return fwdsol
        end
    end

    fwd1d = let inds = inds, args_to_MPara = args_to_MPara, pm = pm,
        theta0 = theta0, gb = gb, pp = pp, pn = pn
        function (theta::Vector)
            # modify pm to allow fitting of its parameters
            for (ias, its) in inds
                args_to_MPara[ias] = theta[its]
            end
            pm_ = MPara(args_to_MPara...)::MPara

            # calculate 1d fields (TODO: maybe use some work-arrays here)
            btilde, temp, fsl = make_fields1d(theta, theta0)

            # TODO: update other parameters?

            # run 1D forward model:
            hs_n, tau_local, tau_mean_n, iv_n, qtot_n = fwdm1d(gb, btilde, fsl, temp, pp, pm_, pn)
            return hs_n, tau_local, tau_mean_n, iv_n, qtot_n, pm_
        end
    end
    return fwd, fwd1d
end

"""
    btilde, fsl, temp = make_fields1d(theta, theta0, names=[:btilde, :temp, :fsl])

Make the 1D fields according to the values in theta0.gb (their
expectation) and theta (their error/correction).

Note that the elevation decreases with parameter index.

TODO: make this more general to allow others to be 1D, say `dist_exp`.
"""
function make_fields1d(theta, theta0::Theta0, names=[:btilde, :temp, :fsl])
    gb = theta0.gb
    out = typeof(gb.bdot)[]
    for name in names
        if haskey(theta0.names, name)
            inds = theta0.names[name]
            li = length(inds)
            ele_theta0 = theta0.th_ele[name]
            f = piecewiselinear(ele_theta0, theta[inds]) # TODO use something other than linear interpolation?
            mean_val = getfield1d(gb, name) # this is the expectation value of the field
            push!(out, f.(gb.ele) .+ mean_val) # add error to expectation
        else # not fitting this parameter
            push!(out, getfield1d(gb, name)) # just return (expectation) value of the field
        end
    end
    return out
end
make_fields1d(theta, theta0::Theta0, name::Symbol) = make_fields1d(theta, theta0, [name])[1]

##################
## Probabilistic model
##################

"""
    get_sigmas(gl::Glacier, pm::MPara) -> sigma_h, sigma_iv

Return total sigmas, i.e. both observational (in gl.iv of gl.h) and
model-sigma (in pm.sigma_*).
"""
function get_sigmas(gl::Glacier, pm::MPara)
    @unpack sigma_h_model, sigma_iv_model = pm
    sigma_h_obs, sigma_iv_obs = gl.h.sigma, gl.iv.sigma
    sigma_h = sqrt(sigma_h_model^2 + sigma_h_obs^2)
    sigma_iv = sqrt(sigma_iv_model^2 + sigma_iv_obs^2)

    return sigma_h, sigma_iv
    # inf = convert(typeof(sigma_iv), Inf)

    # return if fit_target==FitTarget.h_iv
    #     sigma_h, sigma_iv
    # elseif fit_target==FitTarget.h
    #     sigma_h, inf
    # elseif fit_target==FitTarget.iv
    #     inf, sigma_iv
    # elseif fit_target==FitTarget.prior || fit_target==FitTarget.length
    #     inf, inf
    # else
    #     error("Unknown `fit_target`==$fit_target")
    # end
end

############
# Posterior
############
"""
$(SIGNATURES)

Take loglikelihood and logprior and return a function evaluating the
log-posterior and a blob (additional information produced during a
forward model run): `logposterior(theta)`

Note: if the log-prior is -Inf then then log-likelihood is not
evaluated and a empty blob is returned.
"""
function make_logposterior(loglikelihood, loglikelihood1d, logprior)
    f = function (theta)
        lp = logprior(theta)
        # return now if prior==-Inf, this way the forward model is not run
        # with invalid parameters.
        if lp==-Inf
            # return empty blob but needs to be type-stable:
            return -Inf, (NaN, NaN, Array{F}(0), Array{F}(0), Array{F}(0,0), Array{F}(0,0))
        end
        ll, blob = loglikelihood(theta)
        return  ll + lp, blob
    end
    f1d = function (theta)
        lp = logprior(theta)
        # return now if prior==-Inf, this way the forward model is not run
        # with invalid parameters.
        if lp==-Inf
            return -Inf
        end
        ll = loglikelihood1d(theta)
        return  ll + lp
    end
    return f, f1d
end

###
# Likelihood
###

"""
$(SIGNATURES)

Constructs a log-likelihood function: `loglikelihood(theta)` and also
returns the effective sigma_h and sigma_iv

Kwargs:
- iv_comp: how to do the comparison to IV.  Defaults to `:at_bands` for
           2D IV data, and `:at_measurement_points` otherwise.

Does this by making a closure over theta0 (containing gb), pm, pp,
pn.  Do this in a function to avoid global variables.

NOTE: this runs into https://github.com/JuliaLang/julia/issues/15276.
However on Julia 0.6 this particular incarnation of 15276 seems to be
fixed.
"""
function make_loglikelihood(theta0::Theta0, pp::Phys, pm::MPara, pn::Num,
                            fit_target::FitTarget._FitTarget;
                            iv_comp = if theta0.gb.gl.iv.vals isa VAWTools.Gridded
                            :at_bands # this is what is usually used
                            elseif  theta0.gb.gl.iv.vals isa VAWTools.Traj
                            :at_measurement_points
                            elseif theta0.gb.gl.iv.vals isa Void
                            :no_iv
                            else
                            error("Not anticipated")
                            end
                            )
    gb = theta0.gb
    gl = gb.gl
    pn = Num(pn, verbose=false)
    pm = deepcopy(pm) # use a copy as this gets updated in-place

    # # some tests
    # !hasdata(gl.iv) && sigma_iv!=0 &&
    #     warn("$gl has no IV data but sigma_iv is non-zero.")
    # !hasdata(gl.h) && sigma_h!=0 &&
    #     warn("$gl has no thickness data but sigma_h is non-zero.")

    #  forward model function
    fwd_fn, fwd1d_fn = make_fwdm_fn(theta0, pp, pm, pn)

    # https://github.com/JuliaLang/julia/issues/15276
    # "Performance of captured variables in closures"
    f = let gl=gl, gb=gb, fwd_fn=fwd_fn, in_glacier_=in_glacier(gl), fit_target=fit_target, iv_comp=iv_comp,
            gl_h=gl.h, gl_iv=gl.iv  # these two are needed because those fields are not strongly-typed in Glacier

        function (theta::Vector)
            # run forward model:
            fwdsol = fwd_fn(theta)
            pm_ = fwdsol.pm # something might change in pm due to theta
            @unpack vol, vol_above, hs1d, ivs1d, hs2d, ivs2d, ivs2d_band, ivs2d_band_mask = fwdsol
            # make std-dev of h and iv errors:
            sigma_h, sigma_iv = get_sigmas(gl, pm_)

            blob = vol, vol_above, hs1d, ivs1d, hs2d, ivs2d
            # return log-likelihood and blob
            if iv_comp==:at_bands
                return _loglikelihood(gl.dem, gl.iscalving, in_glacier_, gl_h, gl_iv, gl.ivmask, hs2d,
                                      ivs2d_band, ivs2d_band_mask, sigma_h, sigma_iv, gb.x, hs1d, fit_target,
                                      pm_.error_model, pm_.nthin_h, pm_.nthin_iv), blob
            else
                return _loglikelihood(gl.dem, gl.iscalving, in_glacier_, gl_h, gl_iv, gl.ivmask, hs2d,
                                      ivs2d, gl.glaciermask, sigma_h, sigma_iv, gb.x, hs1d, fit_target,
                                      pm_.error_model, pm_.nthin_h, pm_.nthin_iv), blob
            end
        end
    end

    # 1D loglikelihood
    f1d = let gb=gb, fwd1d_fn=fwd1d_fn, fit_target=fit_target, iscalving=gl.iscalving,
              verbose=pn.verbose

        function (theta::Vector)
            # run 1D forward model:
            hs1d = fwd1d_fn(theta)[1]

            # no blob

            # likelihood
            x = gb.x
            if verbose && fit_target!=FitTarget.length
                warn("1D log-likelihood only works with FitTarget.length")
            end
            # TODO add a comparison to Vol-Area scaling?

            # @show length1d(x, hs1d), x[end]
            # @show length_logpdf(length1d(x, hs1d), x, iscalving)
            return length_logpdf(length1d(x, hs1d), x, iscalving)
        end
    end
    return f, f1d, fwd_fn, fwd1d_fn
end

"Return interpolation to test whether an arbitrary point is inside glacier"
in_glacier(gl) = Interp.interpolate((gl.dem.x, gl.dem.y), gl.glaciermask, Interp.Gridded(Interp.Constant()) )

"""
$(SIGNATURES)

Log-Likelihood for given measured-data and given forward model result.
Can be overloaded for specific glaciers.

This assumes that the errors have a normal distribution with std-dev
sigma.

Input
 - gl, ivs1d_measured, ivs1d, hs2d, ivs2d, sigma_h, sigma_iv

If sigmas==Inf then that error is not calculated.  If both sigma_iv and sigma_h
are Inf, then the log-likelihood==0.  This is equivalent of sampling
just the prior (but the blob will still be evaluated and included, i.e. this still runs
the forward model and thus is way slower than just sampling the prior).

Output:
- loglikelihood

Note the log of the product of all the independent Gaussians is:

    -n ( ln(σ) + 1/2*ln(2π) ) - 1/(2σ²) * Σ (x_i - μ_i)²
"""
function _loglikelihood(dem, iscalving, in_glacier, gl_h, gl_iv, gl_ivmask,
                        hs2d, ivs2d, ivs2d_mask, sigma_h, sigma_iv, x, hs1d, fit_target,
                        error_model,
                        nthin_h::Int, nthin_iv::Int) # bug in Julia: the two ::Int are needed for type-inference
    out = 0.0
    if fit_target==FitTarget.h || fit_target==FitTarget.h_iv
        if gl_h isa SomeData{F}
            # For some GlaThiDa glacier there is just a mean thickness value.  This uses this.
            Vmod = sum(hs2d)*step(dem.x)^2
            Vmeas = gl_h.vals
            out += -(Vmod-Vmeas)^2/(2*sigma_h^2) - log(sigma_h*sqrt(2pi))
        elseif gl_h isa SomeData{Traj{F}}
            if error_model==ErrorModel.sqrerr
                ## Normal dist
                te, ne = thick_err(dem, in_glacier, gl_h, hs2d, 1:length(gl_h.vals.splits), nthin_h, sqrerr, +)
                if ne>0
                    out += -te/(2*sigma_h^2) - ne*log(sigma_h*sqrt(2pi))
                end
            elseif error_model==ErrorModel.abserr
                ## Laplace dist (abs-errors), such that sigma is still the stdev
                te, ne = thick_err(dem, in_glacier, gl_h, hs2d, 1:length(gl_h.vals.splits), nthin_h, abserr, +)
                if ne>0
                    out += -te/(sigma_h/sqrt(2)) - ne*log(sigma_h*sqrt(2))
                end
            else
                error("ErrorModel $error_model not implemented")
            end
        else
            error("Not supported gl_h type: $(typeof(gl_h))")
        end
    end
    if fit_target==FitTarget.iv || fit_target==FitTarget.h_iv
        if error_model==ErrorModel.sqrerr
            ## Normal dist
            tiv2d, niv2d = iv2d_err(dem, in_glacier, gl_iv, gl_ivmask, ivs2d, ivs2d_mask, nthin_iv, sqrerr, +)
            if niv2d>0
                out += -tiv2d/(2*sigma_iv^2) - niv2d*log(sigma_iv*sqrt(2pi))
            end
        elseif error_model==ErrorModel.abserr
            ## Laplace dist (abs-errors)
            tiv2d, niv2d = iv2d_err(dem, in_glacier, gl_iv, gl_ivmask, ivs2d, ivs2d_mask, nthin_iv, abserr, +)
            if niv2d>0
                out += -tiv2d/(sigma_iv/sqrt(2)) - niv2d*log(sigma_iv*sqrt(2))
            end
        else
            error("ErrorModel $error_model not implemented")
        end
    end
    ## Make sure the length is correct
    # @show length1d(x, hs1d), x[end]
    # @show length_logpdf(length1d(x, hs1d), x, iscalving)
    out += length_logpdf(length1d(x, hs1d), x, iscalving)

    # TODO add a comparison to V-A scaling?

    return out
end

"""
    length_logpdf(theta::Vector, theta0::Theta0)

This encodes that a glacier should have the length of the domain.

TODO: make parameters adjustable.
"""
function length_logpdf(length, x, iscalving)
    if iscalving
        # then just make sure the glacier reaches the end
        return length<x[end] ? -Inf : 0.0
    else # land terminating, make sure glacier length is close to what it should be
        expectation = x[end]
        sigma = expectation*0.005
        percent = 0.03
        range_abs = [(1-percent) * expectation, (1+percent) * expectation]
        # using the subnormal make the probability drop much faster than Normal beyond range_abs:
        return logprior_scalar_subnormal(length, expectation, sigma, range_abs)
    end
end

# Define a function to quantify a model vs measurements error. The
# standard is to assume a normal distributed mean square error.
"Square error of one observation."
sqrerr(y_measured, y_modelled) = (y_measured-y_modelled)^2

## Note:
# other error functions could be used, as e.g. below.
# Also note the term "heteroscedasticity", which means that the error-sigma
# is dependent on the variable, e.g. relative error.

"Absolute error of one observation."
abserr(y_measured, y_modelled) = abs(y_measured-y_modelled)

"""
Relative square errors [1] using the trick of [2].

[1] https://ieeexplore.ieee.org/abstract/document/159804
[2] https://www.jstor.org/stable/27920135
"""
rel_sqrerr(y_measured, y_modelled) =
    ((y_measured-y_modelled)/y_measured)^2 +
    ((y_measured-y_modelled)/y_modelled)^2

"""
Relative square errors. Not sure this makes sense.
"""
rel_sqrerr_min(y_measured, y_modelled) = min(
    ((y_measured-y_modelled)/y_measured)^2,
    ((y_measured-y_modelled)/y_modelled)^2)

"""
Relative absolute error (LARE)
Chen et al 2010, Eq3
https://www.jstor.org/stable/27920135
"""
rel_abserr(y_measured, y_modelled) =
    abs((y_measured-y_modelled)/y_measured) +
    abs((y_measured-y_modelled)/y_modelled)

rel_abserr_min(y_measured, y_modelled) = min(
    abs((y_measured-y_modelled)/y_measured),
    abs((y_measured-y_modelled)/y_modelled))

"""
$(SIGNATURES)

Reduction (sum) of ice thickness (square) errors. Can be overloaded
for a particular glacier.

Input either:
- dem: gl.dem
- in_glacier : interpolation, true if in glacier
- gl_h: gl.h
- hs2d: 2D thickness model output
or:
- gl
- hs2d

Optional positional args:
- lines=1:length(gl_h.vals.splits) : which radar lines to process
- thin=10 : thin the radar lines.  Needed because of spatial correlation.
- errfn=sqrerr : the error function to use `(value, measurement) -> comparison`
- reduction=+ : how to reduce the errors

Return:
- reduction of errors (sum of square errors)
- number of points used

TODO
- think about weights: maybe two very closely spaced points should
  count less than two far ones?  -> there is some correlation
  between them, what's the length?
"""
function thick_err(dem, in_glacier, gl_h::SomeData{Traj{F}}, hs2d, lines=1:length(gl_h.vals.splits),
                   thin=10, errfn=sqrerr, reduction=+, neutral_element=0.0)
    ihs  = Interp.interpolate((dem.x, dem.y), hs2d, Interp.Gridded(Interp.Linear()) )
    out = neutral_element
    n = 0 # counter for number of points at which an error is evaluated
    start = 3
    step = thin # TODO hack...
    tr = gl_h.vals
    for line in lines
        for i = tr.splits[line][1]:step:tr.splits[line][end]
            x,y,h = tr.x[i],tr.y[i],tr.v[i]
            if in_glacier[x,y]==1 # only inside glacier
                err = errfn(ihs[x,y], h)
                isnan(err) && continue
                out = reduction(out, errfn(ihs[x,y], h))
                n += 1
            end
        end
    end
    return out, n
end
thick_err(gl::Glacier, hs2d::Matrix, lines=1:length(gl.h.vals.splits), thin=10, errfn=sqrerr, reduction=+, neutral_element=0.0) =
    thick_err(gl.dem, in_glacier(gl), gl.h, hs2d, lines, thin, errfn, reduction, neutral_element)
thick_err(gl::Glacier, hs2d::Gridded, lines=1:length(gl.h.vals.splits), thin=10, errfn=sqrerr, reduction=+, neutral_element=0.0) =
    thick_err(gl.dem, in_glacier(gl), gl.h, hs2d.v, lines, thin, errfn, reduction, neutral_element)

"""
$(SIGNATURES)

Sum of 2D ice velocity (square) errors. Can be overloaded for a particular
glacier.

Input:
- gl
- ivs2d: 2D velocity model output
- ivs2d_mask: where 2D velocity model output should be compared to observations
- gl_iv: IV observations
- gl_ivmask: and mask

Optional:
- thin=10: thin the velocity measurements.  Needed because of correlation.

"""
function iv2d_err(dem, in_glacier, gl_iv::SomeData{Gridded{F}}, gl_ivmask, ivs2d, ivs2d_mask,
                  thin=10, errfn=sqrerr, reduction=+)
    iv_ = gl_iv.vals
    iv_obs, flag  = size(iv_.v)==size(ivs2d) ? (iv_.v,true) :
        (Interp.interpolate((iv_.x, iv_.y), iv_.v, Interp.Gridded(Interp.Linear()) ), false)
    iv_obs_mask = flag ? gl_ivmask :
        Interp.interpolate((iv_.x, iv_.y), gl_ivmask, Interp.Gridded(Interp.Constant()) )

    out = 0.0
    n = 0
    step = thin # TODO hack...
    for iy=1:step:length(dem.y)
        y = dem.y[iy]
        for ix=1:step:length(dem.x)
            # skip points not in this mask
            !(ivs2d_mask[ix,iy]) && continue
            x = dem.x[ix]

            # check that observation-mask is true:
            if (flag ? iv_obs_mask[ix,iy] : iv_obs_mask[x,y]==1)
                obs = (flag ? iv_obs[ix,iy] : iv_obs[x,y])::Float64
                out = reduction(out, errfn(obs, ivs2d[ix,iy]))
                n += 1
            end
        end
    end
    return out, n
end
# for trajectory data
function iv2d_err(dem, in_glacier, gl_iv::SomeData{Traj{F}}, gl_ivmask, ivs2d, ivs2d_mask,
                  thin=1, errfn=sqrerr, reduction=+)
    iiv  = Interp.interpolate((dem.x, dem.y), ivs2d, Interp.Gridded(Interp.Linear()) )
    out = 0.0
    n = 0
    step = thin # TODO hack...
    for i=1:step:length(gl_iv.vals.v)
        x = gl_iv.vals.x[i]
        y = gl_iv.vals.y[i]
        iv = gl_iv.vals.v[i]
        if in_glacier[x,y]==1 # only inside glacier
            out = reduction(out, errfn(iiv[x,y], iv))
            n += 1
        end
    end
    return out, n
end
iv2d_err(dem, in_glacier, gl_iv::SomeData{Void}, gl_ivmask, ivs2d, ivs2d_mask,
         thin=10, errfn=sqrerr, reduction=+) = 0.0, 0

iv2d_err(gl, ivs2d::Matrix, ivmask=gl.ivmask, thin=10, errfn=sqrerr, reduction=+) =
    iv2d_err(gl.dem, in_glacier(gl), gl.iv, ivmask, ivs2d, thin, errfn, reduction)
iv2d_err(gl, ivs2d::Gridded, ivmask=gl.ivmask, thin=10, errfn=sqrerr, reduction=+) =
    iv2d_err(gl.dem, in_glacier(gl), gl.iv, ivmask, ivs2d.v, thin, errfn, reduction)

## Not used anymore:
# "Sum of 1D ice velocity errors.  This does not work well."
# function iv1d_err(ivs1d_measured, ivs1d, thin=3)
#     out = 0.0
#     n = 0
#     step = thin
#     for i=1:step:length(ivs1d_measured)
#         iv1 = ivs1d_measured[i]
#         iv2 = ivs1d[i]
#         if iv1>0 # mask
#             out += sqrerr(iv1,iv2)
#             n += 1
#         end
#     end
#     return out, n
# end