text:
"""
```
measurement(m::Model1002{T}, TTT::Matrix{T}, RRR::Matrix{T},
    CCC::Vector{T}; reg::Int = 1) where {T<:AbstractFloat}
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::Model1002{T},
                     TTT::Matrix{T},
                     RRR::Matrix{T},
                     CCC::Vector{T}; reg::Int = 1) where {T<:AbstractFloat}

    endo     = m.endogenous_states
    endo_new = m.endogenous_states_augmented
    exo      = m.exogenous_shocks
    obs      = m.observables

    _n_observables = n_observables(m)
    _n_states = n_states_augmented(m)
    _n_shocks_exogenous = n_shocks_exogenous(m)

    ZZ = zeros(_n_observables, _n_states)
    DD = zeros(_n_observables)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    for para in m.parameters
        if !isempty(para.regimes)
            ModelConstructors.toggle_regime!(para, reg)
        end
    end

    no_integ_inds = inds_states_no_integ_series(m)
    if haskey(m.endogenous_states, :pgap_t)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:pgap_t]])
    end
    if haskey(m.endogenous_states, :ygap_t)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:ygap_t]])
    end
    if haskey(m.endogenous_states, :rw_t)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:rw_t]])
    end
    if haskey(m.endogenous_states, :Rref_t)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:Rref_t]])
    end
    if (get_setting(m, :add_laborproductivity_measurement) || get_setting(m, :add_nominalgdp_level) ||
        get_setting(m, :add_cumulative)) || (haskey(m.endogenous_states, :pgap_t)) || (haskey(m.endogenous_states, :ygap_t)) ||
        (haskey(m.endogenous_states, :rw_t)) || (haskey(m.endogenous_states, :Rref_t))
        # Remove integrated states (e.g. states w/unit roots)
        TTT = @view TTT[no_integ_inds, no_integ_inds]
        CCC = @view CCC[no_integ_inds]
    end

    ## GDP growth - Quarterly!
    ZZ[obs[:obs_gdp], endo[:y_t]]          = 1.0
    ZZ[obs[:obs_gdp], endo_new[:y_t1]]     = -1.0
    ZZ[obs[:obs_gdp], endo[:z_t]]          = 1.0
    ZZ[obs[:obs_gdp], endo_new[:e_gdp_t]]  = 1.0
    ZZ[obs[:obs_gdp], endo_new[:e_gdp_t1]] = -m[:me_level]
    DD[obs[:obs_gdp]]                      = 100*(exp(m[:z_star])-1)

    if haskey(get_settings(m), :add_iid_cond_obs_gdp_meas_err) ?
        get_setting(m, :add_iid_cond_obs_gdp_meas_err) : false
        ZZ[obs[:obs_gdp], endo_new[:e_condgdp_t]] = 1.
    end

    ## GDI growth- Quarterly!
    ZZ[obs[:obs_gdi], endo[:y_t]]          = m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo_new[:y_t1]]     = -m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo[:z_t]]          = m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo_new[:e_gdi_t]]  = 1.0
    ZZ[obs[:obs_gdi], endo_new[:e_gdi_t1]] = -m[:me_level]
    DD[obs[:obs_gdi]]                      = 100*(exp(m[:z_star])-1) + m[:δ_gdi]

    ## Hours growth
    ZZ[obs[:obs_hours], endo[:L_t]] = 1.0
    DD[obs[:obs_hours]]             = m[:Lmean]

    ## Labor Share/real wage growth
    if subspec(m) in ["ss16", "ss17"]
        # log(labor_share) = log(wage) + log(hours) - log(GDP)
        ZZ[obs[:obs_laborshare], endo[:w_t]] = 1.
        ZZ[obs[:obs_laborshare], endo[:L_t]] = 1.
        DD[obs[:obs_laborshare]] = 100. * log(m[:wstar] * m[:Lstar] / m[:ystar])
        ZZ[obs[:obs_laborshare], endo[:y_t]] = -1.
    else
        ZZ[obs[:obs_wages], endo[:w_t]]      = 1.0
        ZZ[obs[:obs_wages], endo_new[:w_t1]] = -1.0
        ZZ[obs[:obs_wages], endo[:z_t]]      = 1.0
        DD[obs[:obs_wages]]                  = 100*(exp(m[:z_star])-1)
    end

    ## Inflation (GDP Deflator)
    ZZ[obs[:obs_gdpdeflator], endo[:π_t]]            = m[:Γ_gdpdef]
    ZZ[obs[:obs_gdpdeflator], endo_new[:e_gdpdef_t]] = 1.0
    DD[obs[:obs_gdpdeflator]]                        = 100*(m[:π_star]-1) + m[:δ_gdpdef]

    ## Inflation (Core PCE)
    ZZ[obs[:obs_corepce], endo[:π_t]]             = 1.0
    ZZ[obs[:obs_corepce], endo_new[:e_corepce_t]] = 1.0
    DD[obs[:obs_corepce]]                         = 100*(m[:π_star]-1)

    ## Nominal interest rate
    ZZ[obs[:obs_nominalrate], endo[:R_t]] = 1.0
    DD[obs[:obs_nominalrate]]             = m[:Rstarn]

    ## Consumption Growth
    ZZ[obs[:obs_consumption], endo[:c_t]]      = 1.0
    ZZ[obs[:obs_consumption], endo_new[:c_t1]] = -1.0
    ZZ[obs[:obs_consumption], endo[:z_t]]      = 1.0
    DD[obs[:obs_consumption]]                  = 100*(exp(m[:z_star])-1)

    ## Investment Growth
    ZZ[obs[:obs_investment], endo[:i_t]]      = 1.0
    ZZ[obs[:obs_investment], endo_new[:i_t1]] = -1.0
    ZZ[obs[:obs_investment], endo[:z_t]]      = 1.0
    DD[obs[:obs_investment]]                  = 100*(exp(m[:z_star])-1)

    ## Spreads
    ZZ[obs[:obs_spread], endo[:ERktil_t]] = 1.0
    ZZ[obs[:obs_spread], endo[:R_t]]       = -1.0
    DD[obs[:obs_spread]]                   = 100*log(m[:spr])

    ## 10 yrs infl exp

    TTT10                                      = (1/40) * ((Matrix{Float64}(I, size(TTT, 1), size(TTT,1))
                                                            - TTT) \ (TTT - TTT^41))
    ZZ[obs[:obs_longinflation], no_integ_inds] = TTT10[endo[:π_t], :]
    DD[obs[:obs_longinflation]]                = 100*(m[:π_star]-1)

    ## Long Rate
    ZZ[obs[:obs_longrate], no_integ_inds]     = ZZ[6, no_integ_inds]' * TTT10
    ZZ[obs[:obs_longrate], endo_new[:e_lr_t]] = 1.0
    DD[obs[:obs_longrate]]                    = m[:Rstarn]

    ## TFP
    ZZ[obs[:obs_tfp], endo[:z_t]] = (1-m[:α])*m[:Iendoα] + 1*(1-m[:Iendoα])
    if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
        ZZ[obs[:obs_tfp], endo_new[:e_tfp_t]]  = 1.0
        ZZ[obs[:obs_tfp], endo_new[:e_tfp_t1]] = -m[:me_level]
    else
        ZZ[obs[:obs_tfp], endo_new[:e_tfp_t]] = 1.0
    end
    if !(subspec(m) in ["ss15", "ss16"])
        ZZ[obs[:obs_tfp], endo[:u_t]]       = m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα] )
        ZZ[obs[:obs_tfp], endo_new[:u_t1]]  = -(m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα]) )
    end

    ## Set up structural shocks covariance matrix
    QQ[exo[:g_sh], exo[:g_sh]]             = m[:σ_g]^2
    QQ[exo[:b_sh], exo[:b_sh]]             = m[:σ_b]^2
    QQ[exo[:μ_sh], exo[:μ_sh]]             = m[:σ_μ]^2
    QQ[exo[:ztil_sh], exo[:ztil_sh]]       = m[:σ_ztil]^2
    QQ[exo[:λ_f_sh], exo[:λ_f_sh]]         = m[:σ_λ_f]^2
    QQ[exo[:λ_w_sh], exo[:λ_w_sh]]         = m[:σ_λ_w]^2
    QQ[exo[:rm_sh], exo[:rm_sh]]           = m[:σ_r_m]^2
    QQ[exo[:σ_ω_sh], exo[:σ_ω_sh]]         = m[:σ_σ_ω]^2
    QQ[exo[:μ_e_sh], exo[:μ_e_sh]]         = m[:σ_μ_e]^2
    QQ[exo[:γ_sh], exo[:γ_sh]]             = m[:σ_γ]^2
    QQ[exo[:π_star_sh], exo[:π_star_sh]]   = m[:σ_π_star]^2
    QQ[exo[:lr_sh], exo[:lr_sh]]           = m[:σ_lr]^2
    QQ[exo[:zp_sh], exo[:zp_sh]]           = m[:σ_z_p]^2
    QQ[exo[:tfp_sh], exo[:tfp_sh]]         = m[:σ_tfp]^2
    QQ[exo[:gdpdef_sh], exo[:gdpdef_sh]]   = m[:σ_gdpdef]^2
    QQ[exo[:corepce_sh], exo[:corepce_sh]] = m[:σ_corepce]^2
    QQ[exo[:gdp_sh], exo[:gdp_sh]]         = m[:σ_gdp]^2
    QQ[exo[:gdi_sh], exo[:gdi_sh]]         = m[:σ_gdi]^2

    if subspec(m) in ["ss59", "ss60", "ss61"]
        QQ[exo[:ziid_sh], exo[:ziid_sh]]   = m[:σ_ziid]^2
        QQ[exo[:biidc_sh], exo[:biidc_sh]] = m[:σ_biidc]^2
        QQ[exo[:φ_sh], exo[:φ_sh]]         = m[:σ_φ]^2
    end

    if haskey(get_settings(m), :add_iid_cond_obs_gdp_meas_err) ?
        get_setting(m, :add_iid_cond_obs_gdp_meas_err) : false
        QQ[exo[:condgdp_sh], exo[:condgdp_sh]] = m[:σ_condgdp] ^ 2
    end

    if haskey(get_settings(m), :add_iid_anticipated_obs_gdp_meas_err) ?
        get_setting(m, :add_iid_anticipated_obs_gdp_meas_err) : false
        QQ[exo[:gdpexp_sh], exo[:gdpexp_sh]] = m[:σ_gdpexp] ^ 2
    end

    # Automated addition of anticipated shocks to QQ
    for (k, v) in get_setting(m, :antshocks)
        for i = 1:v
            QQ[exo[Symbol(k, "_shl$i")], exo[Symbol(k, "_shl$i")]] = m[Symbol("σ_", k, "$i")]^2
        end
    end

    ## Anticipated observables

    # Anticipated monetary policy shocks
    for i = 1:n_mon_anticipated_shocks(m)
        ZZ[obs[Symbol("obs_nominalrate$i")], no_integ_inds] = ZZ[obs[:obs_nominalrate], no_integ_inds]' * (TTT^i)
        DD[obs[Symbol("obs_nominalrate$i")]]    = m[:Rstarn]
        if subspec(m) == "ss11"
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[:σ_r_m]^2 / n_mon_anticipated_shocks(m)
        else
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[Symbol("σ_r_m$i")]^2
        end
    end

    # Anticipated GDP growth
    if haskey(get_settings(m), :add_anticipated_obs_gdp)
        if get_setting(m, :add_anticipated_obs_gdp)
            for i = 1:get_setting(m, :n_anticipated_obs_gdp)
                ZZ_obs_gdp = ZZ[obs[:obs_gdp], :]
                meas_err = haskey(get_settings(m), :meas_err_anticipated_obs_gdp) ?
                    get_setting(m, :meas_err_anticipated_obs_gdp) : 0. # Ignore measurement error for anticipated GDP growth
                ZZ_obs_gdp[endo_new[:e_gdp_t]]  = meas_err
                ZZ_obs_gdp[endo_new[:e_gdp_t1]] = -meas_err * m[:me_level]
                ZZ[obs[Symbol("obs_gdp$i")], no_integ_inds] = ZZ_obs_gdp[no_integ_inds]' * (TTT^i)
                DD[obs[Symbol("obs_gdp$i")]]                = 100. * (exp(m[:z_star]) - 1.)
                if haskey(get_settings(m), :add_iid_anticipated_obs_gdp_meas_err) ?
                    get_setting(m, :add_iid_anticipated_obs_gdp_meas_err) : false
                    ZZ[obs[Symbol("obs_gdp$i")], endo_new[:e_gdpexp_t]] = 1.
                end
            end
        end
    end

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC .!= 0.)
        # Are we using the zero rate alternative policy or not?
        is_zero_rate_rule = (haskey(get_settings(m), :replace_eqcond_func_dict) ?
                             (!isempty(get_setting(m, :replace_eqcond_func_dict)) && get_setting(m, :replace_eqcond)) : false) ||
                             get_setting(m, :alternative_policy).eqcond == zero_rate_eqcond

        # If we are using the zero rate rule, then we don't adjust DD. If we run the code block below,
        # then the zero rate rule will cause CCC to be nonzero in multiple places, leading to the
        # wrong steady state measurement vector DD.
        #
        # The current coding of is_zero_rate_rule assumes that, when using a zero_rate rule (perhaps with gensys2 and temp alt policies)
        # any regimes in which the zero rate rule does not apply also will not have to update DD to reflect nonzero CCC.
        # This is typically the case for m1002.
        #
        # Currently it looks like we modify CCC in augment_states to track expected inflation,
        # but it doesn't look like we use this variable anywhere anyway. So ignoring this code block is fine on that front.
        # However, all the observables and pseudo-observables were not thoroughly checked, so it is possible
        # that the measurement equation is flawed in some way given how this has been coded.
        if !is_zero_rate_rule
            DD += ZZ[:, no_integ_inds]*((UniformScaling(1) - TTT) \ CCC)
        end
    end

    for para in m.parameters
        if !isempty(para.regimes)
            ModelConstructors.toggle_regime!(para, 1)
        end
    end

    return Measurement(ZZ, DD, QQ, EE)
end


text:
###################################################################################################
# FUNCTIONS
###################################################################################################
"""
    `create_Ux([N1,...,NK],[α1,...,αK],β)`

    Function that creates a function to compute potential energy for a region with [N1,...,NK],
    parameters [α1,...,αK], and correlation β.

    # Examples
    ```julia-repl
    julia> Ux_fun = CpelTdm.create_Ux([2,2],[1.0,1.0],1.0);
    julia> Ux_fun([1,1,1,1])
    -7.0
    ```
"""
function create_Ux(n::Vector{Int64},a::Vector{Float64},b::Float64)

    # Define function given n, α, and β
    function Ux(x::Vector{Int64})::Float64
        u = a[1] * sum(x[1:n[1]]) + b * sum(x[1:(end-1)] .* x[2:end])
        @inbounds for i in 2:length(n)
            u += a[i] * sum(x[(sum(n[1:(i-1)])+1):(sum(n[1:i]))])
        end
        return -u
    end

    # Return function
    return Ux

end # end create_Ux
"""
    `get_W(N,α,β)`

    Function that returns W^{N-1} computed efficiently via sum of rank-1 matrices
    assuming α and β.

    # Examples
    ```julia-repl
    julia> CpelTdm.get_W(4,0.0,0.0)
    2×2 Array{Float64,2}:
    4.0  4.0
    4.0  4.0
    ```
"""
function get_W(n::Int64,a::Float64,b::Float64)::Array{Float64,2}

    # if n=1 then return identity
    n==1 && return [1.0 0.0;0.0 1.0]

    # Compute relevant quantities
    exp_a = exp(a)
    exp_b = exp(b)
    cosh_a = 0.5*(exp_a+1.0/exp_a)
    sinh_a = exp_a-cosh_a

    # Eigenvalues
    aux1 = exp_b * cosh_a
    aux2 = sqrt(1.0 + exp_b^4*sinh_a^2)
    lambda1N = (aux1-1.0/exp_b*aux2)^(n-1)
    lambda2N = (aux1+1.0/exp_b*aux2)^(n-1)

    # Eigenvectors
    aux1 = -exp_b^2 * sinh_a
    e1 = [aux1-aux2; 1.0]
    e1 /= sqrt(e1'*e1)
    e2 = [aux1+aux2; 1.0]
    e2 /= sqrt(e2'*e2)

    # Return W^{N-1}
    return e1*e1'*lambda1N + e2*e2'*lambda2N

end
"""
    `get_V([α1,α2],β)`

    Function that returns V matrix assuming [α1,α2] and β.

    # Examples
    ```julia-repl
    julia> CpelTdm.get_V([1.0,1.0],1.0)
    2×2 Array{Float64,2}:
    1.0       0.367879
    0.367879  7.38906
    ```
"""
function get_V(a::Vector{Float64},b::Float64)::Array{Float64,2}

    # Compute auxiliary vars
    exp_b = exp(b)
    exp_a_p = exp(0.5*sum(a))
    exp_a_m = exp(0.5*(a[2]-a[1]))

    # Return V
    return [exp_b/exp_a_p exp_a_m/exp_b;1.0/(exp_b*exp_a_m) exp_b*exp_a_p]

end
"""
    `get_u(α)`

    Function that boundary vector for Z computation.

    # Examples
    ```julia-repl
    julia> CpelTdm.get_u(0.0)
    2-element Vector{Float64}:
     1.0
     1.0
    ```
"""
function get_u(a::Float64)::Vector{Float64}

    # Compute auxiliary
    exp_aux = exp(a/2.0)

    # Return u
    return [1.0/exp_aux; exp_aux]

end
"""
    `comp_Z([N1,...,NK],[α1,...,αK],β)`

    Compute partition function of a model with [N1,...,NK] CpG cites and with
    parameters [α1,...,αK] and β.

    # Examples
    ```julia-repl
    julia> CpelTdm.comp_Z([1,1,1],[1.0,1.0,1.0],1.0)
    155.37102759254836
    ```
"""
function comp_Z(n::Vector{Int64},a::Vector{Float64},b::Float64)::Float64

    # Boundary conditions.
    y = get_u(a[1])'*get_W(n[1],a[1],b)
    if length(n)>1
        y *= prod([get_V(a[(i-1):i],b)*get_W(n[i],a[i],b) for i in 2:length(n)])
    end
    y *= get_u(a[end])

    # Return Z
    return max(1e-100,y)

end
"""
    `get_∇logZ([N1,...,NK],θhat)`

    Numerically computes the gradient of the log partition function of a model with [N1,...,NK]
    CpG cites and estimated parameter vector θhat. This is equivalent to computing the expected
    value of sufficient statistis (SS).

    # Examples
    ```julia-repl
    julia> CpelTdm.get_∇logZ([1,1,1],[1.0,-1.0,1.0,0.0])
    4-element Array{Float64,1}:
     0.7615941559702503
     -0.7615941559702503
     0.7615941559335818
     -1.1600513167692226
    ```
"""
function get_∇logZ(n::Vector{Int64},θhat::Vector{Float64})::Vector{Float64}

    # Define function
    function f(θ::Vector{Float64})
        log(comp_Z(n,θ[1:(end-1)],θ[end]))
    end

    # Return ∇logZ(θ)
    return Calculus.gradient(f,θhat)

end
"""
    `check_boundary(θhat)`

    Function that returns a bool indicating whether model with parameter estimate vector θhat is on the
    boundary of the parameter space in any of its components.

    # Examples
    ```julia-repl
    julia> CpelTdm.check_boundary([1.0,1.0,1.0,1.0])
    false
    julia> CpelTdm.check_boundary([1.0,5.0,1.0,1.0])
    true
    ```
"""
function check_boundary(θhat::Vector{Float64})::Bool

    # Return true θhat on boundary.
    return any(isapprox.(abs.(θhat),ETA_MAX_ABS;atol=5e-2)) || any(abs.(θhat).>ETA_MAX_ABS)

end
"""
    `comp_g([R1,...,RK],[α1,...,αK],β,αp1,αp2)`

    Compute scaling factor in a model with [R1,...,RK] unobserved CpG cites from each block, with
    parameters [α1,...,αK] and β. The transfer-matrix method is used to obtain a computationally
    efficient expression without the need of recursive expressions. αp1 and αp2 are determined by the
    respective kind of boundary.

    # Examples
    ```julia-repl
    julia> CpelTdm.comp_g([1],[1.0],1.0,3.0,3.0)
    20.135323991555527
    ```
"""
function comp_g(r::Vector{Int64},a::Vector{Float64},b::Float64,ap1::Float64,ap2::Float64)::Float64

    # Return one if r=0
    r==0 && return 1.0

    # Boundary conditions
    y = get_u(ap1)'*get_W(r[1],a[1],b)
    if length(r)>1
        y *= prod([get_V(a[(i-1):i],b)*get_W(r[i],a[i],b) for i in 2:length(r)])
    end
    y *= get_u(ap2)

    # Return scaling factor
    return max(1e-100,y)

end
"""
    `comp_lkhd(X,[N1,...,NK],[α1,...,αK],β)`

    Compute likelihood of a partial/full observation X that can be missing values anywhere (given by
    0's).

    # Examples
    ```julia-repl
    julia> CpelTdm.comp_lkhd([1,1,0,1,1],[5],[1.0],1.0)
    0.953646032691218
    ```
"""
function comp_lkhd(x::Vector{Int64},n::Vector{Int64},a::Vector{Float64},b::Float64)::Float64

    # Avoid the rest if N=1
    length(x)==1 && return 0.5 * exp(x[1]*a[1]) / cosh(a[1])

    # Get partition function and energy function
    Z = comp_Z(n,a,b)
    Ux = create_Ux(n,a,b)

    # Find changes to/from 0 in x vector
    zerost = findall(isequal(-1),abs.(x[2:end]) - abs.(x[1:(end-1)])) .+ 1
    zeroend = findall(isequal(1),abs.(x[2:end]) - abs.(x[1:(end-1)])) .+ 1

    # Determine whether it starts/finishes with 0
    x[1]==0 && pushfirst!(zerost,1)
    x[end]==0 && push!(zeroend,length(x)+1)

    # Find subregion label for each CpG site
    subid = vcat([i*ones(Int64,n[i]) for i in 1:length(n)]...)

    # Get overall scaling factor as product of individual factors
    sf = 1.0
    @inbounds for i in 1:length(zerost)
        # Find boundary conditions (other X's or boundary of X)
        ap1 = zerost[i]==1 ? a[1] : 2.0*x[zerost[i]-1]*b+a[subid[zerost[i]]]
        ap2 = zeroend[i]==sum(n)+1 ? a[end] : 2.0*x[zeroend[i]]*b+a[subid[zeroend[i]-1]]

        # Figure out b (block IDs of missing) and r (α indices)
        b_id = subid[zerost[i]:(zeroend[i]-1)]
        n_miss = [count(x->x==n,b_id) for n in unique(b_id)]

        # Call scaling factor function
        sf *= comp_g(n_miss,a[unique(b_id)],b,ap1,ap2)

    end

    # Return energy function properly scaled.
    return exp(-Ux(x)) * sf / Z

end
"""
    `create_Llkhd([N1,...,NK],XOBS)`

    Create function to compute the minus log-likelihood function for a region with N CpG sites given
    the M partial observations XOBS.

    # Examples
    ```julia-repl
    julia> using Random; Random.seed!(1234); n=[4];
    julia> xobs=[CpelTdm.gen_x_mc(n,[0.0,0.0]) for _ in 1:20];
    julia> LogLike=CpelTdm.create_Llkhd(n,xobs);
    julia> LogLike([0.0,0.0,0.0])
    55.45177444479562
    ```
"""
function create_Llkhd(n::Vector{Int64},xobs::Array{Vector{Int64},1})

    # Define minus log-likelihood function
    function Llkhd_fun(theta::Vector{Float64})::Float64
        # Get parameters
        aux = 0.0
        a = theta[1:(end-1)]
        b = theta[end]

        # Get energy function and partition function
        Ux = create_Ux(n,a,b)
        logZ = log(comp_Z(n,a,b))

        # Initialize variables used in for loop
        ap1 = ap2 = 0.0                 # Boundary values left and right
        b_id = n_miss = []              # Block IDs of missing CpG; number missing per ID
        zerost = zeroend = [] # Indices of zeros
        subid = vcat([i*ones(Int64,n[i]) for i in 1:length(n)]...)    # ID's of CpGs

        # Contribution of each observation
        for x in xobs

            # Find changes to/from 0 in x vector
            zerost = findall(isequal(-1),abs.(x[2:end]) - abs.(x[1:(end-1)])) .+ 1
            zeroend = findall(isequal(1),abs.(x[2:end]) - abs.(x[1:(end-1)])) .+ 1

            # Determine whether it starts/finishes with 0
            x[1]==0 && pushfirst!(zerost,1)
            x[end]==0 && push!(zeroend,length(x)+1)

            # Get scaling factors due to all missing data
            @inbounds for i in 1:length(zerost)
                # Find boundary conditions (other X's or boundary of X)
                ap1 = zerost[i]==1 ? a[1] : 2.0*x[zerost[i]-1]*b+a[subid[zerost[i]]]
                ap2 = zeroend[i]==sum(n)+1 ? a[end] : 2.0*x[zeroend[i]]*b+a[subid[zeroend[i]-1]]

                # Figure out b (block IDs of missing) and r (α indices)
                b_id = subid[zerost[i]:(zeroend[i]-1)]
                n_miss = [count(x->x==n,b_id) for n in unique(b_id)]

                # Call scaling factor function
                aux += log(comp_g(n_miss,a[unique(b_id)],b,ap1,ap2))
            end

            # Add log of it to Llkhd
            aux += -Ux(x)
        end

        # Return MINUS log-likelihood. Add as many logZ as samples we have
        -aux + length(xobs) * logZ
    end

    # Return function
    return Llkhd_fun

end
"""
    `est_alpha(XOBS)`

    Estimate parameter α in N=1 case (first entry in returned vector).

    # Examples
    ```julia-repl
    julia> using Random; Random.seed!(1234); n=[1];
    julia> xobs=[CpelTdm.gen_x_mc(n,[0.0,0.0]) for _ in 1:20];
    julia> CpelTdm.est_alpha(xobs)
    2-element Array{Float64,1}:
     0.10033534773107566
     0.0
    ```
"""
function est_alpha(xobs::Array{Vector{Int64},1})::Vector{Float64}
    # Derivation of estimate:
        # log(p)+log(2) = α-log(cosh(α)); log(1-p)+log(2) = -α-log(cosh(α))
        # α = (log(p)-log(1-p))/2

    # Proportion of X=1
    phat = length(findall(x->x==[1],xobs)) / length(xobs)
    a = 0.5*(log(phat)-log(1.0-phat))

    # Return estimate
    return [min(max(-ETA_MAX_ABS,a),ETA_MAX_ABS),0.0]

end
"""
    `est_theta_sa([N1,...,NK],XOBS)`

    Estimate parameter vector [α1,...,αK,β] using simulated annealing.

    # Examples
    ```julia-repl
    julia> using Random; Random.seed!(1234); n=[4];
    julia> xobs=[CpelTdm.gen_x_mc(n,[0.0,0.0]) for _ in 1:100];
    julia> CpelTdm.est_theta_sa(n,xobs)
    2-element Array{Float64,1}:
     -0.08515675262874925
     -0.040190205815463766
    ```
"""
function est_theta_sa(n::Vector{Int64},xobs::Array{Vector{Int64},1})::Vector{Float64}

    # If N=1, then estimate α
    sum(n)==1 && return est_alpha(xobs)

    # Define boundaries and initialization
    L = create_Llkhd(n,xobs)
    init = zeros(Float64,length(n)+1)
    lower = -ETA_MAX_ABS * ones(Float64,length(n)+1)
    upper = ETA_MAX_ABS * ones(Float64,length(n)+1)
    opts = Optim.Options(iterations=10^4,show_trace=false,store_trace=false)

    # Boxed Simulated Annealing (SAMI)
    optim = Optim.optimize(L,lower,upper,init,SAMIN(rt=1e-4;f_tol=1e-3,verbosity=0),opts)

    # Return estimate
    return optim.minimizer

end
###################################################################################################
# OUTPUT QUANTITIES
###################################################################################################
"""
    `comp_mml([N1,...,NK],∇logZ)`

    Function that computes mean methylation level (MML) over the entire region of interest.

    # Examples
    ```julia-repl
    julia> n=[1,1,1]; θ=[1.0,-1.0,1.0,0.0];
    julia> ∇logZ = CpelTdm.get_∇logZ(n,θ);
    julia> CpelTdm.comp_mml(n,∇logZ)
    0.62693236
    ```
"""
function comp_mml(n::Vector{Int64},∇logZ::Vector{Float64})::Float64

    # Return
    return abs(round(0.5*(1.0+1.0/sum(n)*sum(∇logZ[1:length(n)]));digits=8))

end
"""
    `comp_nme([N1,...,NK],θ,∇logZ)`

    Function that computes normalized methylation entropy (NME) using the ∇logZ information.

    # Examples
    ```julia-repl
    julia> n=[5]; α=[0.0]; β=0.0; θ=vcat([α,β]...);
    julia> ∇logZ=CpelTdm.get_∇logZ(n,θ);
    julia> CpelTdm.comp_nme(n,θ,∇logZ)
    1.0
    ```
"""
function comp_nme(n::Vector{Int64},θ::Vector{Float64},∇logZ::Vector{Float64})::Float64

    # Return
    return abs(round(1.0/(sum(n)*LOG2)*(log(comp_Z(n,θ[1:(end-1)],θ[end]))-θ'*∇logZ);digits=8))

end
"""
    `comp_cmd([N1,...,NK],θ1,θ2,∇logZ1s,∇logZ2s)`

    Function that computes the normalized geometric Jensen-Shannon divergence between two Ising models
    parametrized by θ1 and θ2, respectively.

    # Examples
    ```julia-repl
    julia> n=[20]; θ1=θ2=[1.0,1.0];
    julia> ∇logZ1s = CpelTdm.get_∇logZ(n,θ1);
    julia> ∇logZ2s = CpelTdm.get_∇logZ(n,θ2);
    julia> CpelTdm.comp_cmd(n,θ1,θ2,∇logZ1s,∇logZ2s)
    0.0
    julia> n=[20]; θ1=[-5.0,5.0]; θ2=[5.0,5.0];
    julia> ∇logZ1s = CpelTdm.get_∇logZ(n,θ1);
    julia> ∇logZ2s = CpelTdm.get_∇logZ(n,θ2);
    julia> CpelTdm.comp_cmd(n,θ1,θ2,∇logZ1s,∇logZ2s)
    1.0
    ```
"""
function comp_cmd(n::Vector{Int64},θ1::Vector{Float64},θ2::Vector{Float64},∇logZ1::Vector{Float64},∇logZ2::Vector{Float64})::Float64

    # Parameter geometric mixture of Ising
    θγ = 0.5.*(θ1+θ2)

    # Get partition functions
    logZ1 = log(comp_Z(n,θ1[1:(end-1)],θ1[end]))
    logZ2 = log(comp_Z(n,θ2[1:(end-1)],θ2[end]))
    logZγ = log(comp_Z(n,θγ[1:(end-1)],θγ[end]))

    # Compute numerator of 1-CMD
    cmd = logZ1 + logZ2 - (θ1'*∇logZ1 + θ2'*∇logZ2)

    # Normalize 1-CMD
    cmd /= 2.0*logZγ - θγ'* (∇logZ1 + ∇logZ2)

    # Return GJSD
    return 1.0 - cmd

end
"""
    `comp_cmd_unnorm([N1,...,NK],θ1,θ2,∇logZ1s,∇logZ2s)`

    Function that computes the unnormalized geometric Jensen-Shannon divergence between two Ising models
    parametrized by θ1 and θ2, respectively.

    # Examples
    ```julia-repl
    julia> n=[20]; θ1=θ2=[1.0,1.0];
    julia> ∇logZ1s = CpelTdm.get_∇logZ(n,θ1);
    julia> ∇logZ2s = CpelTdm.get_∇logZ(n,θ2);
    julia> CpelTdm.comp_cmd_unnorm(n,θ1,θ2,∇logZ1s,∇logZ2s)
    0.0
    julia> n=[20]; θ1=[-5.0,5.0]; θ2=[5.0,5.0];
    julia> ∇logZ1s = CpelTdm.get_∇logZ(n,θ1);
    julia> ∇logZ2s = CpelTdm.get_∇logZ(n,θ2);
    julia> CpelTdm.comp_cmd_unnorm(n,θ1,θ2,∇logZ1s,∇logZ2s)
    1.3880197058418844
    ```
"""
function comp_cmd_unnorm(n::Vector{Int64},θ1::Vector{Float64},θ2::Vector{Float64},∇logZ1::Vector{Float64},∇logZ2::Vector{Float64})::Float64

    # Parameter geometric mixture of Ising
    θγ = 0.5.*(θ1+θ2)

    # Get partition functions
    logZ1 = log(comp_Z(n,θ1[1:(end-1)],θ1[end]))
    logZ2 = log(comp_Z(n,θ2[1:(end-1)],θ2[end]))
    logZγ = log(comp_Z(n,θγ[1:(end-1)],θγ[end]))

    # Return unnormalized GJSD
    return 2.0*logZγ - θγ'* (∇logZ1 + ∇logZ2) - (logZ1 + logZ2 - (θ1'*∇logZ1 + θ2'*∇logZ2))

end

ext: jl


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(3*(5+-4*cos(x))^-1, x) == :(x+2*arctan((2+-1*cos(x))^-1*sin(x)))
@test integrate((3+cos(x)^2+2*sin(x)+-2*cos(x)*sin(x))^-1*(1+2*sin(x)+cos(x)), x) == :(-1*arctan((2+sin(x))^-1*(-1*sin(x)+2*cos(x))))
@test integrate((-2*sin(x)^2+-2*sin(x)+4*cos(x)+cos(x)*sin(x))^-1*(2+5*sin(x)+cos(x)), x) == :(-1*log(1+-3*cos(x)+sin(x))+log(3+cos(x)+sin(x)))
@test integrate((1+-5*sin(x)+3*cos(x)^2+4*cos(x)+-1*cos(x)*sin(x))^-1*(3+2*sin(x)+7*cos(x)), x) == :(-1*log(1+-2*sin(x)+cos(x))+log(3+cos(x)+sin(x)))
@test integrate((-1+-4*cos(x)+-3*cos(x)^2+4*cos(x)^3)^-1*(-1+4*cos(x)+5*cos(x)^2), x) == :(x+-2*arctan((3+cos(x))^-1*sin(x))+-2*arctan((1+2*cos(x)+5*cos(x)^2)^-1*(3*sin(x)+7*cos(x)*sin(x))))
@test integrate((-1+-9*cos(x)^2+2*cos(x)+4*cos(x)^3)^-1*(-5+2*cos(x)+7*cos(x)^2), x) == :(x+-2*arctan(2*(1+-1*cos(x)+2*cos(x)^2)^-1*cos(x)*sin(x)))
@test integrate(3*(5+4*sin(x))^-1, x) == :(x+2*arctan((2+sin(x))^-1*cos(x)))
@test integrate(2*(1+cos(x)^2)^-1, x) == :(x*2^(1/2)+-1*2^(1/2)*arctan((1+2^(1/2)+cos(x)^2)^-1*cos(x)*sin(x)))
@test integrate((p+q*cos(x)+r*sin(x))^-1, x) == :(2*(p^2+-1*q^2+-1*r^2)^-1//2*arctan((p^2+-1*q^2+-1*r^2)^-1//2*(r+(p+-1q)*tan((1/2)*x))))

max_stars_repo_path: test/0 Independent test suites/Jeffrey Problems.jl


text:
#AudioStick.jl: Copy/synchronize playlist media files to a given drive/folder.
module AudioStick

import Printf: @sprintf
using SHA

include("files.jl")
include("base.jl")
include("console.jl")
include("console_app.jl")
include("winreg.jl")

function run_app()
	try
		run_console()
	catch e
		@error(e)
#		rethrow(e) #debug
	end

	pause() #Allow user to see results before closing terminal.
end

end
#Last line

max_stars_repo_name: ma-laforge/AudioStick

max_issues_repo_name: ma-laforge/AudioStick

max_forks_repo_name: ma-laforge/AudioStick


text:
struct Simulation{DT <: Number, IT <: Integrator{DT}, FPRE <: Union{Function,Nothing}, FPOST <: Union{Function,Nothing}}
    int::IT
    nt::Int

    f_pre::FPRE
    f_post::FPOST

    x::Vector{Vector{DT}}

    function Simulation(x₀::Vector{Vector{DT}}, int::IT, nt::Int, f_preproc::FPRE, f_postproc::FPOST) where {DT, IT <: Integrator{DT}, FPRE, FPOST}
        new{DT, IT, FPRE, FPOST}(int, nt, f_preproc, f_postproc, copy(x₀))
    end
end

function Simulation(equ::Equation, integrator::Type{<:Integrator}, Δt::Real, nt::Int; f_preproc=nothing, f_postproc=nothing, kwargs...)
    int = integrator(equ, Δt; kwargs...)
    Simulation(equ.x₀, int, nt, f_preproc, f_postproc)
end

ndims(sim::Simulation) = ndims(sim.int.equ)
nsamples(sim::Simulation) = nsamples(sim.int.equ)
ntimesteps(sim::Simulation) = sim.nt

eachsample(sim::Simulation) = eachindex(sim.x)
eachtimestep(sim::Simulation) = 1:ntimesteps(sim)

haspreprocessing(::Simulation{<:Any, <:Any, Nothing}) = false
haspreprocessing(::Simulation{<:Any, <:Any, <:Function}) = true

haspostprocessing(::Simulation{<:Any, <:Any, <:Any, Nothing}) = false
haspostprocessing(::Simulation{<:Any, <:Any, <:Any, <:Function}) = true

preprocessing(sim::Simulation, x, n) = haspreprocessing(sim) && sim.f_pre(x, n)
postprocessing(sim::Simulation, x, n) = haspostprocessing(sim) && sim.f_post(x, n)

function (sim::Simulation{DT})() where {DT}
    for n in eachtimestep(sim)
        preprocessing(sim, sim.x, n)
        sim.int(sim.x)
        postprocessing(sim, sim.x, n)
    end
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
ByteLevelBPE(bfile::AbstractString; kws...) = open(io->ByteLevelBPE(io; kws...), bfile)

function ByteLevelBPE(io::IO;
                      glossaries = nothing,
                      merge::Int = -1, sepsym = nothing, endsym = nothing,
                      has_header::Bool = true,
                      codemap = default_codemap(),
                      input_transform = gpt2_tokenizer,
                      output_transform = nothing,
                      normalizer=nothing)

  cache = Dict{String, String}()
  oldendsym, merging_rank = read_bpefile(io, merge=merge, has_header=has_header)
  glossaries !== nothing && (glossaries = Glossary(glossaries))
  return GenericBPE{String}(; sepsym, oldendsym, endsym, input_transform, output_transform, merging_rank, cache, glossaries, codemap, normalizer)
end


text:
using LinearAlgebra, Test
using StaticArrays
using Polyhedra

function test_cartesian_product()
    p1 = polyhedron(hrep([HalfSpace([1.0], 1.0), HalfSpace([-1.0], 0.0)]))
    p2 = polyhedron(hrep([HalfSpace([1.0], 1.0), HalfSpace([-1.0], 0.0)]))
    p = Polyhedra.hcartesianproduct(p1, p2)
    @test p isa DefaultPolyhedron{Float64,
                                  Polyhedra.Intersection{Float64,StaticArrays.SArray{Tuple{2},Float64,1,2},StaticArrays.Size{(2,)}},
                                  Polyhedra.Hull{Float64,StaticArrays.SArray{Tuple{2},Float64,1,2},StaticArrays.Size{(2,)}}}
    expected = HalfSpace((@SVector [1.0, 0.0]), 1.0) ∩ HalfSpace((@SVector [-1.0, 0.0]), 0.0) ∩ HalfSpace((@SVector [0.0, 1.0]), 1.0) ∩ HalfSpace((@SVector [0.0, -1.0]), 0.0)
    inequality_fulltest(p, expected)
end

function _test_sethrep(fh, fv, args...)
    #             x_1 ≤ 2   x_1 ≥ -1 <=> -x_1 ≤ 1
    h1 = HalfSpace([1], 2) ∩ HalfSpace([-1], 1)
    v1 = convexhull([-1], [2])

    ph = polyhedron(h1)
    @test volume(ph) == 3
    inequality_fulltest(ph, h1)
    generator_fulltest(ph, v1)

    pv = polyhedron(v1)
    @test volume(pv) == 3
    inequality_fulltest(pv, h1)
    generator_fulltest(pv, v1)

    #             x_1 ≤ 1   x_1 ≥ -1 <=> -x_1 ≤ 1
    h2 = HalfSpace([1], 1) ∩ HalfSpace([-1], 1)
    v2 = convexhull([-1], [1])

    # Test that the three fields have been modified
    fh(ph, h2, args...)
    @test volume(ph) == 2
    inequality_fulltest(ph, h2)
    generator_fulltest(ph, v2)
    fv(pv, v2, args...)
    @test volume(pv) == 2
    inequality_fulltest(pv, h2)
    generator_fulltest(pv, v2)
end
function test_sethrep()
    _test_sethrep(Polyhedra.sethrep!, Polyhedra.setvrep!)
    _test_sethrep(Polyhedra.sethrep!, Polyhedra.setvrep!, Polyhedra.UNKNOWN_REDUNDANCY)
    _test_sethrep(Polyhedra.resethrep!, Polyhedra.resetvrep!)
    _test_sethrep(Polyhedra.resethrep!, Polyhedra.resetvrep!, Polyhedra.UNKNOWN_REDUNDANCY)
end

function test_redundancy()
    h = HalfSpace([1], 2) ∩ HalfSpace([-1], 1)
    p = polyhedron(h)
    function f()
        removevredundancy!(p)
        removevredundancy!(p, nothing; strongly = false, planar = false)
        detectvlinearity!(p)
        detectvlinearity!(p, nothing; verbose=1)
        removehredundancy!(p)
        removehredundancy!(p, nothing; strongly = false, planar = false)
        detecthlinearity!(p)
        detecthlinearity!(p, nothing; verbose=1)
    end
    alloc_test(f, 0)
end

@testset "Interval tests" begin
    # Closed interval
    h = HalfSpace([-1], 3.) ∩ HalfSpace([1.], 8)
    v = convexhull([-3.], [8.])

    v0 = [[1.], [8], [-2], [4], [-1], [-3], [5]]
    pp = polyhedron(vrep(v0), Polyhedra.DefaultLibrary{Float64}())
    @test !hrepiscomputed(pp)
    @test vrepiscomputed(pp)
    d = StaticArrays.Size{(1,)}()
    for p in (Interval{Float64, SVector{1, Float64}, typeof(d)}(pp),
              Interval{Float64, SVector{1, Float64}, typeof(d)}(d, Polyhedra.vreps(pp)...))
        @test similar_library(pp, 1) == IntervalLibrary{Float64}()
        @test similar_library(pp, StaticArrays.Size((1,))) == IntervalLibrary{Float64}()
        @test library(p) == IntervalLibrary{Float64}()
        @test similar_library(p, 2) == Polyhedra.DefaultLibrary{Float64}()
        @test similar_library(p, StaticArrays.Size((2,))) == Polyhedra.DefaultLibrary{Float64}()
        @test hrepiscomputed(p)
        @test vrepiscomputed(p)
        @test p isa Interval{Float64}
        @test !isempty(p)
        @test volume(p) == 11
        @test center_of_mass(p) ≈ [2.5]
        @test dim(p) == 1
        inequality_fulltest(p, h)
        generator_fulltest(p, v)
    end

    pp = polyhedron(h, Polyhedra.DefaultLibrary{Float64}())
    @test hrepiscomputed(pp)
    @test !vrepiscomputed(pp)
    for p in (Interval{Float64, SVector{1, Float64}, typeof(d)}(pp),
              Interval{Float64, SVector{1, Float64}, typeof(d)}(d, Polyhedra.hreps(pp)...))
        @test hrepiscomputed(p)
        @test vrepiscomputed(p)
        @test p isa Interval{Float64}
        @test !isempty(p)
        @test iszero(surface(p))
        @test volume(p) == 11
        @test center_of_mass(p) ≈ [2.5]
        @test dim(p) == 1
        generator_fulltest(p, v)
        inequality_fulltest(p, h)
    end

    # Singleton
    h = intersect(HyperPlane([1.], 2))
    v = convexhull([2.])

    v0 = reshape([2., 2., 2.], 3, 1)
    p = polyhedron(vrep(v0))
    @test p isa Interval{Float64}
    @test !isempty(p)
    @test iszero(surface(p))
    @test iszero(volume(p))
    @test dim(p) == 0
    inequality_fulltest(p, h)
    generator_fulltest(p, v)

    p = polyhedron(h)
    @test p isa Interval{Float64}
    @test !isempty(p)
    @test iszero(surface(p))
    @test iszero(volume(p))
    @test dim(p) == 0
    generator_fulltest(p, v)
    inequality_fulltest(p, h)

    # Left-open interval
    h = intersect(HalfSpace([1.], 2))
    v = v + conichull([-1.])

    r0 = reshape([-1., -1], 2, 1)
    p = polyhedron(vrep(v0, r0))
    @test p isa Interval{Float64}
    @test !isempty(p)
    @test iszero(surface(p))
    @test dim(p) == 1
    inequality_fulltest(p, h)
    generator_fulltest(p, v)

    p = polyhedron(HalfSpace([1.], 3) ∩ HalfSpace([5.], 10) ∩ HalfSpace([2.], 6))
    @test p isa Interval{Float64}
    @test !isempty(p)
    @test iszero(surface(p))
    @test dim(p) == 1
    generator_fulltest(p, v)
    inequality_fulltest(p, h)

    # Right-open interval
    h = intersect(HalfSpace([-1.], -2))
    v = convexhull([2.]) + conichull([1.])

    p = polyhedron(convexhull([2.], [3.]) + conichull([1.], [3.]))
    @test p isa Interval{Float64}
    @test !isempty(p)
    @test iszero(surface(p))
    @test dim(p) == 1
    inequality_fulltest(p, h)
    generator_fulltest(p, v)

    p = polyhedron(HalfSpace([-3.], -6) ∩ HalfSpace([-2.], -2))
    # TODO test value of p
    p = polyhedron(h)
    @test p isa Interval{Float64}
    @test !isempty(p)
    @test iszero(surface(p))
    @test dim(p) == 1
    generator_fulltest(p, v)
    inequality_fulltest(p, h)

    # Empty
    h = HyperPlane([0], 1) ∩ HyperPlane([1], 0)
    v = vrep(Line{Int, StaticArrays.SVector{1, Int}}[])

    p = polyhedron(v)
    @test p isa Interval{Int}
    @test isempty(p)
    @test iszero(surface(p))
    @test iszero(volume(p))
    @test dim(p) == -1
    inequality_fulltest(p, h)
    generator_fulltest(p, v)

    p = polyhedron(HalfSpace([1], 1) ∩ HalfSpace([-1], -2))
    @test p isa Interval{Int}
    @test isempty(p)
    @test iszero(surface(p))
    @test iszero(volume(p))
    @test dim(p) == -1
    generator_fulltest(p, v)
    inequality_fulltest(p, h)

    p = polyhedron(intersect(HyperPlane([0], 1)))
    @test p isa Interval{Int}
    @test isempty(p)
    @test iszero(surface(p))
    @test iszero(volume(p))
    @test dim(p) == -1
    generator_fulltest(p, v)
    inequality_fulltest(p, h)

    p = polyhedron(intersect(HalfSpace([0], -1)))
    @test p isa Interval{Int}
    @test isempty(p)
    @test iszero(surface(p))
    @test iszero(volume(p))
    @test dim(p) == -1
    generator_fulltest(p, v)
    inequality_fulltest(p, h)

    # Line
    h = hrep(HyperPlane{Int, SVector{1, Int}}[])
    v = conichull(Line([1]))

    p = polyhedron(v)
    @test p isa Interval{Int}
    @test !isempty(p)
    @test iszero(surface(p))
    @test dim(p) == 1
    inequality_fulltest(p, h)
    generator_fulltest(p, v)

    p = polyhedron(h)
    @test p isa Interval{Int}
    @test !isempty(p)
    @test iszero(surface(p))
    @test dim(p) == 1
    generator_fulltest(p, v)
    inequality_fulltest(p, h)

    # Symmetric interval
    h = HalfSpace([1], 1) ∩ HalfSpace([-1], 1)
    v = convexhull([1], [-1])

    p = polyhedron(convexhull([-1], [1], [0]))
    @test p isa Interval{Int}
    @test !isempty(p)
    @test iszero(surface(p))
    @test volume(p) == 2
    @test dim(p) == 1
    inequality_fulltest(p, h)
    generator_fulltest(p, v)

    p = polyhedron(h)
    @test p isa Interval{Int}
    @test !isempty(p)
    @test iszero(surface(p))
    @test volume(p) == 2
    @test dim(p) == 1
    generator_fulltest(p, v)
    inequality_fulltest(p, h)

    @testset "Cartesian product (#132)" begin
        test_cartesian_product()
    end

    @testset "sethrep! (#267)" begin
        test_sethrep()
    end

    @testset "test_redundancy (#298)" begin
        test_redundancy()
    end
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
function get_pv_coordinates(obj::ExtendedPVCoordinatesProvider, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end

function get_pv_coordinates(obj::PVCoordinatesProvider, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function to_field_pv_coordinates_provider(obj::ExtendedPVCoordinatesProvider, arg0::Field)
    return jcall(obj, "toFieldPVCoordinatesProvider", FieldPVCoordinatesProvider, (Field,), arg0)
end


text:
@cache mutable struct GenericImplicitEulerCache{uType,DiffCacheType,uNoUnitsType,rateType,rhsType,nl_rhsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  dual_cache::DiffCacheType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  fsalfirst::rateType
  rhs::rhsType
  nl_rhs::nl_rhsType
end

function alg_cache(alg::GenericImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits)
  k = zero(rate_prototype)
  dual_cache = DiffCache(u,Val{determine_chunksize(u,get_chunksize(alg.nlsolve))})
  rhs = ImplicitRHS(f,tmp,t,t,t,dual_cache,p)
  fsalfirst = zero(rate_prototype)
  nl_rhs = alg.nlsolve(Val{:init},rhs,u)

  GenericImplicitEulerCache{typeof(u),typeof(dual_cache),
                     typeof(atmp),typeof(k),typeof(rhs),typeof(nl_rhs)}(
                     u,uprev,uprev2,dual_cache,tmp,atmp,k,fsalfirst,rhs,nl_rhs)
end

struct GenericImplicitEulerConstantCache{vecuType,rhsType,nl_rhsType} <: OrdinaryDiffEqConstantCache
  uhold::vecuType
  rhs::rhsType
  nl_rhs::nl_rhsType
end

function alg_cache(alg::GenericImplicitEuler,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uBottomEltypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  uhold = Vector{typeof(u)}(undef, 1)
  rhs = ImplicitRHS_Scalar(f,zero(u),t,t,t,p)
  nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
  GenericImplicitEulerConstantCache{typeof(uhold),typeof(rhs),typeof(nl_rhs)}(uhold,rhs,nl_rhs)
end

@cache mutable struct GenericTrapezoidCache{uType,DiffCacheType,uNoUnitsType,rateType,rhsType,nl_rhsType,tType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  fsalfirst::rateType
  dual_cache::DiffCacheType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  rhs::rhsType
  nl_rhs::nl_rhsType
  uprev3::uType
  tprev2::tType
end

function alg_cache(alg::GenericTrapezoid,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  dual_cache = DiffCache(u,Val{determine_chunksize(u,get_chunksize(alg.nlsolve))})
  rhs = ImplicitRHS(f,tmp,t,t,t,dual_cache,p)
  nl_rhs = alg.nlsolve(Val{:init},rhs,u)
  uprev3 = zero(u)
  tprev2 = t
  GenericTrapezoidCache{typeof(u),typeof(dual_cache),typeof(atmp),typeof(k),
                        typeof(rhs),typeof(nl_rhs),typeof(t)}(
                        u,uprev,uprev2,fsalfirst,
                        dual_cache,tmp,atmp,k,rhs,nl_rhs,uprev3,tprev2)
end

@cache mutable struct GenericTrapezoidConstantCache{vecuType,rhsType,nl_rhsType,uType,tType} <: OrdinaryDiffEqConstantCache
  uhold::vecuType
  rhs::rhsType
  nl_rhs::nl_rhsType
  uprev3::uType
  tprev2::tType
end

function alg_cache(alg::GenericTrapezoid,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  uhold = Vector{typeof(u)}(undef, 1)
  rhs = ImplicitRHS_Scalar(f,zero(u),t,t,t,p)
  nl_rhs = alg.nlsolve(Val{:init},rhs,uhold)
  uprev3 = u
  tprev2 = t
  GenericTrapezoidConstantCache(uhold,rhs,nl_rhs,uprev3,tprev2)
end

get_chunksize(cache::GenericImplicitEulerCache{uType,DiffCacheType,rateType,CS}) where {uType,DiffCacheType,rateType,CS} = CS
get_chunksize(cache::GenericTrapezoidCache{uType,DiffCacheType,rateType,CS}) where {uType,DiffCacheType,rateType,CS} = CS


text:
## Tests for correct sketching
using LowRankSketch
using FactCheck

TOL = 1e-1

# we'll sketch rank-1 matrices with a two dimensional sketch
m,n,k,s = 100,200,3,5
# m,n,k,s = 10,20,3,3
srand(1)

facts("Sketch") do

  context("sketch real matrix with GaussDimRedux") do
    sketch = Sketch(m,n,k,s,GaussDimRedux,Float64)
    M = rand(m,1)*rand(1,n)
    linear_update(sketch, M)
    Q,W,P = low_rank_approx(sketch)
    @fact vecnorm(M - Q*W*P')/vecnorm(M) --> roughly(0, TOL)
    U,Σ,V = fixed_rank_approx(sketch,1)
    @fact vecnorm(M - U*Σ*V')/vecnorm(M) --> roughly(0, TOL)
  end

  context("sketch real matrix with SSRFTDimRedux") do
    sketch = Sketch(m,n,k,s,SSRFTDimRedux,Float64)
    M = rand(m,1)*rand(1,n)
    linear_update(sketch, M)
    Q,W,P = low_rank_approx(sketch)
    @fact vecnorm(M - Q*W*P')/vecnorm(M) --> roughly(0, TOL)
    U,Σ,V = fixed_rank_approx(sketch,1)
    @fact vecnorm(M - U*Σ*V')/vecnorm(M) --> roughly(0, TOL)
  end

  context("sketch complex matrix with GaussDimRedux") do
    sketch = Sketch(m,n,k,s,GaussDimRedux,Complex{Float64})
    M = (rand(m,1) + im*rand(m,1))*(rand(1,n) + im*rand(1,n))
    linear_update(sketch, M)
    Q,W,P = low_rank_approx(sketch)
    @fact vecnorm(M - Q*W*P')/vecnorm(M) --> roughly(0, TOL)
    # recovery of rank 1 complex matrix fails with rank 1 fixed rank approx, works with rank 2
    U,Σ,V = fixed_rank_approx(sketch,2)
    @fact vecnorm(M - U*Σ*V')/vecnorm(M) --> roughly(0, TOL)
  end

  context("sketch complex matrix with SSRFTDimRedux") do
    sketch = Sketch(m,n,k,s,SSRFTDimRedux,Complex{Float64})
    M = (rand(m,1) + im*rand(m,1))*(rand(1,n) + im*rand(1,n))
    linear_update(sketch, M)
    Q,W,P = low_rank_approx(sketch)
    @fact vecnorm(M - Q*W*P')/vecnorm(M) --> roughly(0, TOL)
    U,Σ,V = fixed_rank_approx(sketch,1)
    @fact vecnorm(M - U*Σ*V')/vecnorm(M) --> roughly(0, TOL)
  end

  context("sketch integer matrix with SSRFTDimRedux") do
    # (make sure to use sketch compatible with type of matrix to be sketched, eg integer matrix -> use SSRFT)
    # nb linear update doesn't work yet for this (InexactError) b/c output of SSRFT is not integral
    sketch = Sketch(m,n,k,s,SSRFTDimRedux,Complex{Int})
    sketch = Sketch(m,n,k,s,SSRFTDimRedux,Int)
  end

end


text:
##########################################################################
#
# Benchmark optimization algorithms by tracking:
#
# * Run time over 100 runs -- after 1 initial run that forces JIT.
# * Number of iterations
# * Euclidean error of solution
# * Memory requirements (TODO)
#
##########################################################################

using Optim

algorithms = [:gradient_descent,
              :momentum_gradient_descent,
              :bfgs,
              :l_bfgs,
              :cg]

# Print out a header line for the TSV-formatted report
println(join({"Problem",
              "Algorithm",
              "AverageRunTimeInMilliseconds",
              "Iterations",
              "Error",
              "Autodiff"},
              "\t"))

for (name, problem) in Optim.UnconstrainedProblems.examples
    for algorithm in algorithms
        for ad in [false, true]
        if !problem.isdifferentiable &&
              (algorithm == :gradient_descent ||
               algorithm == :momentum_gradient_descent ||
               algorithm == :bfgs ||
               algorithm == :l_bfgs ||
               algorithm == :cg)
            continue
        end

        # Force compilation
        results = optimize(problem.f,
                           problem.initial_x,
                           method = algorithm,
                           g_tol = 1e-8,
                           autodiff = ad)

        # Run each algorithm 1,000 times
        n = 1_000

        # Estimate run time in seconds
        run_time = @elapsed for i = 1:n
            results = optimize(problem.f,
                               problem.initial_x,
                               method = algorithm,
                               g_tol = 1e-8,
                               autodiff = ad)
        end

        # Estimate error in discovered solution
        results = optimize(problem.f,
                           problem.initial_x,
                           method = algorithm,
                           g_tol = 1e-8,
                           autodiff = ad)
        errors = minimum(map(sol -> norm(results.minimum - sol), problem.solutions))

        # Count iterations
        iterations = results.iterations

        # Print out results.
        println(join({problem.name,
                      results.method,
                      run_time,
                      iterations,
                      errors,
                      ad},
                      "\t"))
        end
    end
end


text:
using Trying
using Test

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


text:
using CSV, DataFrames, Statistics
using FortranFiles, MeshArrays, MITgcmTools

##

function get_from_dataverse(nam::String,pth::String)
    tmp = CSV.File("nctiles_climatology.csv") |> DataFrame!
    ii = findall([occursin("$nam", tmp[i,:name]) for i=1:size(tmp,1)])
    !isdir("$pth"*"$nam") ? mkdir("$pth"*"$nam") : nothing
    for i in ii
        ID=tmp[i,:ID]
        nam1=tmp[i,:name]
        nam2=joinpath("$pth"*"$nam/",nam1)
        run(`wget --content-disposition https://dataverse.harvard.edu/api/access/datafile/$ID`);
        run(`mv $nam1 $nam2`);
    end
end

function get_grid_if_needed()
    if !isdir("../inputs/GRID_LLC90")
        run(`git clone https://github.com/gaelforget/GRID_LLC90 ../inputs/GRID_LLC90`)
    end
end

function get_velocity_if_needed()
    pth="../inputs/nctiles_climatology/"
    !isdir("$pth") ? mkdir("$pth") : nothing
    !isdir("$pth"*"UVELMASS") ? get_from_dataverse("UVELMASS",pth) : nothing
    !isdir("$pth"*"VVELMASS") ? get_from_dataverse("VVELMASS",pth) : nothing
end

function read_velocities(mygrid::gcmgrid,t::Int)
    pth="../inputs/nctiles_climatology/"
    u=Main.read_nctiles("$pth"*"UVELMASS/UVELMASS","UVELMASS",mygrid,I=(:,:,:,t))
    v=Main.read_nctiles("$pth"*"VVELMASS/VVELMASS","VVELMASS",mygrid,I=(:,:,:,t))
    return u,v
end

#Convert Velocity (m/s) to transport (m^3/s)
function convert_velocities(U,V,γ)
    for i in eachindex(U)
        tmp1=U[i]; tmp1[(!isfinite).(tmp1)] .= 0.0
        tmp1=V[i]; tmp1[(!isfinite).(tmp1)] .= 0.0
        U[i]=γ["DRF"][i[2]]*U[i].*γ["DYG"][i[1]]
        V[i]=γ["DRF"][i[2]]*V[i].*γ["DXG"][i[1]]
    end
    return U,V
end

##

"""
    trsp_read(myspec::String,mypath::String)

Function that reads files that were generated by `trsp_prep`
"""
function trsp_read(myspec::String,mypath::String)
    mygrid=GridSpec(myspec,mypath)
    TrspX=mygrid.read(mypath*"TrspX.bin",MeshArray(mygrid,Float32))
    TrspY=mygrid.read(mypath*"TrspY.bin",MeshArray(mygrid,Float32))
    TauX=mygrid.read(mypath*"TauX.bin",MeshArray(mygrid,Float32))
    TauY=mygrid.read(mypath*"TauY.bin",MeshArray(mygrid,Float32))
    SSH=mygrid.read(mypath*"SSH.bin",MeshArray(mygrid,Float32))
    return TrspX, TrspY, TauX, TauY, SSH
end

"""
    trsp_prep(mygrid,GridVariables,dirOut)

Function that generates small binary files (2D) from large netcdf ones (4D).

```
using FortranFiles, MeshArrays
!isdir("nctiles_climatology") ? error("missing files") : nothing
include(joinpath(dirname(pathof(MeshArrays)),"gcmfaces_nctiles.jl"))
(TrspX, TrspY, TauX, TauY, SSH)=trsp_prep(mygrid,GridVariables,"GRID_LLC90/");
```
"""
function trsp_prep(mygrid::gcmgrid,GridVariables::Dict,dirOut::String="")

    #wind stress
    fileName="nctiles_climatology/oceTAUX/oceTAUX"
    oceTAUX=read_nctiles(fileName,"oceTAUX",mygrid)
    fileName="nctiles_climatology/oceTAUY/oceTAUY"
    oceTAUY=read_nctiles(fileName,"oceTAUY",mygrid)
    oceTAUX=mask(oceTAUX,0.0)
    oceTAUY=mask(oceTAUY,0.0)

    #sea surface height anomaly
    fileName="nctiles_climatology/ETAN/ETAN"
    ETAN=read_nctiles(fileName,"ETAN",mygrid)
    fileName="nctiles_climatology/sIceLoad/sIceLoad"
    sIceLoad=read_nctiles(fileName,"sIceLoad",mygrid)
    rhoconst=1029.0
    myssh=(ETAN+sIceLoad./rhoconst)
    myssh=mask(myssh,0.0)

    #seawater transports
    fileName="nctiles_climatology/UVELMASS/UVELMASS"
    U=read_nctiles(fileName,"UVELMASS",mygrid)
    fileName="nctiles_climatology/VVELMASS/VVELMASS"
    V=read_nctiles(fileName,"VVELMASS",mygrid)
    U=mask(U,0.0)
    V=mask(V,0.0)

    #time averaging and vertical integration
    TrspX=similar(GridVariables["DXC"])
    TrspY=similar(GridVariables["DYC"])
    TauX=similar(GridVariables["DXC"])
    TauY=similar(GridVariables["DYC"])
    SSH=similar(GridVariables["XC"])

    for i=1:mygrid.nFaces
        tmpX=mean(U.f[i],dims=4)
        tmpY=mean(V.f[i],dims=4)
        for k=1:length(GridVariables["RC"])
            tmpX[:,:,k]=tmpX[:,:,k].*GridVariables["DYG"].f[i]
            tmpX[:,:,k]=tmpX[:,:,k].*GridVariables["DRF"][k]
            tmpY[:,:,k]=tmpY[:,:,k].*GridVariables["DXG"].f[i]
            tmpY[:,:,k]=tmpY[:,:,k].*GridVariables["DRF"][k]
        end
        TrspX.f[i]=dropdims(sum(tmpX,dims=3),dims=(3,4))
        TrspY.f[i]=dropdims(sum(tmpY,dims=3),dims=(3,4))
        TauX.f[i]=dropdims(mean(oceTAUX.f[i],dims=3),dims=3)
        TauY.f[i]=dropdims(mean(oceTAUY.f[i],dims=3),dims=3)
        SSH.f[i]=dropdims(mean(myssh.f[i],dims=3),dims=3)
    end

    if !isempty(dirOut)
        write_bin(TrspX,dirOut*"TrspX.bin")
        write_bin(TrspY,dirOut*"TrspY.bin")
        write_bin(TauX,dirOut*"TauX.bin")
        write_bin(TauY,dirOut*"TauY.bin")
        write_bin(SSH,dirOut*"SSH.bin")
    end

    return TrspX, TrspY, TauX, TauY, SSH
end

"""
    trsp_prep(mygrid,GridVariables,dirOut)

Function that writes a `MeshArray` to a binary file using `FortranFiles`.
"""
function write_bin(inFLD::MeshArray,filOut::String)
    recl=prod(inFLD.grid.ioSize)*4
    tmp=Float32.(convert2gcmfaces(inFLD))
    println("saving to file: "*filOut)
    f =  FortranFile(filOut,"w",access="direct",recl=recl,convert="big-endian")
    write(f,rec=1,tmp)
    close(f)
end

##

"""
    rotate_uv(uv,γ)

    1. Convert to `Sv` units and mask out land
    2. Interpolate `x/y` transport to grid cell center
    3. Convert to `Eastward/Northward` transport
    4. Display Subdomain Arrays (optional)
"""
function rotate_uv(uv,γ)
    u=1e-6 .*uv["U"]; v=1e-6 .*uv["V"];
    u[findall(γ["hFacW"][:,1].==0)].=NaN
    v[findall(γ["hFacS"][:,1].==0)].=NaN;

    nanmean(x) = mean(filter(!isnan,x))
    nanmean(x,y) = mapslices(nanmean,x,dims=y)
    (u,v)=exch_UV(u,v); uC=similar(u); vC=similar(v)
    for iF=1:u.grid.nFaces
        tmp1=u[iF][1:end-1,:]; tmp2=u[iF][2:end,:]
        uC[iF]=reshape(nanmean([tmp1[:] tmp2[:]],2),size(tmp1))
        tmp1=v[iF][:,1:end-1]; tmp2=v[iF][:,2:end]
        vC[iF]=reshape(nanmean([tmp1[:] tmp2[:]],2),size(tmp1))
    end

    cs=γ["AngleCS"]
    sn=γ["AngleSN"]
    u=uC.*cs-vC.*sn
    v=uC.*sn+vC.*cs;

    return u,v,uC,vC
end

"""
    interp_uv(u,v)
"""
function interp_uv(u,v)
    mypath="../inputs/GRID_LLC90/"
    SPM,lon,lat=read_SPM(mypath) #interpolation matrix (sparse)
    uI=MatrixInterp(write(u),SPM,size(lon)) #interpolation itself
    vI=MatrixInterp(write(v),SPM,size(lon)); #interpolation itself
    return transpose(uI),transpose(vI),vec(lon[:,1]),vec(lat[1,:])
end

"""
    read_llc90_grid()
"""
function read_llc90_grid()
    mypath="../inputs/GRID_LLC90/"
    mygrid=GridSpec("LatLonCap",mypath)
    γ=GridLoad(mygrid)
end

"""
    initialize_locations(uv_etc::Dict)

Define `uInitS` as an array of initial conditions
"""
function initialize_locations(uv_etc::Dict,n_subset::Int=1)
    msk=uv_etc["msk"]
    uInitS = Array{Float64,2}(undef, 3, prod(msk.grid.ioSize))

    kk = 0
    for fIndex = 1:5
        nx, ny = msk.fSize[fIndex]
        ii1 = 0.5:1.0:nx
        ii2 = 0.5:1.0:ny
        n1 = length(ii1)
        n2 = length(ii2)
        for i1 in eachindex(ii1)
            for i2 in eachindex(ii2)
                if msk[fIndex][Int(round(i1+0.5)),Int(round(i2+0.5))]
                    kk += 1
                    let kk = kk
                        uInitS[1, kk] = ii1[i1]
                        uInitS[2, kk] = ii2[i2]
                        uInitS[3, kk] = fIndex
                    end
                end
            end
        end
    end

    uInitS=uInitS[:,1:kk]
    du=fill(0.0,size(uInitS));

    uInitS=uInitS[:,1:n_subset:end]
    du=du[:,1:n_subset:end]
    return uInitS,du
end

"""
    postprocess_ODESolution()

Copy `sol` to a `DataFrame` & map position to lon,lat coordinates
"""
function postprocess_ODESolution(sol,uv_etc::Dict)
    XC=uv_etc["XC"]; YC=uv_etc["YC"]
    ID=collect(1:size(sol,2))*ones(1,size(sol,3))
    x=sol[1,:,:]
    y=sol[2,:,:]
    fIndex=sol[3,:,:]
    df = DataFrame(ID=Int.(ID[:]), x=x[:], y=y[:], fIndex=fIndex[:])

    lon=Array{Float64,1}(undef,size(df,1)); lat=similar(lon)

    for ii=1:length(lon)
        #get location in grid index space
        x=df[ii,:x]; y=df[ii,:y]; fIndex=Int(df[ii,:fIndex])
        dx,dy=[x - floor(x),y - floor(y)]
        i_c,j_c = Int32.(floor.([x y])) .+ 2
        #interpolate lon and lat to position
        tmp=view(YC[fIndex],i_c:i_c+1,j_c:j_c+1)
        lat[ii]=(1.0-dx)*(1.0-dy)*tmp[1,1]+dx*(1.0-dy)*tmp[2,1]+(1.0-dx)*dy*tmp[1,2]+dx*dy*tmp[2,2]

        tmp=view(XC[fIndex],i_c:i_c+1,j_c:j_c+1)
        if (maximum(tmp)>minimum(tmp)+180)&&(lat[ii]<88)
            tmp1=deepcopy(tmp)
            tmp1[findall(tmp.<maximum(tmp)-180)] .+= 360.
            tmp=tmp1
        end
        #kk=findall(tmp.<maximum(tmp)-180); tmp[kk].=tmp[kk].+360.0
        lon[ii]=(1.0-dx)*(1.0-dy)*tmp[1,1]+dx*(1.0-dy)*tmp[2,1]+(1.0-dx)*dy*tmp[1,2]+dx*dy*tmp[2,2]
    end

    df.lon=lon; df.lat=lat; #show(df[end-3:end,:])
    return df
end

##

function read_uv_etc(k::Int,γ::Dict)
    nt=12; msk=(γ["hFacC"][:,k] .> 0.) #select depth

    u=0. *γ["XC"]; v=0. *γ["XC"];
    for t=1:nt
        (U,V)=read_velocities(γ["XC"].grid,t)
        for i=1:size(u,1)
            u[i]=u[i] + U[i,k]
            v[i]=v[i] + V[i,k]
        end
    end
    u=u ./ nt
    v=v ./ nt #time average

    u[findall(isnan.(u))]=0.0; v[findall(isnan.(v))]=0.0 #mask with 0s rather than NaNs
    u=u./γ["DXC"]; v=v./γ["DYC"]; #normalize to grid units

    (u,v)=exchange(u,v,1) #add 1 point at each edge for u and v
    XC=exchange(γ["XC"]) #add 1 lon point at each edge
    YC=exchange(γ["YC"]) #add 1 lat point at each edge

    t0=0.0; t1=86400*366*10.0; dt=10*86400.0;
    uv_etc = Dict("u0" => u, "u1" => u, "v0" => v, "v1" => v,
    "t0" => t0, "t1" => t1, "dt" => dt, "msk" => msk, "XC" => XC, "YC" => YC)
    uv_etc=merge(uv_etc,IndividualDisplacements.NeighborTileIndices_cs(γ));

    return uv_etc
end

max_stars_repo_name: abelolivelanna/GlobalOceanNotebooks

max_issues_repo_name: abelolivelanna/GlobalOceanNotebooks

max_forks_repo_name: abelolivelanna/GlobalOceanNotebooks


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

# TODO: Fix tests.

function test_solver_multiple_dirichlet_bc()

    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]

    e1 = Quad4([1, 2, 4, 3])
    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
    e1["youngs modulus"] = 900.0
    e1["poissons ratio"] = 0.25
    e1["yield stress"] = 100.0
    e1["material model"] = :vonMises
    b1 = Seg2([3, 4])
    b1["geometry"] = Vector[N[3], N[4]]
    b1["displacement traction force"] = (
        0.0 => Vector[[0.0,    0.0], [0.0,    0.0]],
        1.0 => Vector[[0.0, -100.0], [0.0, -100.0]])

    #problem = PlaneStressElasticityProblem()
    problem = PlaneStressElasticPlasticProblem()
    push!(problem, e1)
    push!(problem, b1)

    # boundary elements for dirichlet dx=0
    dx = Seg2([1, 3])
    dx["geometry"] = Vector[N[1], N[3]]
    dx["displacement 1"] = 0.0

    # boundary elements for dirichlet dy=0
    dy = Seg2([1, 2])
    dy["geometry"] = Vector[N[1], N[2]]
    dy["displacement 2"] = 0.0

    problem2 = DirichletProblem("displacement", 2)
    push!(problem2, dx)

    problem3 = DirichletProblem("displacement", 2)
    push!(problem3, dy)

    solver = DirectSolver()
    #solver.dump_matrices = true
    solver.name = "test_solver_multiple_dirichlet_bc"
    push!(solver, problem)
    push!(solver, problem2)
    push!(solver, problem3)

    # launch solver
    #norm = solver(0.0)
    norm = solver(1.0)
    disp = e1("displacement", [1.0, 1.0], 1.0)
    info("displacement at tip: $disp")
    #@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])

end

function test_direct_cholesky_with_non_homogeneous_dirichlet_conditions()

    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]

    e1 = Quad4([1, 2, 4, 3])
    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
    e1["youngs modulus"] = 900.0
    e1["poissons ratio"] = 0.25

    problem = PlaneStressElasticityProblem()
    push!(problem, e1)

    # left boundary: dx=-0.1, dy=0.1
    bc1 = Seg2([1, 3])
    bc1["geometry"] = Vector[N[1], N[3]]
    bc1["displacement 1"] = -0.1
    bc1["displacement 2"] =  0.1

    # right boundary: dx=0.2, dy=-0.2
    bc2 = Seg2([2, 4])
    bc2["geometry"] = Vector[N[2], N[4]]
    bc2["displacement 1"] =  0.2
    bc2["displacement 2"] = -0.2

    boundary = DirichletProblem("displacement", 2)
    push!(boundary, bc1)
    push!(boundary, bc2)

    solver = DirectSolver("test_direct_cholesky_with_non_homogeneous_dirichlet_boundary_conditions")
    push!(solver, problem)
    push!(solver, boundary)

    # launch solver
    solver.method = :UMFPACK
    solver.dump_matrices = true
    solver.max_iterations = 1
    iters, status = solver(0.0)
#   FIXME: solver gives no convergence warning when all dofs are fixed.
    n1disp = e1("displacement", [-1.0, -1.0], 0.0)
    n2disp = e1("displacement", [ 1.0, -1.0], 0.0)
    n3disp = e1("displacement", [-1.0,  1.0], 0.0)
    n4disp = e1("displacement", [ 1.0,  1.0], 0.0)
    udisp = [n1disp n2disp n3disp n4disp]
    info("nodal disp = ", udisp)
    @test isapprox(n1disp, [-0.1,  0.1])
    @test isapprox(n3disp, [-0.1,  0.1])
    @test isapprox(n2disp, [ 0.2, -0.2])
    @test isapprox(n4disp, [ 0.2, -0.2])
    @test status == true
end

function test_solver_no_convergence()

    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]

    e1 = Quad4([1, 2, 4, 3])
    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
    e1["youngs modulus"] = 900.0
    e1["poissons ratio"] = 0.25
    b1 = Seg2([3, 4])
    b1["geometry"] = Vector[N[3], N[4]]
    b1["displacement traction force"] = Vector[[100.0, 100.0], [100.0, 100.0]]

    problem = PlaneStressElasticityProblem()
    push!(problem, e1)
    push!(problem, b1)

    # boundary elements for dirichlet dx=0
    dx = Seg2([1, 3])
    dx["geometry"] = Vector[N[1], N[3]]
    dx["displacement 1"] = 0.0

    # boundary elements for dirichlet dy=0
    dy = Seg2([1, 2])
    dy["geometry"] = Vector[N[1], N[2]]
    dy["displacement 2"] = 0.0

    problem2 = DirichletProblem("displacement", 2)
    push!(problem2, dx)

    problem3 = DirichletProblem("displacement", 2)
    push!(problem3, dy)

    solver = DirectSolver()
    solver.max_iterations = 1
    push!(solver, problem)
    push!(solver, problem2)
    push!(solver, problem3)

    # launch solver
    iterations, status = solver(0.0)
    @test status == false
end

function test_solver_multiple_bodies_multiple_dirichlet_bc()
    N = Vector[
        [0.0, 0.0], [1.0, 0.0],
        [0.0, 1.0], [1.0, 1.0],
        [0.0, 2.0], [1.0, 2.0]]

    e1 = Quad4([1, 2, 4, 3])
    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
    e2 = Quad4([3, 4, 6, 5])
    e2["geometry"] = Vector[N[3], N[4], N[6], N[5]]
    for el in [e1, e2]
        el["youngs modulus"] = 900.0
        el["poissons ratio"] = 0.25
    end
    b1 = Seg2([5, 6])
    b1["geometry"] = Vector[N[5], N[6]]
    b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]

    body1 = PlaneStressElasticityProblem()
    push!(body1, e1)

    body2 = PlaneStressElasticityProblem()
    push!(body2, e2)
    push!(body2, b1)

    # boundary elements for dirichlet dx=0
    dx1 = Seg2([1, 3])
    dx1["geometry"] = Vector[N[1], N[3]]
    dx2 = Seg2([3, 5])
    dx2["geometry"] = Vector[N[3], N[5]]
    for dx in [dx1, dx2]
        dx["displacement 1"] = 0.0
    end

    boundary1 = DirichletProblem("displacement", 2)
    push!(boundary1, dx1)
    push!(boundary1, dx2)

    # boundary elements for dirichlet dy=0
    dy1 = Seg2([1, 2])
    dy1["geometry"] = Vector[N[1], N[2]]
    dy1["displacement 2"] = 0.0

    boundary2 = DirichletProblem("displacement", 2)
    push!(boundary2, dy1)

    solver = DirectSolver()
    push!(solver, body1)
    push!(solver, body2)
    push!(solver, boundary1)
    push!(solver, boundary2)

    # launch solver
    norm = solver(0.0)

    disp = e2("displacement", [1.0, 1.0], 0.0)
    info("displacement at tip: $disp")
    # code aster verification, two_elements.comm
    @test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])

end


text:
using Test, YaoBlocks, BitBasis
using YaoBlocks: check_dumpload

@testset "check_dumpload" begin
    @test check_dumpload(X)
    @test check_dumpload(X + Y)
    @test check_dumpload(+(put(5, 4 => X), put(5, 5 => X)))
    @test check_dumpload(kron(put(1, 1 => X), Y))
    @test check_dumpload(kron(5, 2 => X, 4 => Y))
    @test check_dumpload(shift(0.5))
    @test check_dumpload(phase(0.5))
    @test check_dumpload(time_evolve(X, 0.5))
    @test check_dumpload(put(5, 2 => X))
    @test check_dumpload(chain(put(5, 2 => X)))
    @test check_dumpload(put(5, 2 => rot(X, 0.5)))
    @test check_dumpload(control(5, 1, 2 => rot(X, 0.5)))
    @test check_dumpload(control(5, (1, -4), 2 => rot(X, 0.5)))
    @test check_dumpload(subroutine(5, rot(SWAP, 0.5), (2, 5)))
    @test check_dumpload(repeat(5, put(1, 1 => X), (2, 5)))
    @test check_dumpload(Measure(5))
    @test check_dumpload(Measure(5, operator = put(5, 2 => X)))
    @test check_dumpload(Measure(5, locs = (3, 1), resetto = bit"01"))
    @test check_dumpload(
        Measure(5, locs = (3, 2), operator = put(2, 2 => X), resetto = bit"11"),
    )
    @test check_dumpload(Daggered(X))
    @test check_dumpload(2 * X)
    @test check_dumpload(cache(2 * X))
    @test_throws ErrorException check_dumpload(kron(5, 2:3 => SWAP))
end

@testset "yao macro" begin
    c = yao"""
    let nqubits = 5
        3 => rot(X, 0.3)
        2 => X
    end
    """
    c2 = @yaoscript let nqubits = 5
        3 => rot(X, 0.3)
        2 => X
    end

    y = chain(5, put(3 => rot(X, 0.3)), put(2 => X))
    @test c == y
    @test c == c2
    @test check_dumpload(y)
    yaotofile("_test.yao", y)
    yy = @eval $(yaofromfile("_test.yao"))
    @test y == yy

    g = eval(yaofromfile(joinpath(dirname(@__FILE__), "yaoscript.yao")))
    s = string(yaotoscript(g))
    g1 = eval(yaofromstring(s))
    @test g == g1
end


text:
module SuyamaBayes

using StatsFuns
using SpecialFunctions
using Distributions
using PyPlot
using ArgCheck
using BenchmarkTools

import Distributions: rand

export
    # pmm.jl
    # gibbs sampling
    PriorParameters,
    PoissionMixtureModel,
    init_model,
    generate_toy_data,
    sample_s,
    sample_λ,
    sample_p,
    gibbs!,
    # variational inference
    update_s,
    update_λ,
    update_p,
    vi!

include("pmm.jl")
# Write your package code here.

end

max_stars_repo_name: SosUts/SuyamaBayes

max_issues_repo_name: SosUts/SuyamaBayes

max_forks_repo_name: SosUts/SuyamaBayes


text:
struct FaceToCellGlue{A,B,C,D} <: GridapType
  face_to_bgface::A
  bgface_to_lcell::B
  face_to_cell::Vector{Int32}
  face_to_lface::Vector{Int8}
  face_to_lcell::Vector{Int8}
  face_to_ftype::C
  cell_to_ctype::D
  cell_to_lface_to_pindex::Table{Int8,Vector{Int8},Vector{Int32}}
  ctype_to_lface_to_ftype::Vector{Vector{Int8}}
end

function FaceToCellGlue(
  topo::GridTopology,
  cell_grid::Grid,
  face_grid::Grid,
  face_to_bgface::AbstractVector,
  bgface_to_lcell::AbstractVector)

  D = num_cell_dims(cell_grid)
  cD = num_cell_dims(face_grid)
  bgface_to_cells = get_faces(topo,cD,D)
  cell_to_bgfaces = get_faces(topo,D,cD)
  cell_to_lface_to_pindex = Table(get_cell_permutations(topo,cD))

  bgface_to_cell = lazy_map(getindex,bgface_to_cells, bgface_to_lcell)
  bgface_to_lface = find_local_index(bgface_to_cell, cell_to_bgfaces)

  face_to_cell = collect(Int32,lazy_map(Reindex(bgface_to_cell), face_to_bgface))
  face_to_lface = collect(Int8,lazy_map(Reindex(bgface_to_lface), face_to_bgface))
  face_to_lcell = collect(Int8,lazy_map(Reindex(bgface_to_lcell), face_to_bgface))

  f = (p)->fill(Int8(UNSET),num_faces(p,cD))
  ctype_to_lface_to_ftype = map( f, get_reffes(cell_grid) )
  face_to_ftype = get_cell_type(face_grid)
  cell_to_ctype = get_cell_type(cell_grid)

  _fill_ctype_to_lface_to_ftype!(
    ctype_to_lface_to_ftype,
    face_to_cell,
    face_to_lface,
    face_to_ftype,
    cell_to_ctype)

  FaceToCellGlue(
    face_to_bgface,
    bgface_to_lcell,
    face_to_cell,
    face_to_lface,
    face_to_lcell,
    face_to_ftype,
    cell_to_ctype,
    cell_to_lface_to_pindex,
    ctype_to_lface_to_ftype)
end

function _fill_ctype_to_lface_to_ftype!(
  ctype_to_lface_to_ftype,
  face_to_cell,
  face_to_lface,
  face_to_ftype,
  cell_to_ctype)
  for (face, cell) in enumerate(face_to_cell)
    ctype = cell_to_ctype[cell]
    lface = face_to_lface[face]
    ftype = face_to_ftype[face]
    ctype_to_lface_to_ftype[ctype][lface] = ftype
  end
end

function get_children(n::TreeNode, a::FaceToCellGlue)
  (
    similar_tree_node(n,a.face_to_cell),
    similar_tree_node(n,a.face_to_lface),
    similar_tree_node(n,a.face_to_lcell),
    similar_tree_node(n,a.face_to_ftype),
    similar_tree_node(n,a.cell_to_ctype),
    similar_tree_node(n,a.cell_to_lface_to_pindex),
    similar_tree_node(n,a.ctype_to_lface_to_ftype)
  )
end

"""
"""
struct BoundaryTriangulation{Dc,Dp,A,B} <: Triangulation{Dc,Dp}
  trian::A
  glue::B

  function BoundaryTriangulation(
    trian::BodyFittedTriangulation,
    glue::FaceToCellGlue)

    Dc = num_cell_dims(trian)
    Dp = num_point_dims(trian)
    A = typeof(trian)
    B = typeof(glue)
    new{Dc,Dp,A,B}(trian,glue)
  end
end

function Boundary(args...;kwargs...)
  BoundaryTriangulation(args...;kwargs...)
end

# Constructors

"""
    BoundaryTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
    BoundaryTriangulation(model::DiscreteModel)
"""
function BoundaryTriangulation(
  model::DiscreteModel,
  face_to_bgface::AbstractVector{<:Integer},
  bgface_to_lcell::AbstractVector{<:Integer})

  D = num_cell_dims(model)
  topo = get_grid_topology(model)
  bgface_grid = Grid(ReferenceFE{D-1},model)

  face_grid = view(bgface_grid,face_to_bgface)
  cell_grid = get_grid(model)
  glue = FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell)
  trian = BodyFittedTriangulation(model,face_grid,face_to_bgface)

  BoundaryTriangulation(trian,glue)
end

function BoundaryTriangulation(
  model::DiscreteModel,
  face_to_bgface::AbstractVector{<:Integer})
  BoundaryTriangulation(model,face_to_bgface,Fill(1,num_facets(model)))
end

function BoundaryTriangulation(
  model::DiscreteModel,
  bgface_to_mask::AbstractVector{Bool},
  bgface_to_lcell::AbstractVector{<:Integer})

  face_to_bgface = findall(bgface_to_mask)
  BoundaryTriangulation(model,face_to_bgface,bgface_to_lcell)
end

function BoundaryTriangulation(
  model::DiscreteModel, bgface_to_mask::AbstractVector{Bool}, lcell::Integer=1)
  BoundaryTriangulation(model,bgface_to_mask, Fill(lcell,num_facets(model)) )
end

"""
    BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling;tags::Vector{Int})
    BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling;tags::Vector{String})
    BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling;tag::Int)
    BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling;tag::String)
"""
function BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling;tags=nothing)
  D = num_cell_dims(model)
  if tags == nothing
    topo = get_grid_topology(model)
    face_to_mask = get_isboundary_face(topo,D-1)
  else
    face_to_mask = get_face_mask(labeling,tags,D-1)
  end
  BoundaryTriangulation(model,face_to_mask)
end

"""
    BoundaryTriangulation(model::DiscreteModel,tags::Vector{Int})
    BoundaryTriangulation(model::DiscreteModel,tags::Vector{String})
    BoundaryTriangulation(model::DiscreteModel,tag::Int)
    BoundaryTriangulation(model::DiscreteModel,tag::String)
"""
function BoundaryTriangulation(model::DiscreteModel;tags=nothing)
  labeling = get_face_labeling(model)
  BoundaryTriangulation(model,labeling,tags=tags)
end

function BoundaryTriangulation(rtrian::Triangulation,args...;kwargs...)
  rmodel = get_active_model(rtrian)
  dtrian = BoundaryTriangulation(rmodel,args...;kwargs...)
  CompositeTriangulation(rtrian,dtrian)
end

# API

get_background_model(t::BoundaryTriangulation) = get_background_model(t.trian)
get_grid(t::BoundaryTriangulation) = get_grid(t.trian)
get_glue(t::BoundaryTriangulation{D},::Val{D}) where D = get_glue(t.trian,Val(D))

function get_glue(trian::BoundaryTriangulation,::Val{Dp}) where Dp
  model = get_background_model(trian)
  Dm = num_cell_dims(model)
  get_glue(trian,Val(Dp),Val(Dm))
end

function get_glue(trian::BoundaryTriangulation,::Val{Dp},::Val{Dm}) where {Dp,Dm}
  nothing
end

function get_glue(trian::BoundaryTriangulation,::Val{D},::Val{D}) where D
  tface_to_mface = trian.glue.face_to_cell
  face_to_q_vertex_coords = _compute_face_to_q_vertex_coords(trian)
  f(p) = get_shapefuns(LagrangianRefFE(Float64,get_polytope(p),1))
  ftype_to_shapefuns = map( f, get_reffes(trian) )
  face_to_shapefuns = expand_cell_data(ftype_to_shapefuns,trian.glue.face_to_ftype)
  face_s_q = lazy_map(linear_combination,face_to_q_vertex_coords,face_to_shapefuns)
  tface_to_mface_map = face_s_q
  mface_to_tface = nothing
  FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
end

function get_facet_normal(trian::BoundaryTriangulation)

  glue = trian.glue
  cell_grid = get_grid(get_background_model(trian.trian))

  ## Reference normal
  function f(r)
    p = get_polytope(r)
    lface_to_n = get_facet_normal(p)
    lface_to_pindex_to_perm = get_face_vertex_permutations(p,num_cell_dims(p)-1)
    nlfaces = length(lface_to_n)
    lface_pindex_to_n = [ fill(lface_to_n[lface],length(lface_to_pindex_to_perm[lface])) for lface in 1:nlfaces ]
    lface_pindex_to_n
  end
  ctype_lface_pindex_to_nref = map(f, get_reffes(cell_grid))
  face_to_nref = FaceCompressedVector(ctype_lface_pindex_to_nref,glue)
  face_s_nref = lazy_map(constant_field,face_to_nref)

  # Inverse of the Jacobian transpose
  cell_q_x = get_cell_map(cell_grid)
  cell_q_Jt = lazy_map(∇,cell_q_x)
  cell_q_invJt = lazy_map(Operation(pinvJt),cell_q_Jt)
  face_q_invJt = lazy_map(Reindex(cell_q_invJt),glue.face_to_cell)

  # Change of domain
  D = num_cell_dims(cell_grid)
  glue = get_glue(trian,Val(D))
  face_s_q = glue.tface_to_mface_map
  face_s_invJt = lazy_map(∘,face_q_invJt,face_s_q)
  face_s_n = lazy_map(Broadcasting(Operation(push_normal)),face_s_invJt,face_s_nref)
  Fields.MemoArray(face_s_n)
end

@inline function push_normal(invJt,n)
  v = invJt⋅n
  m = sqrt(inner(v,v))
  if m < eps()
    return zero(n)
  else
    return v/m
  end
end

function _compute_face_to_q_vertex_coords(trian::BoundaryTriangulation)
    d = num_cell_dims(trian)
    cell_grid = get_grid(get_background_model(trian.trian))
    polytopes = map(get_polytope, get_reffes(cell_grid))
    cell_to_ctype = trian.glue.cell_to_ctype
    ctype_to_lvertex_to_qcoords = map(get_vertex_coordinates, polytopes)
    ctype_to_lface_to_lvertices = map((p)->get_faces(p,d,0), polytopes)
    ctype_to_lface_to_pindex_to_perm = map( (p)->get_face_vertex_permutations(p,d), polytopes)

    P = eltype(eltype(ctype_to_lvertex_to_qcoords))
    D = num_components(P)
    T = eltype(P)
    ctype_to_lface_to_pindex_to_qcoords = Vector{Vector{Vector{Point{D,T}}}}[]

    for (ctype, lface_to_pindex_to_perm) in enumerate(ctype_to_lface_to_pindex_to_perm)
      lvertex_to_qcoods = ctype_to_lvertex_to_qcoords[ctype]
      lface_to_pindex_to_qcoords = Vector{Vector{Point{D,T}}}[]
      for (lface, pindex_to_perm) in enumerate(lface_to_pindex_to_perm)
        cfvertex_to_lvertex = ctype_to_lface_to_lvertices[ctype][lface]
        nfvertices = length(cfvertex_to_lvertex)
        pindex_to_qcoords = Vector{Vector{Point{D,T}}}(undef,length(pindex_to_perm))
        for (pindex, cfvertex_to_ffvertex) in enumerate(pindex_to_perm)
          ffvertex_to_qcoords = zeros(Point{D,T},nfvertices)
          for (cfvertex, ffvertex) in enumerate(cfvertex_to_ffvertex)
            lvertex = cfvertex_to_lvertex[cfvertex]
            qcoords = lvertex_to_qcoods[lvertex]
            ffvertex_to_qcoords[ffvertex] = qcoords
          end
          pindex_to_qcoords[pindex] = ffvertex_to_qcoords
        end
        push!(lface_to_pindex_to_qcoords,pindex_to_qcoords)
      end
      push!(ctype_to_lface_to_pindex_to_qcoords,lface_to_pindex_to_qcoords)
    end

    FaceCompressedVector(ctype_to_lface_to_pindex_to_qcoords,trian.glue)
end

struct FaceCompressedVector{T,G<:FaceToCellGlue} <: AbstractVector{T}
  ctype_lface_pindex_to_value::Vector{Vector{Vector{T}}}
  glue::G
end

Base.size(a::FaceCompressedVector) = (length(a.glue.face_to_cell),)

Base.IndexStyle(::Type{<:FaceCompressedVector}) = IndexLinear()

function Base.getindex(a::FaceCompressedVector,face::Integer)
  cell = a.glue.face_to_cell[face]
  ctype = a.glue.cell_to_ctype[cell]
  lface = a.glue.face_to_lface[face]
  p = a.glue.cell_to_lface_to_pindex.ptrs[cell]-1
  pindex = a.glue.cell_to_lface_to_pindex.data[p+lface]
  value = a.ctype_lface_pindex_to_value[ctype][lface][pindex]
  value
end

function get_children(n::TreeNode, a::FaceCompressedVector)
  (similar_tree_node(n,a.ctype_lface_pindex_to_value),similar_tree_node(n,a.glue))
end

function lazy_map(k::Fields.LinearCombinationMap,::Type{T},b::FaceCompressedVector,c::Fill) where T
  d = CompressedArray([c.value,],Fill(1,length(c)))
  lazy_map(k,T,a,b,d)
end

function lazy_map(k::Fields.LinearCombinationMap,::Type{T},b::FaceCompressedVector,c::CompressedArray) where T
  if c.ptrs === b.glue.face_to_ftype || c.ptrs == b.glue.face_to_ftype

    ctype_lface_pindex_to_r = Vector{Vector{Vector{T}}}(undef,length(b.ctype_lface_pindex_to_value))
    for (ctype, lface_pindex_to_value) in enumerate(b.ctype_lface_pindex_to_value)
      lface_pindex_to_r = Vector{Vector{T}}(undef,length(lface_pindex_to_value))
      for (lface, pindex_to_value) in enumerate(lface_pindex_to_value)
        pindex_to_r = Vector{T}(undef,length(pindex_to_value))
        ftype = b.glue.ctype_to_lface_to_ftype[ctype][lface]
        if ftype != UNSET
          for (pindex, value) in enumerate(pindex_to_value)
            pindex_to_r[pindex] = k(value,c.values[ftype])
          end
        end
        lface_pindex_to_r[lface] = pindex_to_r
      end
      ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
    end
    FaceCompressedVector(ctype_lface_pindex_to_r,b.glue)

  else
    @notimplemented
  end
end

function lazy_map(k::typeof(evaluate),::Type{T},a::Fill,b::FaceCompressedVector) where T
  @check length(a) == length(b)
  ctype_lface_pindex_to_r = Vector{Vector{Vector{T}}}(undef,length(b.ctype_lface_pindex_to_value))
  for (ctype, lface_pindex_to_value) in enumerate(b.ctype_lface_pindex_to_value)
    lface_pindex_to_r = Vector{Vector{T}}(undef,length(lface_pindex_to_value))
    for (lface, pindex_to_value) in enumerate(lface_pindex_to_value)
      pindex_to_r = Vector{T}(undef,length(pindex_to_value))
      ftype = b.glue.ctype_to_lface_to_ftype[ctype][lface]
      if ftype != UNSET
        for (pindex, value) in enumerate(pindex_to_value)
          pindex_to_r[pindex] = evaluate(a.value,value)
        end
      end
      lface_pindex_to_r[lface] = pindex_to_r
    end
    ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
  end
  FaceCompressedVector(ctype_lface_pindex_to_r,b.glue)
end

function lazy_map(
  k::typeof(evaluate),
  ::Type{T},
  a::CompressedArray,
  b::FaceCompressedVector) where T

  @check length(a) == length(b)
  ctype_lface_pindex_to_r = Vector{Vector{Vector{T}}}(undef,length(b.ctype_lface_pindex_to_value))
  for (ctype, lface_pindex_to_value) in enumerate(b.ctype_lface_pindex_to_value)
    lface_pindex_to_r = Vector{Vector{T}}(undef,length(lface_pindex_to_value))
    for (lface, pindex_to_value) in enumerate(lface_pindex_to_value)
      pindex_to_r = Vector{T}(undef,length(pindex_to_value))
      ftype = b.glue.ctype_to_lface_to_ftype[ctype][lface]
      if ftype != UNSET
        for (pindex, value) in enumerate(pindex_to_value)
          pindex_to_r[pindex] = evaluate(a.values[ctype],value)
        end
      end
      lface_pindex_to_r[lface] = pindex_to_r
    end
    ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
  end
  FaceCompressedVector(ctype_lface_pindex_to_r,b.glue)
end

function lazy_map(k::typeof(evaluate),::Type{T},b::FaceCompressedVector,a::Fill) where T
  @check length(a) == length(b)
  ctype_lface_pindex_to_r = Vector{Vector{Vector{T}}}(undef,length(b.ctype_lface_pindex_to_value))
  for (ctype, lface_pindex_to_value) in enumerate(b.ctype_lface_pindex_to_value)
    lface_pindex_to_r = Vector{Vector{T}}(undef,length(lface_pindex_to_value))
    for (lface, pindex_to_value) in enumerate(lface_pindex_to_value)
      pindex_to_r = Vector{T}(undef,length(pindex_to_value))
      ftype = b.glue.ctype_to_lface_to_ftype[ctype][lface]
      if ftype != UNSET
        for (pindex, value) in enumerate(pindex_to_value)
          pindex_to_r[pindex] = evaluate(value,a.value)
        end
      end
      lface_pindex_to_r[lface] = pindex_to_r
    end
    ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
  end
  FaceCompressedVector(ctype_lface_pindex_to_r,b.glue)
end

function lazy_map(k::typeof(evaluate),::Type{T},a::Fill,b::FaceCompressedVector,c::FaceCompressedVector) where T
  if b.glue !== c.glue
    return LazyArray(T,a,b,c)
  end
  @check length(a) == length(b)
  ctype_lface_pindex_to_r = Vector{Vector{Vector{T}}}(undef,length(b.ctype_lface_pindex_to_value))
  for (ctype, lface_pindex_to_value) in enumerate(b.ctype_lface_pindex_to_value)
    lface_pindex_to_r = Vector{Vector{T}}(undef,length(lface_pindex_to_value))
    for (lface, pindex_to_value) in enumerate(lface_pindex_to_value)
      pindex_to_r = Vector{T}(undef,length(pindex_to_value))
      ftype = b.glue.ctype_to_lface_to_ftype[ctype][lface]
      if ftype != UNSET
        for (pindex, bvalue) in enumerate(pindex_to_value)
          cvalue = c.ctype_lface_pindex_to_value[ctype][lface][pindex]
          pindex_to_r[pindex] = evaluate(a.value,bvalue,cvalue)
        end
      end
      lface_pindex_to_r[lface] = pindex_to_r
    end
    ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
  end
  FaceCompressedVector(ctype_lface_pindex_to_r,b.glue)
end


text:
# --------------------------------------------------
# HYPERDIMENSIONAL COMPUTING IN JULIA: HOST CLASSIFICATION
#
# @author: dimiboeckaerts
# --------------------------------------------------

# LIBRARIES
# --------------------------------------------------
using CSV
using MLJ
using DataFrames
using ProgressMeter
using ScikitLearn.CrossValidation: KFold
include("hyperdimensional_computing.jl")

# DATA PROCESSING
# --------------------------------------------------
fibers = DataFrame(CSV.File("/Users/Dimi/Documents/GitHub_Local/BacteriophageHostPrediction/RBP_database887.csv"))
sequences = fibers.protein_sequence;
classes = fibers.class;
classes = coerce(classes, Multiclass);

# HYPERDIMENSIONAL COMPUTING
# --------------------------------------------------
#dimensions = LinRange(5000, 10000, 6)
dimension = 10000
nfolds = 5
kvalues = LinRange(5, 10, 6)
encoded_alphabet = encode_items(["A", "G", "I", "L", "P", "V", "F", "W", "Y", "D",
                                    "E", "R", "H", "K", "S", "T", "C", "M", "N", "Q"],
                                    dim=dimension)

# kfold over the sequences: use 80% to build the class HVs and 20% to make predictions for
for kval in kvalues
    predictions = zeros(length(sequences))
    pb = Progress(nfolds)
    for (train_indices, test_indices) in KFold(length(sequences), n_folds=nfolds)
        # define matrices for hypervectors & labels
        x_train = zeros(length(train_indices), dimension)
        x_test = zeros(length(test_indices), dimension)
        y_train = classes[train_indices]; y_test = classes[test_indices]

        # encode sequences
        for (i, train_index) in enumerate(train_indices)
            x_train[i,:] = encode_sequence(sequences[train_index], encoded_alphabet, k=Int64.(round(kval, digits=1)))
        end
        for (i, test_index) in enumerate(test_indices)
            x_test[i,:] = encode_sequence(sequences[test_index], encoded_alphabet, k=Int64.(round(kval, digits=1)))
        end

        # encode class vectors
        class_vectors = encode_associative_memory(x_train, y_train)

        # make predictions
        preds = make_predictions(x_test, class_vectors)

        # save predictions
        for (i, test_index) in enumerate(test_indices)
            predictions[test_index] = preds[i]
        end
        next!(pb)
    end

    # compute performance
    predictions = Int64.(predictions)
    predictions = coerce(predictions, Multiclass)
    F1_it = MLJ.macro_f1score(predictions, classes)
    println("F1 score for k-mer value of ", kval, " : ", F1_it)
end

# compute performance
predictions = Int64.(predictions)
predictions = coerce(predictions, Multiclass)
MLJ.macro_f1score(predictions, classes)
MLJ.confmat(predictions, classes)

re = MLJ.multiclass_recall(predictions, classes)
pre = MLJ.multiclass_positive_predictive_value(predictions, classes)
Fm = (2*pre*re)/(pre+re)

# compute weighted multiclass F1 score
class_weights = [count(i->(i==x), classes) for x in levels(classes)]
total_F1 = 0
for (i, level) in enumerate(levels(classes))
    # binarize arrays
    binary_predictions = coerce(Int64.(predictions .== level), Multiclass)
    binary_classes = coerce(Int64.(classes .== level), Multiclass)
    # compute binary F1
    F1 = MLJ.f1score(binary_predictions, binary_classes)
    println("F1 score for class ", level, " : ", F1)
    # refactor for class weights
    total_F1 += F1*(class_weights[i]/sum(class_weights))
end
println("Weighted average F1 score: ", total_F1) # 80.05%