text:
@testset "link" begin
    x = rand()
    xs = rand(3)

    # Generic link
    f = sin
    l = Link(f)
    @test l(x) == f(x)

    # Log
    l = LogLink()
    @test l(x) == log(x)
    @test inv(l) == ExpLink()
    @test inv(inv(l)) == l

    # Exp
    l = ExpLink()
    @test l(x) == exp(x)
    @test inv(l) == LogLink()
    @test inv(inv(l)) == l

    # Sqrt
    l = SqrtLink()
    @test l(x) == sqrt(x)
    @test inv(l) == SquareLink()
    @test inv(inv(l)) == l

    # Square
    l = SquareLink()
    @test l(x) == x^2
    @test inv(l) == SqrtLink()
    @test inv(inv(l)) == l

    # Logit
    l = LogitLink()
    @test l(x) == logit(x)
    @test inv(l) isa LogisticLink
    @test inv(inv(l)) == l

    # Logistic
    l = LogisticLink()
    @test l(x) == logistic(x)
    @test inv(l) isa LogitLink
    @test inv(inv(l)) == l

    # Probit
    l = ProbitLink()
    @test l(x) == norminvcdf(x)
    @test inv(l) == NormalCDFLink()
    @test inv(inv(l)) == l

    # NormalCDF
    l = NormalCDFLink()
    @test l(x) == normcdf(x)
    @test inv(l) == ProbitLink()
    @test inv(inv(l)) == l

    # SoftMax
    l = SoftMaxLink()
    @test l(xs) == softmax(xs)
end


text:
@with_kw struct CollisionAvoidance
    ddh_max::Float64 = 1.0 # vertical acceleration limit [m/s²]
    collision_threshold::Float64 = 50.0 # collision threshold [m]
    reward_collision::Float64 = -1.0 # reward obtained if collision occurs
    reward_change::Float64 = -0.01 # reward obtained if action changes
    𝒜::Vector{Float64} = [-5.0, 0.0, 5.0] # Available actions to command for vertical rate [m/s]
    # The transition distribution produces successor states
    # with dh noise with the below probabilities.
    Pν::SetCategorical{Float64} = SetCategorical([2.0, 0.0, -2.0], [0.25, 0.5, 0.25])
end

struct CollisionAvoidanceState
    h::Float64 # vertical separation [m]
    dh::Float64 # rate of change of h [m/s]
    a_prev::Float64 # previous acceleration [m/s²]
    τ::Float64 # horizontal separation time [s]
end

Base.vec(s::CollisionAvoidanceState) = [s.h, s.dh, s.a_prev, s.τ]

function transition(𝒫::CollisionAvoidance, s::CollisionAvoidanceState, a::Float64)
    h = s.h + s.dh
    dh = s.dh
    if a != 0.0
        if abs(a - dh) < 𝒫.ddh_max
            dh += a
        else
            dh += sign(a - dh)*𝒫.ddh_max
        end
    end
    a_prev = a
    τ = max(s.τ - 1.0, -1.0)
    states = [
        CollisionAvoidanceState(h, dh + ν, a_prev, τ) for ν in 𝒫.Pν.elements
    ]
    return SetCategorical(states, 𝒫.Pν.distr.p)
end

is_terminal(𝒫::CollisionAvoidance, s::CollisionAvoidanceState) = s.τ < 0.0

function reward(𝒫::CollisionAvoidance, s::CollisionAvoidanceState, a::Float64)
    r = 0.0
    if abs(s.h) < 𝒫.collision_threshold && abs(s.τ) < eps()
        # We collided
        r += 𝒫.reward_collision
    end
    if a != s.a_prev
        # We changed our action
        r += 𝒫.reward_change
    end
    return r
end

@with_kw struct CollisionAvoidanceStateDistribution
    a_prev = 0.0
    h = Uniform(-200, 200)
    dh = Uniform(-10, 10)
    tau = 40.0
end

function rand(b::CollisionAvoidanceStateDistribution)
    return CollisionAvoidanceState(Distributions.rand(b.h), Distributions.rand(b.dh), b.a_prev, b.tau)
end

@with_kw struct SimpleCollisionAvoidancePolicy
    thresh_h = 50
    thresh_τ = 30
    𝒜 = (up = 5.0, down = -5.0, noalert = 0.0)
end

struct OptimalCollisionAvoidancePolicy
    𝒜
    grid
    Q
end

function OptimalCollisionAvoidancePolicy(mdp = CollisionAvoidance())
    𝒜 = mdp.𝒜

    hs = range(-200, 200, length=21) # discretization of h in m
    dhs = range(-10, 10, length=21) # discretization of dh in m/s
    τs = range(0, 40, length=41) # discretization of τ in s

    # Discretization grid
    grid = GridInterpolations.RectangleGrid(hs, dhs, 𝒜, τs)

    # State space
    𝒮 = [CollisionAvoidanceState(h, dh, a_prev, τ) for h in hs, dh in dhs, a_prev in mdp.𝒜, τ in τs]

    # State value function
    U = zeros(length(𝒮))

    # State-action value function
    Q = [zeros(length(𝒮)) for a in mdp.𝒜]

    # Solve with backwards induction value iteration
    for (si, s) in enumerate(𝒮)
        for (ai, a) in enumerate(mdp.𝒜)
            Tsa = transition(mdp, s, a)
            Q[ai][si] = reward(mdp, s, a)
            Q[ai][si] += sum(is_terminal(mdp, s′) ? 0.0 : Tsa.distr.p[j]*GridInterpolations.interpolate(grid, U, vec(s′)) for (j, s′) in enumerate(Tsa.elements))
        end
        U[si] = maximum(q[si] for q in Q)
    end
    return OptimalCollisionAvoidancePolicy(mdp.𝒜, grid, Q)
end

function action(policy::OptimalCollisionAvoidancePolicy, s::CollisionAvoidanceState)
    vec_s = vec(s)
    a_best = first(policy.𝒜)
    q_best = -Inf
    for (a,q) in zip(policy.𝒜, policy.Q)
        q = GridInterpolations.interpolate(policy.grid, q, vec_s)
        if q > q_best
            a_best, q_best = a, q
        end
    end
    return a_best
end

function (policy::OptimalCollisionAvoidancePolicy)(s::CollisionAvoidanceState)
    return action(policy, s)
end

function action(policy::SimpleCollisionAvoidancePolicy, s::CollisionAvoidanceState)
    if abs(s.h) < policy.thresh_h && s.τ < policy.thresh_τ
        return (s.h > 0.0) ? policy.𝒜.up : policy.𝒜.down
    end
    return policy.𝒜.noalert
end

function (policy::SimpleCollisionAvoidancePolicy)(s::CollisionAvoidanceState)
    return action(policy, s)
end

struct CollisionAvoidanceValueFunction
    grid
    U
end

function CollisionAvoidanceValueFunction(𝒫::CollisionAvoidance, policy)
    𝒜 = 𝒫.𝒜

    hs = range(-200, 200, length=21) # discretization of h in m
    dhs = range(-10, 10, length=21) # discretization of dh in m/s
    τs = range(0, 40, length=41) # discretization of τ in s

    # Discretization grid
    grid = GridInterpolations.RectangleGrid(hs, dhs, 𝒜, τs)

    # State space
    𝒮 = [CollisionAvoidanceState(h, dh, a_prev, τ) for h in hs, dh in dhs, a_prev in 𝒫.𝒜, τ in τs]

    # State value function
    U = zeros(length(𝒮))

    # Solve with backwards induction value iteration
    for (si, s) in enumerate(𝒮)
        a = action(policy, s)
        Tsa = transition(𝒫, s, a)
        q = reward(𝒫, s, a)
        q += sum(is_terminal(𝒫, s′) ? 0.0 : Tsa.distr.p[j]*GridInterpolations.interpolate(grid, U, vec(s′)) for (j, s′) in enumerate(Tsa.elements))
        U[si] = q
    end
    return CollisionAvoidanceValueFunction(grid, U)
end

function (U::CollisionAvoidanceValueFunction)(s)
    return GridInterpolations.interpolate(U.grid, U.U, vec(s))
end

function MDP(mdp::CollisionAvoidance; γ::Float64=1.0)
    return MDP(
            γ,
            nothing, # no ordered states
            mdp.𝒜,
            (s,a) -> transition(mdp,s,a), # no probabilistic transition function
            (s,a) -> reward(mdp, s, a),
            (s,a)->begin
                s′ = rand(transition(mdp,s,a))
                r = reward(mdp, s, a)
                return (s′, r)
            end
        )
end

ext: jl


text:
abstract type AbstractActivation end

struct Sigmoid <: AbstractActivation end
struct SoftPlus <: AbstractActivation end
struct Relu <: AbstractActivation end
struct Tanh <: AbstractActivation end
struct Identity <: AbstractActivation end

activate(::Sigmoid, x::Number) = 1 ./ (1 + exp(-x))
activate(::SoftPlus, x::Number) = log(1 + exp(x))
activate(::Relu, x::Number) = max(x, 0)
activate(::Tanh, x::Number) = tanh(x)
activate(::Identity, x::Number) = x

activate(a::AbstractActivation, x::AbstractArray) =
    map(_1->activate(a, _1), x)

activate!(out::AbstractArray, a::AbstractActivation, x::AbstractArray) =
    map!(_1->activate(a, _1), out, x)

for T in InteractiveUtils.subtypes(AbstractActivation)
    @eval (ta::$(T))(x) where T = activate(ta, x)
end

deriv(::Sigmoid, x) = exp(-x) ./((1 + exp(-x)).^2)
deriv(::SoftPlus, x) = exp(x) ./ (1 + exp(x))
deriv(::Relu, x::Number) = x >= 0.0 ? 1.0 : 0.0
deriv(r::Relu, x::AbstractArray) = map(_1->deriv(r, _1), x)
deriv(::Tanh, x) = 1 - tanh(x).^2
deriv(::Identity, x::Number) = one(x)
deriv(::Identity, x::AbstractArray) = ones(x)


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# DynamicBoundspODEIneq.jl
# See https://github.com/PSORLab/DynamicBoundspODEIneq.jl
#############################################################################
# src/DynamicBoundspODEsIneq.jl
# main module file
#############################################################################

module DynamicBoundspODEsIneq

using DynamicBoundsBase, McCormick, ElasticArrays, DocStringExtensions
using ForwardDiff: Chunk, Dual, Partials, construct_seeds, single_seed
using DiffEqSensitivity: extract_local_sensitivities, ODEForwardSensitivityProblem
using DiffEqBase: remake, AbstractODEProblem, AbstractContinuousCallback, solve
using Sundials

const DBB = DynamicBoundsBase

import DynamicBoundsBase: relax!, set!, setall!, get, getall!, relax!
export DifferentialInequality, set!, setall!, get, getall!, relax!,
       DifferentialInequalityCond, DifferentialInequalityAffect,
       DifferentialInequalityAffectNeg, DifferentialInequalityf,
       LocalProblemStorage

import Base.size

include("contractors.jl")
include("DifferentialInequality.jl")

end


text:
##################################################################################################
## Transfer matrix method functions for conditional expectations E[⋅|y]
##################################################################################################
"""

    `get_uc(α,call)`

    Generates vector u1 and uN used in transfer matrix method for p(x̄|y).

    # Examples
    ```julia-repl
    julia> call = CpelNano.MethCallCpgGrp(-50.0,-45.0);
    julia> CpelNano.get_uc(0.0,call)
    2-element Array{Float64,1}:
     1.3887943864964021e-11
     1.6918979226151304e-10
    ```
"""
function get_uc(α::Float64,call::MethCallCpgGrp)::Vector{Float64}

    # Return u for CpG passed
    return [exp(get_γc(false,α,call)); exp(get_γc(true,α,call))]

end
"""

    `get_Wc(α1,α2,β,call1,call2)`

    Generates transition matrix W used in transfer matrix method for p(x̄|y).

    # Examples
    ```julia-repl
    julia> call = CpelNano.MethCallCpgGrp(-50.0,-45.0);
    julia> Wc = CpelNano.get_Wc(0.0,0.0,0.0,call,call)
    2×2 Array{Float64,2}:
     1.92875e-22  2.3497e-21
     2.3497e-21   2.86252e-20
    ```
"""
function get_Wc(α1::Float64,α2::Float64,β::Float64,call1::MethCallCpgGrp,call2::MethCallCpgGrp)::Array{Float64,2}

    # Store intermediate quantities
    γ1_f = get_γc(false,α1,call1)
    γ1_t = get_γc(true,α1,call1)
    γ2_f = get_γc(false,α2,call2)
    γ2_t = get_γc(true,α2,call2)
    exp_γ1_f = exp(γ1_f)
    exp_γ1_t = exp(γ1_t)
    exp_γ2_f = exp(γ2_f)
    exp_γ2_t = exp(γ2_t)
    exp_β = exp(β)

    # Define transition matrix
    W = zeros(2,2)
    W[1,1] = exp_γ1_f * exp_γ2_f * exp_β
    W[2,1] = exp_γ1_t * exp_γ2_f / exp_β
    W[1,2] = exp_γ1_f * exp_γ2_t / exp_β
    W[2,2] = exp_γ1_t * exp_γ2_t * exp_β

    # Return W
    return W

end
"""

    `get_Zc(u1,uL,Ws)`

    Computes partition function of p(x̄_m|y_m;ϕ) using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 2; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-100.0,-100.0),L);
    julia> u1 = CpelNano.get_uc(αs[1],obs[1]); uL = CpelNano.get_uc(αs[end],obs[end]);
    julia> Ws = [CpelNano.get_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> CpelNano.get_Zc(u1,uL,Ws)
    5.53558610694695e-87
    ```
"""
function get_Zc(u1::Vector{Float64},uL::Vector{Float64},Ws::Vector{Array{Float64,2}})::Float64

    # Return u_1⊺ W_1 ⋯ W_L u_L for m-th observation
    return u1' * prod(Ws) * uL

end
"""

    `get_Ec_X(u1,uL,Ws)`

    Computes E[X̄|y_m;ϕ] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 2; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-175.0,-170.0),L);
    julia> u1 = CpelNano.get_uc(αs[1],obs[1]); uL = CpelNano.get_uc(αs[end],obs[end]);
    julia> Ws = [CpelNano.get_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> Zc = CpelNano.get_Zc(u1,uL,Ws);
    julia> CpelNano.get_Ec_X(u1,uL,Ws,Zc)
    2-element Array{Float64,1}:
     0.9866142981514304
     0.9866142981514304
    ```
"""
function get_Ec_X(u1::Vector{Float64},uL::Vector{Float64},Ws::Vector{Array{Float64,2}},Zc::Float64)::Vector{Float64}

    # Compute E[X̄_l|y_m] for l=1,2,…,L
    prod_Ws = prod(Ws)
    exs = fill(NaN,length(Ws)+1)
    @inbounds for l=1:length(exs)
        if l==1
            exs[1] = u1' * D_EX * prod_Ws * uL
        elseif l==length(exs)
            exs[end] = u1' * prod_Ws * D_EX * uL
        else
            exs[l] = u1' * prod(Ws[1:(l-1)]) * D_EX * prod(Ws[l:end]) * uL
        end
    end

    # Return vector of Ec[X̄_l|y_m] for m-th observation
    return exs ./ Zc

end
"""

    `get_Ec_XX(u1,uL,Ws)`

    Computes E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-50.0,-50.0),L);
    julia> u1 = CpelNano.get_uc(αs[1],obs[1]); uL = CpelNano.get_uc(αs[end],obs[end]);
    julia> Ws = [CpelNano.get_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> Zc = CpelNano.get_Zc(u1,uL,Ws);
    julia> CpelNano.get_Ec_XX(u1,uL,Ws,Zc)
    9-element Array{Float64,1}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
    ```
"""
function get_Ec_XX(u1::Vector{Float64},uL::Vector{Float64},Ws::Vector{Array{Float64,2}},Zc::Float64)::Vector{Float64}

    # Init output
    exxs = fill(NaN,length(Ws))

    # If only one W
    if length(Ws)==1
        exxs[1] = u1'*(D_EXX.* Ws[1]) * uL
        return exxs ./ Zc
    end

    # Compute E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] for l=1,2,…,L-1
    @inbounds for l=1:length(exxs)
        if l==1
            exxs[1] = u1' * (D_EXX .* Ws[1]) * prod(Ws[2:end]) * uL
        elseif l==length(Ws)
            exxs[end] = u1' * prod(Ws[1:(end-1)]) * (Ws[end] .* D_EXX) * uL
        else
            exxs[l] = u1' * prod(Ws[1:(l-1)]) * (Ws[l] .* D_EXX) * prod(Ws[(l+1):end]) * uL
        end
    end

    # Return vector of E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] for m-th observation
    return exxs ./ Zc

end
"""

    `get_cond_exs(αs,βs,obs)`

    Computes E[X̄_{m,l}|Y_m=y_m;ϕ] & E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] using the transfer matrix method. This function
    does not use the log-sum trick, however. Thus, it is not as numerically stable as `get_cond_exs_log()`.

    # Examples
    ```julia-repl
    julia> L = 5; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-50.0,-45.0),L);
    julia> ex,exx = CpelNano.get_cond_exs(αs,βs,obs);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-45.0,-50.0),L);
    julia> ex,exx = CpelNano.get_cond_exs(αs,βs,obs);
    ```
"""
function get_cond_exs(αs::Vector{Float64},βs::Vector{Float64},obs::Vector{MethCallCpgGrp})::NTuple{2,Vector{Float64}}

    # Get u's
    u1 = get_uc(αs[1],obs[1])
    uL = get_uc(αs[end],obs[end])

    # Get W's
    Ws = [get_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:length(βs)]
    # print_log("$(prod(Ws))")

    # Scale Ws for numerical stability
    Ws = UniformScaling(1.0/maximum(maximum.(Ws))) * Ws
    # print_log("Scaled: $(prod(Ws))")
    # readline()

    # Get Zc
    Zc = get_Zc(u1,uL,Ws)

    # Ex
    ex = get_Ec_X(u1,uL,Ws,Zc)

    # Exx
    exx = get_Ec_XX(u1,uL,Ws,Zc)

    # Return
    return ex,exx

end
"""

    `get_Ec_logpyx(log_u1,log_uN,log_Ws,log_Zc)`

    Computes E[log p(y|x)|ỹ_m] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-50.0,-50.0),L);
    julia> u1 = CpelNano.get_uc(αs[1],obs[1]);
    julia> uL = CpelNano.get_uc(αs[end],obs[end]);
    julia> Ws = [CpelNano.get_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> Zc = CpelNano.get_Zc(u1,uL,Ws);
    julia> CpelNano.get_Ec_logpyx_log(logu1,loguL,logWs,logZc,obs)

    ```
"""
function get_Ec_logpyx(u1::Vector{Float64},uL::Vector{Float64},Ws::Vector{Array{Float64,2}},Zc::Float64,call_m::Vector{MethCallCpgGrp})::Vector{Float64}

    # Compute E[log p(y|x)|ỹ]
    prod_Ws = prod(Ws)
    ex_log_pyxs = zeros(length(call_m))
    @inbounds for l=1:length(ex_log_pyxs)

        # Skip if no data
        call_m[l].obs || continue

        # Diagonal matrix to introduce
        Dl = [call_m[l].log_pyx_u 0.0; 0.0 call_m[l].log_pyx_m]

        # Get l-th term
        if l==1
            ex_log_pyxs[1] = u1' * Dl * prod_Ws * uL
        elseif l==length(ex_log_pyxs)
            ex_log_pyxs[end] = u1' * prod_Ws * Dl * uL
        else
            ex_log_pyxs[l] = u1' * prod(Ws[1:(l-1)]) * Dl * prod(Ws[l:end]) * uL
        end

    end

    # Return E[log p(y|x)|ỹ]
    return ex_log_pyxs ./ Zc

end
##################################################################################################
## Transfer matrix method functions for expectations E[X̄] && E[X̄X̄]
##################################################################################################
"""

    `get_u(α)`

    Generates vector u1 and uL used in transfer matrix method.

    # Examples
    ```julia-repl
    julia> CpelNano.get_u(0.0)
    2-element Array{Float64,1}:
     1.0
     1.0
    ```
"""
function get_u(α::Float64)::Vector{Float64}

    # Compute intermediate
    aux = exp(0.5 * α)

    # Return u for CpG passed
    return [1.0/aux; aux]

end
"""

    `get_W(α1,α2,β)`

    Generates transition matrix W used in transfer matrix method.

    # Examples
    ```julia-repl
    julia> CpelNano.get_W(0.0,0.0,0.0)
    2×2 Array{Float64,2}:
     1.0  1.0
     1.0  1.0
    ```
"""
function get_W(α1::Float64,α2::Float64,β::Float64)::Array{Float64,2}

    # Intermediate quantities
    exp_α1 = exp(0.5*α1)
    exp_α2 = exp(0.5*α2)
    exp_β = exp(β)

    # Define transition matrix
    W = zeros(2,2)
    W[1,1] = exp_β / (exp_α1 * exp_α2)
    W[2,1] = exp_α1 / (exp_α2 * exp_β)
    W[1,2] = exp_α2 / (exp_α1 * exp_β)
    W[2,2] = exp_α1 * exp_α2 * exp_β

    # Return W
    return W

end
"""

    `get_Z(u1,uL,Ws)`

    Computes partition function of p(x̄;ϕ) using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> u1 = CpelNano.get_u(αs[1]); uL = CpelNano.get_u(αs[end]);
    julia> Ws = [CpelNano.get_W(αs[l],αs[l+1],βs[l]) for l=1:length(βs)];
    julia> CpelNano.get_Z(u1,uL,Ws)
    1024.0
    ```
"""
function get_Z(u1::Vector{Float64},uL::Vector{Float64},Ws::Vector{Array{Float64,2}})::Float64

    # Return u_1⊺ W_1 ⋯ W_L uL of p(x̄;ϕ)
    return u1' * prod(Ws) * uL

end
"""

    `get_E_X(u1,uN,Ws,Z)`

    Computes E[X̄;ϕ] using the transfer matrix method. CAUTION: This function is not as numerically
    stable as `get_E_X_log(logu1,loguN,logWs,logZ)`, albeit significantly faster.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> u1 = CpelNano.get_u(αs[1]); uL = CpelNano.get_u(αs[end]);
    julia> Ws = [CpelNano.get_W(αs[l],αs[l+1],βs[l]) for l=1:length(βs)];
    julia> Ws = UniformScaling(1.0/maximum(maximum.(Ws))) * Ws
    julia> Z = CpelNano.get_Z(u1,uL,Ws);
    julia> CpelNano.get_E_X(u1,uL,Ws,Z)
    10-element Array{Float64,1}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
    ```
"""
function get_E_X(u1::Vector{Float64},uL::Vector{Float64},Ws::Vector{Array{Float64,2}},Z::Float64)::Vector{Float64}

    # Compute E[X̄_l;ϕ] for l=1,2,…,L
    prod_Ws = prod(Ws)
    exs = fill(NaN,length(Ws)+1)
    @inbounds for l=1:length(exs)
        if l==1
            exs[1] = u1' * D_EX * prod_Ws * uL
        elseif l==length(exs)
            exs[end] = u1' * prod_Ws * D_EX * uL
        else
            exs[l] = u1' * prod(Ws[1:(l-1)]) * D_EX * prod(Ws[l:end]) * uL
        end
    end

    # Return vector of E[X̄_l;ϕ] for l=1,2,…,L
    return exs ./ Z

end
"""

    `get_E_XX(u1,uL,Ws,Z)`

    Computes E[X̄_{l}X̄_{l+1};ϕ] using the transfer matrix method. CAUTION: This function is not as numerically
    stable as `get_E_XX_log(logu1,loguN,logWs,logZ)`, albeit significantly faster.

    # Examples
    ```julia-repl
    julia> using LinearAlgebra
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> u1 = CpelNano.get_u(αs[1]); uL = CpelNano.get_u(αs[end]);
    julia> Ws = [CpelNano.get_W(αs[l],αs[l+1],βs[l]) for l=1:length(βs)];
    julia> Z = CpelNano.get_Z(u1,uL,Ws);
    julia> Ws = UniformScaling(1.0/maximum(maximum.(Ws))) * Ws
    julia> CpelNano.get_E_XX(u1,uL,Ws,Z)
    9-element Array{Float64,1}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
    ```
"""
function get_E_XX(u1::Vector{Float64},uL::Vector{Float64},Ws::Vector{Array{Float64,2}},Z::Float64)::Vector{Float64}

    # Init output
    exxs = fill(NaN,length(Ws))

    # If only one W
    if length(Ws)==1
        exxs[1] = u1'*(D_EXX.* Ws[1]) * uL
        return exxs ./ Z
    end

    # Compute E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] for l=1,2,…,L-1
    @inbounds for l=1:length(exxs)
        if l==1
            exxs[1] = u1' * (D_EXX .* Ws[1]) * prod(Ws[2:end]) * uL
        elseif l==length(Ws)
            exxs[end] = u1' * prod(Ws[1:(end-1)]) * (Ws[end] .* D_EXX) * uL
        else
            exxs[l] = u1' * prod(Ws[1:(l-1)]) * (Ws[l] .* D_EXX) * prod(Ws[(l+1):end]) * uL
        end
    end

    # Return vector of E[X̄_{l}X̄_{l+1};ϕ] for l=1,2,…,L-1
    return exxs ./ Z

end
"""

    `get_marg_px(u1,uL,Ws,Z)`

    Computes marginal probability P(X̄_l=1;ϕ) using the transfer matrix method. CAUTION: This function is not as
    numerically stable as `get_marg_px_log(logu1,loguN,logWs,logZ)`, albeit significantly faster.

    # Examples
    ```julia-repl
    julia> using LinearAlgebra
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> u1 = CpelNano.get_u(αs[1]);
    julia> uL = CpelNano.get_u(αs[end]);
    julia> Ws = [CpelNano.get_W(αs[l],αs[l+1],βs[l]) for l=1:(L-1)];
    julia> Ws = UniformScaling(1.0/maximum(maximum.(Ws))) * Ws
    julia> Z = CpelNano.get_Z(u1,uL,Ws);
    julia> CpelNano.get_marg_px(u1,uL,Ws,Z)
    10-element Array{Float64,1}:
     0.5
     0.5
     0.5
     0.5
     0.5
     0.5
     0.5
     0.5
     0.5
    ```
"""
function get_marg_px(u1::Vector{Float64},uL::Vector{Float64},Ws::Vector{Array{Float64,2}},Z::Float64)::Vector{Float64}

    # Compute P(X̄_l=1;ϕ) for l=1,2,…,L
    prod_Ws = prod(Ws)
    exs = fill(NaN,length(Ws)+1)
    @inbounds for l=1:length(exs)
        if l==1
            exs[1] = u1' * D_PX * prod_Ws * uL
        elseif l==length(exs)
            exs[end] = u1' * prod_Ws * D_PX * uL
        else
            exs[l] = u1' * prod(Ws[1:(l-1)]) * D_PX * prod(Ws[l:end]) * uL
        end
    end

    # Return vector of P(X̄_l=1;ϕ) for l=1,2,…,L
    return exs ./ Z

end
##################################################################################################
## Transfer matrix base functions
##################################################################################################
"""

    `log_sum_exp(x)`

    Computes log sum using vector x as log(1+sum(exp(x-x_max)))+xmax.

    # Examples
    ```julia-repl
    julia> CpelNano.log_sum_exp([1.0e-2,1.0e-5,1.0e-10])
    0.6981596805576123
    ```
"""
function log_sum_exp(x::Vector{Float64})::Float64

    # Get max and min
    x_max = maximum(x)
    x_min = minimum(x)

    # Return
    return log(1.0+exp(x_min-x_max))+x_max

end
"""

    `log_mat_mult(logW1,logW2)`

    Computes log(W1*W2) in a numeric stable way.

    # Examples
    ```julia-repl
    julia> CpelNano.log_mat_mult(rand(2,2),rand(2,2))
    ```
"""
function log_mat_mult(logW1::Array{Float64,2},logW2::Array{Float64,2})::Array{Float64,2}

    # Init output
    logW1W2 = fill(NaN,(2,2))

    # Compute each entry
    logW1W2[1,1] = log_sum_exp(logW1[1,:]+logW2[:,1])
    logW1W2[1,2] = log_sum_exp(logW1[1,:]+logW2[:,2])
    logW1W2[2,1] = log_sum_exp(logW1[2,:]+logW2[:,1])
    logW1W2[2,2] = log_sum_exp(logW1[2,:]+logW2[:,2])

    # Return log(W1*W2)
    return logW1W2

end
"""

    `log_vec_mat_mult(logu,logW)`

    Computes log(u'*W) in a numeric stable way.

    # Examples
    ```julia-repl
    julia> CpelNano.log_vec_mat_mult(rand(2),rand(2,2))
    ```
"""
function log_vec_mat_mult(logu::Vector{Float64},logW::Array{Float64,2})::Vector{Float64}

    # Init output
    loguW = fill(NaN,2)

    # Compute each entry
    loguW[1] = log_sum_exp(logu+logW[:,1])
    loguW[2] = log_sum_exp(logu+logW[:,2])

    # Return log(u'*W)
    return loguW

end
"""

    `log_vec_vec_mult(logu1,logu2)`

    Computes log(u1'*u2) in a numeric stable way.

    # Examples
    ```julia-repl
    julia> CpelNano.log_vec_vec_mult(rand(2),rand(2))
    ```
"""
function log_vec_vec_mult(logu1::Vector{Float64},logu2::Vector{Float64})::Float64

    # Return log(u1'*u2)
    return log_sum_exp(logu1+logu2)

end
"""

    `mult_log_mats([log(M1),log(M2),...])`

    Computes log(M1*M2*...) in a numeric stable way.

    # Examples
    ```julia-repl
    julia> CpelNano.mult_log_mats(fill(rand(2,2),10))
    ```
"""
function mult_log_mats(logMs::Vector{Array{Float64,2}})::Array{Float64,2}

    # Check if only one matrix
    length(logMs)==1 && return logMs[1]

    # Perform first multiplication
    mult_logMs = log_mat_mult(logMs[1],logMs[2])
    length(logMs)==2 && return mult_logMs

    # Do the rest
    @inbounds for i=3:length(logMs)
        mult_logMs = log_mat_mult(mult_logMs,logMs[i])
    end

    # Return log of the multiplication
    return mult_logMs

end
##################################################################################################
## Transfer matrix method functions for conditional expectations E[⋅|y] using log trick
##################################################################################################
"""

    `get_γc(x,α,call)`

    Computes function γ in transfer matrix method for p(x|y).

    # Examples
    ```julia-repl
    julia> call = CpelNano.MethCallCpgGrp(-50.0,-45.0);
    julia> CpelNano.get_γc(false,0.0,call)
    -25.0
    julia> CpelNano.get_γc(true,0.0,call)
    -22.5
    ```
"""
function get_γc(x::Bool,α::Float64,call::MethCallCpgGrp)::Float64

    # Compute γ_{m,l} = 0.5⋅(α⋅x + 1(obs)⋅log p(y|x))
    γ = x ? α : - α
    if call.obs
        γ += x ? call.log_pyx_m : call.log_pyx_u
    end

    # Return
    return 0.5 * γ

end
"""

    `get_log_uc(α,call)`

    Generates vector log(u1) and log(uN) used in transfer matrix method for p(x|y).

    # Examples
    ```julia-repl
    julia> call = CpelNano.MethCallCpgGrp(-50.0,-45.0);
    julia> CpelNano.get_log_uc(0.0,call)
     -25.0
     -22.5
    ```
"""
function get_log_uc(α::Float64,call::MethCallCpgGrp)::Vector{Float64}

    # Return u for CpG group passed
    return [get_γc(false,α,call); get_γc(true,α,call)]

end
"""

    `get_log_Wc(α1,α2,β,call1,call2)`

    Generates transition matrix log(W) used in transfer matrix method for p(x|y).

    # Examples
    ```julia-repl
    julia> call = CpelNano.MethCallCpgGrp(-50.0,-45.0);
    julia> Wc = CpelNano.get_log_Wc(0.0,0.0,0.0,call,call)
    2×2 Array{Float64,2}:
     -50.0  -47.5
     -47.5  -45.0
    ```
"""
function get_log_Wc(α1::Float64,α2::Float64,β::Float64,call1::MethCallCpgGrp,call2::MethCallCpgGrp)::Array{Float64,2}

    # Define transition matrix
    log_W = zeros(2,2)
    log_W[1,1] = get_γc(false,α1,call1) + get_γc(false,α2,call2) + β
    log_W[2,1] = get_γc(true,α1,call1) + get_γc(false,α2,call2) - β
    log_W[1,2] = get_γc(false,α1,call1) + get_γc(true,α2,call2) - β
    log_W[2,2] = get_γc(true,α1,call1) + get_γc(true,α2,call2) + β

    # Return log_W
    return log_W

end
"""

    `get_log_Zc(log_u1,log_uN,log_Ws)`

    Computes log of partition function of P(X|Y;ϕ) using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-100.0,-100.0),L);
    julia> logu1 = CpelNano.get_log_uc(αs[1],obs[1]);
    julia> loguL = CpelNano.get_log_uc(αs[end],obs[end]);
    julia> logWs = [CpelNano.get_log_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> logZc = CpelNano.get_log_Zc(logu1,loguL,logWs)
    -993.06
    ```
"""
function get_log_Zc(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}})::Float64

    # Return log u_1⊺ W_1 ⋯ W_N u_N for m-th observation
    return log_vec_vec_mult(log_vec_mat_mult(logu1,mult_log_mats(logWs)),loguL)

end
"""

    `get_Ec_X_log(log_u1,log_uN,log_Ws,log_Zc)`

    Computes E[X̄_l|ỹ_m] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-45.0,-45.0),L);
    julia> logu1 = CpelNano.get_log_uc(αs[1],obs[1]);
    julia> loguL = CpelNano.get_log_uc(αs[end],obs[end]);
    julia> logWs = [CpelNano.get_log_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> logZc = CpelNano.get_log_Zc(logu1,loguL,logWs);
    julia> CpelNano.get_Ec_X_log(logu1,loguL,logWs,logZc)
    10-element Array{Float64,1}:
     3.552713678800501e-15
     3.552713678800501e-15
     3.552713678800501e-15
     3.552713678800501e-15
     3.552713678800501e-15
     3.552713678800501e-15
     3.552713678800501e-15
     3.552713678800501e-15
     3.552713678800501e-15
     3.552713678800501e-15
    ```
"""
function get_Ec_X_log(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZc::Float64)::Vector{Float64}

    # Compute E[X̄_{m,l}|ỹ_m] for l=1,2,…,L
    log_exs = Vector{Float64}()
    @inbounds for l=1:(length(logWs)+1)
        if l==1
            log_DWs = mult_log_mats(vcat([logD1],logWs))
            log_u1DWs = log_vec_mat_mult(logu1,log_DWs)
            push!(log_exs,log_vec_vec_mult(log_u1DWs,loguL))
        elseif l==(length(logWs)+1)
            log_WsD = mult_log_mats(vcat(logWs,[logD1]))
            log_u1WsD = log_vec_mat_mult(logu1,log_WsD)
            push!(log_exs,log_vec_vec_mult(log_u1WsD,loguL))
        else
            log_WsDWs = mult_log_mats(vcat(logWs[1:(l-1)],[logD1],logWs[l:end]))
            log_u1WsDWs = log_vec_mat_mult(logu1,log_WsDWs)
            push!(log_exs,log_vec_vec_mult(log_u1WsDWs,loguL))
        end
    end

    # Return E[X̄_m|ỹ_m]
    return exp.(log_exs .- logZc) .- 2.0

end
"""

    `get_Ec_logpyx_log(log_u1,log_uN,log_Ws,log_Zc)`

    Computes E[log p(y|x)|ỹ_m] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(-1.0,L); βs = fill(1.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-105.0,-95.0),L);
    julia> logu1 = CpelNano.get_log_uc(αs[1],obs[1]);
    julia> loguL = CpelNano.get_log_uc(αs[end],obs[end]);
    julia> logWs = [CpelNano.get_log_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> logZc = CpelNano.get_log_Zc(logu1,loguL,logWs);
    julia> CpelNano.get_Ec_logpyx_log(logu1,loguL,logWs,logZc,obs)
    10-element Array{Float64,1}:
     -95.00045412823272
     -95.00006161154741
     -95.00006148227914
     -95.00006148224674
     -95.00006148224674
     -95.00006148223594
     -95.00006148223594
     -95.00006148226834
     -95.00006161153661
     -95.00045412823272
    ```
"""
function get_Ec_logpyx_log(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZc::Float64,call_m::Vector{MethCallCpgGrp})::Vector{Float64}

    # Compute E[log p(y|x)|ỹ]
    log_pyxs = zeros(length(call_m))
    @inbounds for l=1:(length(logWs)+1)
        call_m[l].obs || continue
        min_logpyx = min(call_m[l].log_pyx_u,call_m[l].log_pyx_m) - 1.0 - eps()
        # println(min_logpyx)
        # logD = log.(max.([call_m[l].log_pyx_u+250.0 0.0; 0.0 call_m[l].log_pyx_m+250.0],1.0))
        logD = log.([call_m[l].log_pyx_u-min_logpyx 0.0; 0.0 call_m[l].log_pyx_m-min_logpyx])
        if l==1
            log_DWs = mult_log_mats(vcat([logD],logWs))
            log_u1DWs = log_vec_mat_mult(logu1,log_DWs)
            log_pyxs[l] = log_vec_vec_mult(log_u1DWs,loguL)
        elseif l==(length(logWs)+1)
            log_WsD = mult_log_mats(vcat(logWs,[logD]))
            log_u1WsD = log_vec_mat_mult(logu1,log_WsD)
            log_pyxs[l] = log_vec_vec_mult(log_u1WsD,loguL)
        else
            log_WsDWs = mult_log_mats(vcat(logWs[1:(l-1)],[logD],logWs[l:end]))
            log_u1WsDWs = log_vec_mat_mult(logu1,log_WsDWs)
            log_pyxs[l] = log_vec_vec_mult(log_u1WsDWs,loguL)
        end
        log_pyxs[l] = exp(log_pyxs[l] - logZc) - min_logpyx
    end

    # Return E[log p(y|x)|ỹ]
    return log_pyxs

end
"""

    `get_logpy_lgtrck(log_u1,log_uN,log_Ws,log_Z,obs)`

    Computes log p(y1,...,yL) using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(-1.0,L); βs = fill(1.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-105.0,-95.0),L);
    julia> logu1 = CpelNano.get_log_u(αs[1]);
    julia> loguL = CpelNano.get_log_u(αs[end]);
    julia> logWs = [CpelNano.get_log_W(αs[l],αs[l+1],βs[l]) for l=1:(L-1)];
    julia> logZc = CpelNano.get_log_Z(logu1,loguL,logWs);
    julia> CpelNano.get_logpy_lgtrck(logu1,loguL,logWs,logZ,obs)

    ```
"""
function get_logpy_lgtrck(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZ::Float64,call_m::Vector{MethCallCpgGrp})::Float64

    # Compute log p(y)
    log_mats = fill(0.0,(2,2))
    @inbounds for l=1:length(logWs)
        if call_m[l].obs
            logD = [call_m[l].log_pyx_u -Inf; -Inf call_m[l].log_pyx_m]
            log_DWs = mult_log_mats(vcat([logD],[logWs[l]]))
            log_mats = l==1 ? log_DWs : mult_log_mats(vcat([log_mats],[log_DWs]))
        else
            log_mats = l==1 ? logWs[l] : mult_log_mats(vcat([log_mats],[logWs[l]]))
        end
    end

    # Take care of last group
    if call_m[end].obs
        logD = [call_m[end].log_pyx_u -Inf; -Inf call_m[end].log_pyx_m]
        log_mats = mult_log_mats(vcat([log_mats],[logD]))
    end

    # Compute
    log_mats = log_vec_mat_mult(logu1,log_mats)
    log_mats = log_vec_vec_mult(log_mats,loguL)

    # Return log p(y)
    return log_mats-logZ

end
"""

    `get_Ec_XX_log(log_u1,log_uN,log_Ws,log_Zc)`

    Computes vector of E[X̄_{m,l}X̄_{m,l+1}|ỹ_m] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-100.0,-100.0),L);
    julia> logu1 = CpelNano.get_log_uc(αs[1],obs[1]);
    julia> loguL = CpelNano.get_log_uc(αs[end],obs[end]);
    julia> logWs = [CpelNano.get_log_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> logZc = CpelNano.get_log_Zc(logu1,loguL,logWs);
    julia> CpelNano.get_Ec_XX_log(logu1,loguL,logWs,logZc)
    9-element Array{Float64,1}:
     -1.099120794378905e-13
     -1.099120794378905e-13
     -1.099120794378905e-13
     -1.099120794378905e-13
     -1.099120794378905e-13
     1.1723955140041653e-13
     1.1723955140041653e-13
     1.1723955140041653e-13
     1.1723955140041653e-13
    ```
"""
function get_Ec_XX_log(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZc::Float64)::Vector{Float64}

    # Compute E[X̄_{m,l}X̄_{m,l+1}|ỹ_m] for l=1,2,…,L-1
    log_exxs = Vector{Float64}()
    @inbounds for l=1:length(logWs)
        if l==1
            aux = logWs[1]+logD2
            log_DWs = mult_log_mats(vcat([aux],logWs[2:end]))
            log_u1DWs = log_vec_mat_mult(logu1,log_DWs)
            push!(log_exxs,log_vec_vec_mult(log_u1DWs,loguL))
        elseif l==length(logWs)
            aux = logWs[end]+logD2
            log_WsD = mult_log_mats(vcat(logWs[1:(end-1)],[aux]))
            log_u1WsD = log_vec_mat_mult(logu1,log_WsD)
            push!(log_exxs,log_vec_vec_mult(log_u1WsD,loguL))
        else
            aux = logWs[l]+logD2
            log_WsDWs = mult_log_mats(vcat(logWs[1:(l-1)],[aux],logWs[(l+1):end]))
            log_u1WsDWs = log_vec_mat_mult(logu1,log_WsDWs)
            push!(log_exxs,log_vec_vec_mult(log_u1WsDWs,loguL))
        end

    end

    # Return vector of E[X̄_{m,l}X̄_{m,l+1}|ỹ_m]
    return exp.(log_exxs .- logZc) .- 2.0

end
"""

    `get_cond_exs_log(αs,βs,obs)`

    Computes conditional expectations using the transfer matrix method. This function uses the
    log-sum trick to cope with the numerical instability.

    # Examples
    ```julia-repl
    julia> L = 2; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-50.0,-45.0),L);
    julia> ex,exx = CpelNano.get_cond_exs_log(αs,βs,obs)
    ([0.9866142981513866, 0.9866142981513866], [0.9734077733168003])
    julia> obs = fill(CpelNano.MethCallCpgGrp(-200.0,-200.0),L);
    julia> ex,exx = CpelNano.get_cond_exs_log(αs,βs,obs)
    ([3.552713678800501e-15, 3.552713678800501e-15], [3.552713678800501e-15])
    ```
"""
function get_cond_exs_log(αs::Vector{Float64},βs::Vector{Float64},obs::Vector{MethCallCpgGrp})::NTuple{2,Vector{Float64}}

    # Get log u's
    logu1 = get_log_uc(αs[1],obs[1])
    loguL = get_log_uc(αs[end],obs[end])

    # Get log W's
    logWs = [get_log_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:length(βs)]

    # Get Zc
    logZc = get_log_Zc(logu1,loguL,logWs)

    # Ex
    ex = get_Ec_X_log(logu1,loguL,logWs,logZc)

    # Exx
    exx = get_Ec_XX_log(logu1,loguL,logWs,logZc)

    # Return E[X̄] and E[X̄X̄]
    return ex,exx

end
##################################################################################################
## Transfer matrix method functions for WGBS inference using log sum trick (not used)
##################################################################################################
"""

    `get_Ec_X_wgbs_log(u1,uL,Ws,Zc,obs)`

    Computes E[X̄|y_m;ϕ] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 2; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-2.0,0.0),L);
    julia> logu1 = CpelNano.get_log_uc(αs[1],obs[1]); loguL = CpelNano.get_log_uc(αs[end],obs[end]);
    julia> logWs = [CpelNano.get_log_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> logZc = CpelNano.get_log_Zc(logu1,loguL,logWs);
    julia> CpelNano.get_Ec_X_wgbs_log(logu1,loguL,logWs,logZc,obs)
    2-element Array{Float64,1}:
     1.0
     1.0
    ```
"""
function get_Ec_X_wgbs_log(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZc::Float64,obs::Vector{MethCallCpgGrp})::Vector{Float64}

    # Compute E[X̄_l|y_m] for l=1,2,…,L
    log_exs = fill(NaN,length(obs))
    @inbounds for l=1:length(obs)

        # Check if observed
        if obs[l].obs

            # If observed just take observed value & continue
            log_exs[l] = obs[l].log_pyx_m > obs[l].log_pyx_u ? log(3.0) : log(1.0)

        else

            # Compute E[X̄_{m,l}|ỹ_m] for l=1,2,…,L
            if l==1
                log_DWs = mult_log_mats(vcat([logD1],logWs))
                log_u1DWs = log_vec_mat_mult(logu1,log_DWs)
                log_exs[l] = log_vec_vec_mult(log_u1DWs,loguL) - logZc
            elseif l==(length(logWs)+1)
                log_WsD = mult_log_mats(vcat(logWs,[logD1]))
                log_u1WsD = log_vec_mat_mult(logu1,log_WsD)
                log_exs[l] = log_vec_vec_mult(log_u1WsD,loguL) - logZc
            else
                log_WsDWs = mult_log_mats(vcat(logWs[1:(l-1)],[logD1],logWs[l:end]))
                log_u1WsDWs = log_vec_mat_mult(logu1,log_WsDWs)
                log_exs[l] = log_vec_vec_mult(log_u1WsDWs,loguL) - logZc
            end

        end

    end

    # Return vector of Ec[X̄_l|y_m] for m-th observation
    return exp.(log_exs) .- 2.0
    # return exs

end
"""

    `get_Ec_XX_wgbs_log(u1,uL,Ws,obs)`

    Computes E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 3; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-5.0,0.0),L);
    julia> logu1 = CpelNano.get_log_uc(αs[1],obs[1]); loguL = CpelNano.get_log_uc(αs[end],obs[end]);
    julia> logWs = [CpelNano.get_log_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:(L-1)];
    julia> logZc = CpelNano.get_log_Zc(logu1,loguL,logWs);
    julia> CpelNano.get_Ec_XX_wgbs_log(logu1,loguL,logWs,logZc,obs)
    2-element Array{Float64,1}:
     1.0
     1.0
    ```
"""
function get_Ec_XX_wgbs_log(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZc::Float64,obs::Vector{MethCallCpgGrp})::Vector{Float64}

    # Init output
    log_exxs = fill(NaN,length(logWs))

    # Compute E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] for l=1,2,…,L-1
    @inbounds for l=1:length(logWs)

        if obs[l].obs && obs[l+1].obs

            # If both are observed just take observed value & continue
            xl = obs[l].log_pyx_m > obs[l].log_pyx_u ? 1.0 : -1.0
            xlp1 = obs[l+1].log_pyx_m > obs[l+1].log_pyx_u ? 1.0 : -1.0
            log_exxs[l] = xl==xlp1 ? log(3.0) : log(1.0)

        else

            # If not observed compute E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ]
            if l==1
                aux = logWs[1]+logD2
                log_DWs = mult_log_mats(vcat([aux],logWs[2:end]))
                log_u1DWs = log_vec_mat_mult(logu1,log_DWs)
                log_exxs[l] = log_vec_vec_mult(log_u1DWs,loguL) - logZc
            elseif l==length(logWs)
                aux = logWs[end]+logD2
                log_WsD = mult_log_mats(vcat(logWs[1:(end-1)],[aux]))
                log_u1WsD = log_vec_mat_mult(logu1,log_WsD)
                log_exxs[l] = log_vec_vec_mult(log_u1WsD,loguL) - logZc
            else
                aux = logWs[l]+logD2
                log_WsDWs = mult_log_mats(vcat(logWs[1:(l-1)],[aux],logWs[(l+1):end]))
                log_u1WsDWs = log_vec_mat_mult(logu1,log_WsDWs)
                log_exxs[l] = log_vec_vec_mult(log_u1WsDWs,loguL) - logZc
            end

        end

    end

    # Return vector of E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] for m-th observation
    return exp.(log_exxs) .- 2.0

end
"""

    `get_cond_exs_wgbs_log(αs,βs,obs)`

    Computes E[X̄_{m,l}|Y_m=y_m;ϕ] & E[X̄_{m,l}X̄_{m,l+1}|Y_m=y_m;ϕ] using the transfer matrix method. This function
    does not use the log-sum trick, however. Thus, it is not as numerically stable as `get_cond_exs_log()`.

    # Examples
    ```julia-repl
    julia> L = 5; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> obs = fill(CpelNano.MethCallCpgGrp(-2.0,0.0),L);
    julia> ex,exx = CpelNano.get_cond_exs_wgbs_log(αs,βs,obs);
    julia> obs = fill(CpelNano.MethCallCpgGrp(0.0,-2.0),L);
    julia> ex,exx = CpelNano.get_cond_exs_wgbs_log(αs,βs,obs);
    ```
"""
function get_cond_exs_wgbs_log(αs::Vector{Float64},βs::Vector{Float64},obs::Vector{MethCallCpgGrp})::NTuple{2,Vector{Float64}}

    ## Get matrices

    # Get u's
    logu1 = get_log_uc(αs[1],obs[1])
    loguL = get_log_uc(αs[end],obs[end])

    # Get W's
    logWs = [get_log_Wc(αs[l],αs[l+1],βs[l],obs[l],obs[l+1]) for l=1:length(βs)]

    # Get Zc
    logZc = get_log_Zc(logu1,loguL,logWs)

    ## Compute expectations

    # Ex
    ex = get_Ec_X_wgbs_log(logu1,loguL,logWs,logZc,obs)

    # Exx
    exx = get_Ec_XX_wgbs_log(logu1,loguL,logWs,logZc,obs)

    # Return
    return ex,exx

end
##################################################################################################
## Transfer matrix method functions for expectations E[X̄] && E[X̄X̄] using log-sum trick
##################################################################################################
"""

    `get_log_u(α)`

    Generates vector log(u) used in transfer matrix method.

    # Examples
    ```julia-repl
    julia> CpelNano.get_log_u(0.0)
    2-element Array{Float64,1}:
     0.0
     0.0
    ```
"""
function get_log_u(α::Float64)::Vector{Float64}

    # Compute intermediate quantity
    aux = 0.5 * α

    # Return u for CpG passed
    return [-aux;aux]

end
"""

    `get_log_W(α1,α2,β)`

    Generates transition matrix log(W) used in transfer matrix method.

    # Examples
    ```julia-repl
    julia> CpelNano.get_log_W(0.0,0.0,0.0)
    2×2 Array{Float64,2}:
     0.0  0.0
     0.0  0.0
    ```
"""
function get_log_W(α1::Float64,α2::Float64,β::Float64)::Array{Float64,2}

    # Define transition matrix
    logW = zeros(2,2)
    logW[1,1] = -0.5*(α1+α2) + β
    logW[2,1] = 0.5*(α1-α2) - β
    logW[1,2] = 0.5*(α2-α1) - β
    logW[2,2] = 0.5*(α1+α2) + β

    # Return logW
    return logW

end
"""

    `get_log_Z(u1,uN,Ws)`

    Computes log of partition function of q(x̄;ϕ) using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> logu1 = CpelNano.get_log_u(αs[1]);
    julia> loguL = CpelNano.get_log_u(αs[end]);
    julia> logWs = [CpelNano.get_log_W(αs[l],αs[l+1],βs[l]) for l=1:length(βs)];
    julia> logZ = CpelNano.get_log_Z(logu1,loguL,logWs)
    6.931471805599453
    ```
"""
function get_log_Z(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}})::Float64

    # Compute intermediate quantities
    aux = mult_log_mats(logWs)
    aux = log_vec_mat_mult(logu1,aux)

    # Return log u_1⊺ W_1 ⋯ W_L u_L
    return log_vec_vec_mult(aux,loguL)

end
"""
    `get_∇logZ(REG_ST,ϕhat)`

    Numerically computes the gradient of the log partition function of a model with estimated parameter
    vector ϕhat=[a,b,c]. This is equivalent to computing the expected value of sufficient statistis (SS).

    # Examples
    ```julia-repl
    julia> x=CpelNano.RegStruct(); x.L=10; x.Nl=fill(5.0,x.L); x.ρl=fill(0.1,x.L); x.dl=fill(10.0,x.L);
    julia> CpelNano.get_∇logZ(x,[0.0,0.0,0.0])
    3-element Array{Float64,1}:
     0.0
     0.0
     0.0
    ```
"""
function get_∇logZ(rg::RegStruct,ϕhat::Vector{Float64})::Vector{Float64}

    # Define function
    function f(ϕ::Vector{Float64})

        # Get αs & βs
        αs,βs = get_αβ_from_ϕ(ϕ,rg)

        # Get matrices
        logu1 = get_log_u(αs[1])
        loguL = get_log_u(αs[end])
        logWs = [get_log_W(αs[l],αs[l+1],βs[l]) for l=1:length(βs)]

        # Return log ζ(ϕ)
        return get_log_Z(logu1,loguL,logWs)

    end

    # Return ∇logζ(ϕ)
    return Calculus.gradient(f,ϕhat)

end
"""

    `get_E_X_log(log_u1,log_uL,log_Ws,log_Z)`

    Computes E[X̄_l;ϕ] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> logu1 = CpelNano.get_log_u(αs[1]);
    julia> loguL = CpelNano.get_log_u(αs[end]);
    julia> logWs = [CpelNano.get_log_W(αs[l],αs[l+1],βs[l]) for l=1:(L-1)];
    julia> logZ = CpelNano.get_log_Z(logu1,loguL,logWs);
    julia> CpelNano.get_E_X_log(logu1,loguL,logWs,logZ)
    10-element Array{Float64,1}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
    ```
"""
function get_E_X_log(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZ::Float64)::Vector{Float64}

    # Compute E[X̄_l;ϕ] for l=1,2,…,L
    log_exs = Vector{Float64}()
    @inbounds for l=1:(length(logWs)+1)
        if l==1
            log_DWs = mult_log_mats(vcat([logD1],logWs))
            log_u1DWs = log_vec_mat_mult(logu1,log_DWs)
            push!(log_exs,log_vec_vec_mult(log_u1DWs,loguL))
        elseif l==(length(logWs)+1)
            log_WsD = mult_log_mats(vcat(logWs,[logD1]))
            log_u1WsD = log_vec_mat_mult(logu1,log_WsD)
            push!(log_exs,log_vec_vec_mult(log_u1WsD,loguL))
        else
            log_WsDWs = mult_log_mats(vcat(logWs[1:(l-1)],[logD1],logWs[l:end]))
            log_u1WsDWs = log_vec_mat_mult(logu1,log_WsDWs)
            push!(log_exs,log_vec_vec_mult(log_u1WsDWs,loguL))
        end
    end

    # Return vector of E[X̄_l;ϕ] for l=1,2,…,L
    return exp.(log_exs .- logZ) .- 2.0

end
"""

    `get_E_XX_log(log_u1,log_uL,log_Ws,log_Z)`

    Computes E[X̄_{l}X̄_{l+1};ϕ] using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> logu1 = CpelNano.get_log_u(αs[1]);
    julia> loguL = CpelNano.get_log_u(αs[end]);
    julia> logWs = [CpelNano.get_log_W(αs[l],αs[l+1],βs[l]) for l=1:(L-1)];
    julia> logZ = CpelNano.get_log_Z(logu1,loguL,logWs);
    julia> CpelNano.get_E_XX_log(logu1,loguL,logWs,logZ)
    9-element Array{Float64,1}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
    ```
"""
function get_E_XX_log(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZ::Float64)::Vector{Float64}

    # Compute E[X̄_{l}X̄_{l+1};ϕ] for l=1,2,…,L-1
    log_exxs = Vector{Float64}()
    @inbounds for l=1:length(logWs)
        if l==1
            aux = logWs[1]+logD2
            log_DWs = mult_log_mats(vcat([aux],logWs[2:end]))
            log_u1DWs = log_vec_mat_mult(logu1,log_DWs)
            push!(log_exxs,log_vec_vec_mult(log_u1DWs,loguL))
        elseif l==length(logWs)
            aux = logWs[end]+logD2
            log_WsD = mult_log_mats(vcat(logWs[1:(end-1)],[aux]))
            log_u1WsD = log_vec_mat_mult(logu1,log_WsD)
            push!(log_exxs,log_vec_vec_mult(log_u1WsD,loguL))
        else
            aux = logWs[l]+logD2
            log_WsDWs = mult_log_mats(vcat(logWs[1:(l-1)],[aux],logWs[(l+1):end]))
            log_u1WsDWs = log_vec_mat_mult(logu1,log_WsDWs)
            push!(log_exxs,log_vec_vec_mult(log_u1WsDWs,loguL))
        end

    end

    # Return vector of E[X̄_{l}X̄_{l+1};ϕ] for l=1,2,…,L-1
    return exp.(log_exxs .- logZ) .- 2.0

end
"""

    `get_rs_exps!(REG_STRUCT)`

    Stores E[X] and E[XX] in estimation region structure.

    # Examples
    ```julia-repl
    julia> get_rs_exps!(rs)
    ```
"""
function get_rs_exps!(rs::RegStruct)::Nothing

    # Get E[X]
    rs.exps.ex = get_E_X_log(rs.tm.log_u1,rs.tm.log_uN,rs.tm.log_Ws,rs.logZ)

    # Get E[XX]
    rs.exps.exx = get_E_XX_log(rs.tm.log_u1,rs.tm.log_uN,rs.tm.log_Ws,rs.logZ)

    # Return nothing
    return nothing

end
"""

    `get_marg_px_log(log_u1,log_uL,log_Ws,log_Z)`

    Computes marginal probability p(x_l) using the transfer matrix method.

    # Examples
    ```julia-repl
    julia> L = 10; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> logu1 = CpelNano.get_log_u(αs[1]);
    julia> loguL = CpelNano.get_log_u(αs[end]);
    julia> logWs = [CpelNano.get_log_W(αs[l],αs[l+1],βs[l]) for l=1:(L-1)];
    julia> logZ = CpelNano.get_log_Z(logu1,loguL,logWs);
    julia> CpelNano.get_marg_px_log(logu1,loguL,logWs,logZ)
    10-element Array{Float64,1}:
     0.5
     0.5
     0.5
     0.5
     0.5
     0.5
     0.5
     0.5
     0.5
    ```
"""
function get_marg_px_log(logu1::Vector{Float64},loguL::Vector{Float64},logWs::Vector{Array{Float64,2}},logZ::Float64)::Vector{Float64}

    # Compute p(xl) for l=1,2,…,L
    log_exs = Vector{Float64}()
    @inbounds for l=1:(length(logWs)+1)
        if l==1
            log_DWs = mult_log_mats(vcat([logD3],logWs))
            log_u1DWs = log_vec_mat_mult(logu1,log_DWs)
            push!(log_exs,log_vec_vec_mult(log_u1DWs,loguL))
        elseif l==(length(logWs)+1)
            log_WsD = mult_log_mats(vcat(logWs,[logD3]))
            log_u1WsD = log_vec_mat_mult(logu1,log_WsD)
            push!(log_exs,log_vec_vec_mult(log_u1WsD,loguL))
        else
            log_WsDWs = mult_log_mats(vcat(logWs[1:(l-1)],[logD3],logWs[l:end]))
            log_u1WsDWs = log_vec_mat_mult(logu1,log_WsDWs)
            push!(log_exs,log_vec_vec_mult(log_u1WsDWs,loguL))
        end
    end

    # Return vector of p(xl) for l=1,2,…,L
    return exp.(log_exs .- logZ) .- 2.0

end
"""

    `get_rs_px!(REG_DATA)`

    Computes marginal probability p(x_l) from E[X].

    # Examples
    ```julia-repl
    julia>
    ```
"""
function get_rs_px!(rs::RegStruct)::Nothing

    # Compute p(xl) for l=1,2,…,L
    rs.px = 0.5 .* (rs.exps.ex .+ 1.0)

    # Return nothing
    return nothing

end


text:
function formdataseries(;
  refreshcomp::Bool = true,
  refreshtretcrsp::Bool = true,
  refreshtvolmcapcrsp::Bool = true,
  datapath::String=DATA_PATH,
  levname::String = LEV_NAME,
  injlsstream::Function = IN_JLS_STREAM,
  outjlsstream::Function = OUT_JLS_STREAM,
  outpath::String = OUT_PATH,
  refreshdataseries::Bool = true,
  trialrunonly::Bool = false,
  archive::Bool = true,
  retaineddataseriescolumns = RETAINED_COLS_DATA_SERIES)::DataFrame

  local lev::DataFrame
  local fyear::Int

  local dailyunlesstrial::Symbol = trialrunonly ? :monthly : :daily
  refreshdataseries = (refreshdataseries || refreshcomp)

  levname = "$levname-$dailyunlesstrial"

  if refreshdataseries
    lev = prepcomp(refreshcomp=refreshcomp)

    #supplements the market value and constructs the trailing daily vol
    lev = maketvolmcapseries!(
      lev, refreshtvolmcapcrsp=refreshtvolmcapcrsp, crsptype=dailyunlesstrial)

    lev = maketretseries!(lev, refreshtretcrsp = refreshtretcrsp, crsptype=:monthly)

    winsorizekeydata!(lev)
    constructportfolios!(lev)

    if TEST_OUTPUT
      fyear = 2015
      CSV.write("output\\$(levname)_$fyear.csv", lev[(y->y==fyear).(lev.fyear), :])
    end

    if TEST_OUTPUT
      fyear = 2014
      CSV.write("output\\$(levname)_$fyear.csv", lev[(y->y==fyear).(lev.fyear), :])
    end

    select!(lev, RETAINED_COLS_DATA_SERIES)

    outjlsstream("$datapath\\$levname.jls.lz4") do s
      serialize(s, lev)
    end

    archive && archivefile("$levname.jls.lz4") #never hurts to have a backup
  else

    try
      injlsstream("$datapath\\$levname.jls.lz4") do s
        lev = deserialize(s) #NOTE: DEBUG ONLY
      end
    catch err
      @warn "Could not load $levname. Error: $err"
      @info "Attempting to restore back-up and load from archive..."
      unarchivefile("$levname.jls.lz4")
      injlsstream("$datapath\\$levname.jls.lz4") do s
        lev = deserialize(s) #NOTE: DEBUG ONLY
      end
      @info "Successfully restored from archive."
    end
  end

  return lev
end

function winsorizekeydata!(lev::DataFrame;
    Fvals2winsorize=CQUART_FVALS, winsorizeprop=WINSORIZE_PROP, groupfield=:fyear)

  #dups will lead to race errors
  @assert length(Fvals2winsorize) == length(unique(Fvals2winsorize))

  @mpar for s ∈ Fvals2winsorize
    local nbefore::Int = length(unique(lev[!, s]))
    winsorize!(lev[!, s], prop = winsorizeprop)
    (length(unique(lev[!, s])) == nbefore) && (@warn "winsorization didn't do anything.
      something is probably wrong.")
  end

  return nothing
end

# a couple of helper functions
@inline function rowbeforedatesorted(df::AbstractDataFrame, dt::Date)::Union{Nothing, DataFrameRow}
  local r::Union{DataFrameRow,Nothing}

  i::Int = searchsortedfirst(df.date, dt)

  #check if the index is valid for our purposes
  r = (
    (i>size(df,1)) || #this means the date or a larger date wasn't found
    (i==1)
    ? nothing : df[i-1,:])

  return r
end

function maketretseries!(lev::DataFrame;
  refreshtretcrsp::Bool = true,
  crsptype::Symbol = :daily,
  months2stale::Int = MONTHS_2_STALE_SHORT)

  crspcolumns::Vector{Symbol} = [:date, :permno, :ret, :lret,
    :ret12m, :lret12m, :net12m, :lnet12m]

  periods2stale = Month(months2stale)

  validpermnos::Vector{Int} = unique(lev.lpermno)

  crsp::DataFrame = prepcrsp(
    crspcolumns = crspcolumns,
    crsptype = crsptype,
    refreshcrsp=refreshtretcrsp,
    validpermnos = validpermnos,
    prefix = "tret") do df::DataFrame
    computenetreturns!(df)
    computetrailingreturns!(df, 12, Flret = :lret, Fnewret = :ret12m,
      months2stale = MONTHS_2_STALE_SHORT + 12)
    computetrailingreturns!(df, 12, Flret = :lnet, Fnewret = :net12m,
      months2stale = MONTHS_2_STALE_SHORT + 12)
    computetrailingreturns!(df, 24, Flret = :lret, Fnewret = :ret24m,
      months2stale = MONTHS_2_STALE_SHORT + 24)
  end

  lev = mergecrsp(crsp, lev, refreshmerge=true, cachedata = false, crsptype=crsptype, prefix="tret")
  obliterate!(crsp)

  #supplement the market value
  sort!(lev, [:permno, :date])
  lev.last = falses(size(lev,1))

  #now collapse the return data
  levgrp::GroupedDataFrame = groupby(lev, :permno)
  @mpar for i ∈ 1:length(levgrp)
    slev::SubDataFrame = levgrp[i]

    for sslev ∈ groupby(slev, :compid)
      sslev.last[end] = true
    end
  end

  #now that we have patched the market value, drop all other rows
  #can either free the memory or do this in place
  lev = filter!(r->r.last, lev)

  #now sort and lag
  sort!(lev, [:permno, :fyear])
  lagwithin!(lev, [:ret12m, :lret12m, :net12m, :lnet12m, :ret24m],
    :permno, :fyear, sorted=true)

  lagwithin!(lev, [:Lret12m, :Llret12m, :Lnet12m, :Llnet12m, :Lret24m],
    :permno, :fyear, sorted=true)

  return lev
end

function maketvolmcapseries!(lev::DataFrame;
  refreshtvolmcapcrsp::Bool = true,
  crsptype::Symbol = :daily,
  months2stale::Int = MONTHS_2_STALE_SHORT)

  crspcolumns::Vector{Symbol} = [:date, :permno, :mktcap, :vol252d]
  periods2stale = Month(months2stale)

  validpermnos::Vector{Int} = unique(lev.lpermno)

  crsp::DataFrame = prepcrsp(
    crspcolumns = crspcolumns,
    crsptype = crsptype,
    refreshcrsp=refreshtvolmcapcrsp,
    validpermnos = validpermnos,
    prefix = "tvolmcap") do df::DataFrame
    computecrspmktcap!(df)
    computenetreturns!(df)
    computetrailingvol!(df, crsptype==:daily ? 252 : 12,
      Fret = :net, Fvol = :vol252d, months2stale = MONTHS_2_STALE_SHORT + 12)
    computetrailingvol!(df, crsptype==:daily ? 504 : 24, Fret = :net, Fvol = :vol504d,
    months2stale = MONTHS_2_STALE_SHORT + 24)
  end

  lev = mergecrsp(crsp, lev,
    refreshmerge=true,
    cachedata = false,
    crsptype=crsptype,
    prefix="tvolmcap")

  obliterate!(crsp)

  #supplement the market value
  sort!(lev, [:permno, :date])
  lev.last = falses(size(lev,1))

  #now supplement the market data
  levgrp::GroupedDataFrame = groupby(lev, :permno)
  mkequityadded::Int = sum((ismissing).(lev.mkequity))
  @mpar for i ∈ 1:length(levgrp)
    slev::SubDataFrame = levgrp[i]

    for sslev ∈ groupby(slev, :compid)
      sslev.last[end] = true

      local Nmissing::Int = sum((ismissing).(sslev.mkequity))
      if Nmissing > 1
        (Nmissing < size(sslev,1)) && error("Inconsistent presence of mkequity")

        dt::Date = eom(sslev.fyrmonth[1])

        #find the first date before the end of the fiscal year if available
        #need to get this from suniv
        local r::Union{Nothing, DataFrameRow} = rowbeforedatesorted(slev, dt)
        (!isnothing(r)) && ( dt - periods2stale ≤ r.date) && (sslev.mkequity .= r.mktcap)
      end
    end
  end

  mkequityadded -= sum((ismissing).(lev.mkequity))
  @info "Supplemented $mkequityadded mkequity records using $crsptype data"

  #now that we have patched the market value, drop all other rows
  #can either free the memory or do this in place
  lev = filter!(r->r.last, lev)

  #make the cash field
  lev.mkat = lev.lt .+ lev.mkequity
  lev.mknegcash = 1.0 .- (lev.cash ./ lev.mkat)

  #compute #flev field
  lev.mkflev = lev.netdebt ./ (lev.netdebt .+ lev.mkequity)

  #compute liab
  lev.mkliab = lev.lt ./ lev.mkat

  #make finite (note we winsorize later on)
  @mpar for s ∈ [:mknegcash, :mkflev, :mkliab]
    lev[!,s] = (f::MFloat64->coalesce(isfinite(f),false) ? f : missing).(lev[!,s])
  end

  sort!(lev, [:permno, :fyear]) #not completely sure we need this sort

  #make lagged change in assets as a control variable
  lev.lat = (log).(lev.at)
  differencewithin!(compa, :lat, :gvkey, :fyear, sorted=true)
  lev.Dat = (exp).(lev.Dlat)

  lev.bm = lev.bkequity ./ lev.mkequity

  lagwithin!(lev, [:vol252d, :vol504d, :mkat, :mkequity, :bm, :Dat], :permno, :fyear, sorted=true)
  lagwithin!(lev, [:Lvol252d, :Lvol504d], :permno, :fyear, sorted=true)
  #lagwithin!(lev, [:LLvol252d, :LLvol504d], :fyear, sorted=true)

  rename!(lev, [:date=>:dailyenddate, :compid=>:compidtvolmcap])
  select!(lev, Not([:permno, :last])) #dup of lpermno

  return lev
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
envpath = joinpath(ENV["HOME"], ".julia/environments")
if isdir(envpath) == false
	mkpath(envpath) # Make the environment directory
	cp("/home/jovyan/.julia/environments", envpath, force = true)
end

max_stars_repo_name: ubcecon/syzygy

max_issues_repo_name: ubcecon/syzygy

max_forks_repo_name: ubcecon/syzygy


text:
@testset "recover camera distortions" begin
    focal_length = 6e-3
    pixelspermeter = 181818.0

    trackprop = TrackProperties(track_gauge = 1.435)

    cameraposition_reference = [0.0,-0.2,-2.4];

    ψθφ_ref = [0.4,-5.0,0.5] .* (pi/180.0);
    camera_reference = VideoCamera(cameraposition_reference;
        focal_length = focal_length,
        pixelspermeter = pixelspermeter,
        ψθφ = ψθφ_ref
    )

    left_track_ref  = [
        [x, -trackprop.track_gauge/2.0, 0.0]
    for x = LinRange(4.0,40.0,40)];

    right_track_ref = [
        [x, trackprop.track_gauge/2.0,0.0]
    for x = LinRange(4.0,40.0,40)];

    left_uvs_ref = camera_image(camera_reference, left_track_ref);
    right_uvs_ref = camera_image(camera_reference, right_track_ref);

    vLs_ref = [uv[2] for uv in left_uvs_ref];
    uLs_ref = [uv[1] for uv in left_uvs_ref];
    vRs_ref = [uv[2] for uv in right_uvs_ref];
    uRs_ref = [uv[1] for uv in right_uvs_ref];

    N = 600;
    dYs = (rand(N) .- 0.5) .* 1e-2;
    dZs = (rand(N) .- 0.5) .* 1e-2;
    dXYZs = [[0,dYs[i],dZs[i]] for i in 1:N];
    XYZ_mean = mean(abs.(d) for d in dXYZs)

    dψθφs = [ [(rand() - 0.5), (rand() - 0.5),0.0] .* (pi/180.0) .* 0.1 for i in 1:N];
    ψθφ_mean = mean(abs.(d) for d in dψθφs)

    # curvatures
    βs = (rand(N) .- 0.5) .* 1e-4 .* 1e-1;
    β_mean = mean(abs.(βs))

    αs = (rand(N) .- 0.5) .* 1e-4.* 1e-1;
    α_mean = mean(abs.(αs))

    δs = [[-dXYZs[i][2:3]; dψθφs[i][1:2]; αs[i]; βs[i]] for i in 1:N];
    δ_mean = [XYZ_mean[2:3]; ψθφ_mean[1:2]; α_mean; β_mean];
    choose_distortions = [:Y,:Z,:ψ,:θ,:α,:β]

    ΔuR = right_track_Δu(camera_reference, trackprop, vRs_ref; choose_distortions = choose_distortions);
    ΔuL = left_track_Δu(camera_reference, trackprop, vLs_ref; choose_distortions = choose_distortions);

    exact_δs = map(1:N) do i

        uLs = uLs_ref + ΔuL * δs[i];
        uRs = uRs_ref + ΔuR * δs[i];

        left_uvs = [ [uLs[j], vLs_ref[j]] for j in eachindex(uLs)];
        right_uvs = [ [uRs[j], vRs_ref[j]] for j in eachindex(uRs)];

        distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop; choose_distortions = choose_distortions)

        [distortion[s] for s in choose_distortions]
    end;

    error_max = maximum(hcat([abs.(d) for d in δs - exact_δs]...), dims=2)
    error_mean = mean(abs.(d) for d in δs - exact_δs)

    @test maximum(error_max ./ δ_mean) < 1e-11
    @test maximum(error_mean ./ δ_mean) < 1e-12

    # use yaw φ instead of roll ψ
    dψθφs = [ [0.0, (rand() - 0.5),(rand() - 0.5)] .* (pi/180.0) .* 0.1 for i in 1:N];
    ψθφ_mean = mean(abs.(d) for d in dψθφs)

    δs = [[-dXYZs[i][2:3]; dψθφs[i][2:3]; αs[i]; βs[i]] for i in 1:N];
    δ_mean = [XYZ_mean[2:3]; ψθφ_mean[2:3]; α_mean; β_mean];
    choose_distortions = [:Y,:Z,:θ,:φ,:α,:β]

    ΔuR = right_track_Δu(camera_reference, trackprop, vRs_ref; choose_distortions = choose_distortions);
    ΔuL = left_track_Δu(camera_reference, trackprop, vLs_ref; choose_distortions = choose_distortions);

    exact_δs = map(1:N) do i

        uLs = uLs_ref + ΔuL * δs[i];
        uRs = uRs_ref + ΔuR * δs[i];

        left_uvs = [ [uLs[j], vLs_ref[j]] for j in eachindex(uLs)];
        right_uvs = [ [uRs[j], vRs_ref[j]] for j in eachindex(uRs)];

        distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop; choose_distortions = choose_distortions)

        [distortion[s] for s in choose_distortions]
    end;

    error_max = maximum(hcat([abs.(d) for d in δs - exact_δs]...), dims=2)
    error_mean = mean(abs.(d) for d in δs - exact_δs)

    @test maximum(error_max ./ δ_mean) < 1e-11
    @test maximum(error_mean ./ δ_mean) < 1e-12

    estimate_δs = map(1:N) do i
        left_track  = [
            [x, -trackprop.track_gauge/2.0 + βs[i] * x^2, αs[i] *x^2]
        for x = LinRange(4.0,40.0,40)];

        right_track = [
            [x, trackprop.track_gauge/2.0 + βs[i] * x^2, αs[i] *x^2]
        for x = LinRange(4.0,40.0,40)];

        camera = VideoCamera(camera_reference.xyz + dXYZs[i], camera_reference.ψθφ + dψθφs[i], camera_reference.opticalproperties)

        left_uvs = camera_image(camera, left_track);
        right_uvs = camera_image(camera, right_track);

        distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop;
            choose_distortions = choose_distortions
            , iterations = 2
        );

        [distortion[s] for s in choose_distortions]
    end;

    error_max = maximum(hcat([abs.(d) for d in δs - estimate_δs]...), dims=2);
    error_mean = mean(abs.(d) for d in δs - estimate_δs);

    # To improve this, I need to add a reference that includes curvature. At the moment, the reference always
    @test maximum(error_max ./ δ_mean) < 0.02
    @test maximum(error_mean ./ δ_mean) < 0.005

    # change angles used
    dψθφs = [ [(rand() - 0.5), (rand() - 0.5),0.0] .* (pi/180.0) .* 0.1 for i in 1:N];
    ψθφ_mean = mean(abs.(d) for d in dψθφs)

    δs = [[-dXYZs[i][2:3]; dψθφs[i][1:2]; αs[i]; βs[i]] for i in 1:N];
    δ_mean = [XYZ_mean[2:3]; ψθφ_mean[1:2]; α_mean; β_mean];
    choose_distortions = [:Y,:Z,:ψ,:θ,:α,:β]

    estimate_δs = map(1:N) do i
        left_track  = [
            [x, -trackprop.track_gauge/2.0 + βs[i] * x^2, αs[i] *x^2]
        for x = LinRange(4.0,40.0,40)];

        right_track = [
            [x, trackprop.track_gauge/2.0 + βs[i] * x^2, αs[i] *x^2]
        for x = LinRange(4.0,40.0,40)];

        camera = VideoCamera(camera_reference.xyz + dXYZs[i], camera_reference.ψθφ + dψθφs[i], camera_reference.opticalproperties)

        left_uvs = camera_image(camera, left_track);
        right_uvs = camera_image(camera, right_track);

        distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop;
            choose_distortions = choose_distortions
            , iterations = 2
        );

        [distortion[s] for s in choose_distortions]
    end;

    error_max = maximum(hcat([abs.(d) for d in δs - estimate_δs]...), dims=2);
    error_mean = mean(abs.(d) for d in δs - estimate_δs);

    # To improve this, I need to add a reference that includes curvature. At the moment, the reference always
    @test maximum(error_max ./ δ_mean) < 0.02
    @test maximum(error_mean ./ δ_mean) < 0.005

end

@testset "recover train distortions" begin
    f = 6e-3
    pixelspermeter = 181818.0

    trackprop = TrackProperties(track_gauge = 1.435)

    cameraposition_reference = [0.0,-0.2,-2.4];

    ψθφ_ref = [0.4,-5.0,0.5] .* (pi/180.0);
    camera_reference = VideoCamera(cameraposition_reference;
        focal_length = f,
        pixelspermeter = pixelspermeter,
        ψθφ = ψθφ_ref
    )

    left_track_ref  = [
        [x, -trackprop.track_gauge/2.0, 0.0]
    for x = LinRange(4.0,40.0,40)];

    right_track_ref = [
        [x, trackprop.track_gauge/2.0,0.0]
    for x = LinRange(4.0,40.0,40)];

    left_uvs_ref = camera_image(camera_reference, left_track_ref);
    right_uvs_ref = camera_image(camera_reference, right_track_ref);

    vLs_ref = [uv[2] for uv in left_uvs_ref];
    uLs_ref = [uv[1] for uv in left_uvs_ref];
    vRs_ref = [uv[2] for uv in right_uvs_ref];
    uRs_ref = [uv[1] for uv in right_uvs_ref];

    N = 600;
    dYs = (rand(N) .- 0.5) .* 1e-2;
    dZs = (rand(N) .- 0.5) .* 1e-2;
    dXYZs = [[0,dYs[i],dZs[i]] for i in 1:N];
    XYZ_mean = mean(abs.(d) for d in dXYZs)

    # dψθφs = [ [(rand() - 0.5), (rand() - 0.5),0.0] .* (pi/180.0) .* 0.1 for i in 1:N];
    dψθφs = [ [0.0, (rand() - 0.5), (rand() - 0.5)] .* (pi/180.0) .* 0.1 for i in 1:N];
    ψθφ_mean = mean(abs.(d) for d in dψθφs)

    # curvatures
    βs = (rand(N) .- 0.5) .* 1e-4 .* 1e-1;
    β_mean = mean(abs.(βs))

    αs = (rand(N) .- 0.5) .* 1e-4.* 1e-1;
    α_mean = mean(abs.(αs))

    δs = [[-dXYZs[i][2:3]; dψθφs[i][2:3]; αs[i]; βs[i]] for i in 1:N];
    δ_mean = [XYZ_mean[2:3]; ψθφ_mean[2:3]; α_mean; β_mean];
    choose_distortions = [:Y,:Z,:θT,:φT,:α,:β]

    ΔuR = right_track_Δu(camera_reference, trackprop, vRs_ref; choose_distortions = choose_distortions);
    ΔuL = left_track_Δu(camera_reference, trackprop, vLs_ref; choose_distortions = choose_distortions);

    exact_δs = map(1:N) do i
        uLs = uLs_ref + ΔuL * δs[i];
        uRs = uRs_ref + ΔuR * δs[i];

        left_uvs = [ [uLs[j], vLs_ref[j]] for j in eachindex(uLs)];
        right_uvs = [ [uRs[j], vRs_ref[j]] for j in eachindex(uRs)];

        distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop; choose_distortions = choose_distortions)

        [distortion[s] for s in choose_distortions]
    end;

    error_max = maximum(hcat([abs.(d) for d in δs - exact_δs]...), dims=2)
    error_mean = mean(abs.(d) for d in δs - exact_δs)

    @test maximum(error_max[[1:3;5:6]] ./ δ_mean[[1:3;5:6]]) < 1e-11
    @test maximum(error_max) < 1e-14

    @test maximum(error_mean[[1:3;5:6]] ./ δ_mean[[1:3;5:6]]) < 1e-12
    @test maximum(error_mean) < 1e-15

    estimate_δs = map(1:N) do i
        left_track  = [
            [x, -trackprop.track_gauge/2.0 + βs[i] * x^2, αs[i] *x^2]
        for x = LinRange(4.0,40.0,40)];

        right_track = [
            [x, trackprop.track_gauge/2.0 + βs[i] * x^2, αs[i] *x^2]
        for x = LinRange(4.0,40.0,40)];

        camera = VideoCamera(camera_reference.xyz + dXYZs[i], camera_reference.ψθφ, camera_reference.opticalproperties)

        traincar = TrainCar(dψθφs[i])

        left_uvs = camera_image(camera, left_track, traincar);
        right_uvs = camera_image(camera, right_track, traincar);

        distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop;
            choose_distortions = choose_distortions
            , iterations = 2
        );

        [distortion[s] for s in choose_distortions]
    end;

    error_max = maximum(hcat([abs.(d) for d in δs - estimate_δs]...), dims=2);
    error_mean = mean(abs.(d) for d in δs - estimate_δs);

    # To improve this, I need to add a reference that includes both the traincar position and curvature.

    @test maximum(error_max[[1:3;5:6]] ./ δ_mean[[1:3;5:6]]) < 0.07
    @test maximum(error_max) < 1e-5

    @test maximum(error_mean[[1:3;5:6]] ./ δ_mean[[1:3;5:6]]) < 0.02
    @test maximum(error_mean) < 1e-5

end


text:
module FESpaceConstructorsTests

using Test
using Gridap

model = CartesianDiscreteModel(partition=(2,2))

V = FESpace(
  reffe=:Lagrangian,
  conformity=:H1,
  valuetype = Float64,
  order = 2,
  model = model)

@test isa(V,CLagrangianFESpace)

V = FESpace(
  reffe=:Lagrangian,
  conformity=:H1,
  valuetype = Float64,
  order = 2,
  model = model,
  diritags = 1)

@test num_diri_dofs(V) == 1
@test isa(V,CLagrangianFESpace)

labels = FaceLabels(model)
add_tag_from_tags!(labels,"diri0",[1,2])

V = FESpace(
  reffe=:Lagrangian,
  conformity=:H1,
  valuetype = Float64,
  order = 1,
  model = model,
  labels = labels,
  diritags = "diri0")

@test num_diri_dofs(V) == 2
@test isa(V,CLagrangianFESpace)

V = FESpace(
  reffe=:Lagrangian,
  conformity=:H1,
  valuetype = Float64,
  order = 2,
  model = model,
  diritags = [1,2])

@test isa(V,CLagrangianFESpace)

V = FESpace(
  reffe=:QLagrangian,
  conformity=:H1,
  valuetype = Float64,
  order = 2,
  model = model)

@test isa(V,CLagrangianFESpace)

V = FESpace(
  reffe=:Lagrangian,
  conformity=:L2,
  valuetype = Float64,
  order = 2,
  model = model)

@test isa(V,DLagrangianFESpace)

V = FESpace(
  reffe=:QLagrangian,
  conformity=:L2,
  valuetype = Float64,
  order = 2,
  model = model)

@test isa(V,DLagrangianFESpace)

V = FESpace(
  reffe=:PLagrangian,
  conformity=:L2,
  valuetype = Float64,
  order = 2,
  model = model)

@test isa(V,DiscFESpace)

V = FESpace(
  reffe=:RaviartThomas,
  conformity=:HDiv,
  order = 2,
  model = model)

@test isa(V,ConformingFESpace)

V = FESpace(
  reffe=:RaviartThomas,
  conformity=:L2,
  order = 2,
  model = model)

@test isa(V,DiscFESpace)

V = FESpace(
  reffe=:Nedelec,
  conformity=:HCurl,
  order = 2,
  model = model)

@test isa(V,ConformingFESpace)

V = FESpace(
  reffe=:Nedelec,
  conformity=:L2,
  order = 2,
  model = model)

@test isa(V,DiscFESpace)

V = FESpace(
  reffe=:PLagrangian,
  conformity=:L2,
  valuetype = Float64,
  order = 2,
  model = model,
  constraint = :zeromean)

@test isa(V,ZeroMeanFESpace)

end # module


text:
# To use this package
using CapacityExpansion

#########################
#= Load Time Series Data
#########################
How to load data provided with the package:
The data is for a Capacity Expansion Problem "CEP"
and for the single node representation of Germany "GER_1"
The original timeseries has 8760 entries (one for each hour of the year)
It should be cut into K=365 periods (365 days) with T=24 timesteps per period (24h per day) =#
ts_input_data = load_timeseries_data_provided("GER_1"; T=24)

#= ClustData
How the struct is setup:
    ClustData{region::String,K::Int,T::Int,data::Dict{String,Array},weights::Array{Float64},mean::Dict{String,Array},sdv::Dict{String,Array}} <: TSData
-region: specifies region data belongs to
-K: number of periods
-T: time steps per period
-data: Data in form of a dictionary for each attribute `"[file name]-[column name]"`
-weights: this is the absolute weight. E.g. for a year of 365 days, sum(weights)=365
-mean: For normalized data: The shift of the mean as a dictionary for each attribute
-sdv: For normalized data: Standard deviation as a dictionary for each attribute

How to access a struct:
    [object].[fieldname]                                                      =#
number_of_periods=ts_input_data.K
# How to access a dictionary:
data_solar_germany=ts_input_data.data["solar-germany"]
# How to plot data
using Plots
# plot(Array of our data, no legend, dotted lines, label on the x-Axis, label on the y-Axis)
plot_input_solar=plot(ts_input_data.data["solar-germany"], legend=false, linestyle=:dot, xlabel="Time [h]", ylabel="Solar availability factor [%]")

# How to load your own data:
# put your data into your homedirectory into a folder called tutorial
# The data should have the following structure: see TimeSeriesClustering/data folder
#=
- Loading all `*.csv` files in the folder or the file `data_path`
The `*.csv` files shall have the following structure and must have the same length:
|Timestamp |[column names...]|
|[iterator]|[values]         |
The first column should be called `Timestamp` if it contains a time iterator
The other columns can specify the single timeseries like specific geolocation.
Each column in `[file name].csv` file will be added to the ClustData.data called `"[file name]-[column name]"`
- region is an additional String to specify the loaded time series data
- K describes the number of periods in the input data
- T describes the length of each period                             =#
load_your_own_data=false
if load_your_own_data
    # Single file at the path e.g. homedir/tutorial/solar.csv
    # It will automatically call the data 'solar' within the datastruct
    my_path=joinpath(homedir(),"tutorial","solar.csv")
    your_data_1=load_timeseries_data(my_path; T=24)
    # Multiple files in the folder e.g. homedir/tutorial/
    # Within the data struct, it will automatically call the data the names of the csv filenames
    my_path=joinpath(homedir(),"tutorial")
    your_data_2=load_timeseries_data(my_path; T=24)
end

#############
# Clustering
#############
# Quick example and investigation of the best result:
ts_clust_result = run_clust(ts_input_data; method="kmeans", representation="centroid", n_init=5, n_clust=5) # note that you should use n_init=1000 at least for kmeans.
ts_clust_data = ts_clust_result.clust_data
# And some plotting:
plot_comb_solar=plot!(plot_input_solar, ts_clust_data.data["solar-germany"], linestyle=:solid, width=3)
plot_clust_soar=plot(ts_clust_data.data["el_demand-germany"], legend=false, linestyle=:solid, width=3, xlabel="Time [h]", ylabel="Solar availability factor [%]")

#= Clustering options:
`run_clust()` takes the full `data` and gives a struct with the clustered data as the output.

## Supported clustering methods
The following combinations of clustering method and representations are supported by `run_clust`:

Name                                                | method            | representation
----------------------------------------------------|-------------------|----------------
k-means clustering                                  | `<kmeans>`        | `<centroid>`
k-means clustering with medoid representation       | `<kmeans>`        | `<medoid>`
k-medoids clustering (partitional)                  | `<kmedoids>`      | `<medoid>`
k-medoids clustering (exact) [requires Gurobi]      | `<kmedoids_exact>`| `<medoid>`
hierarchical clustering with centroid representation| `<hierarchical>`  | `<centroid>`
hierarchical clustering with medoid representation  | `<hierarchical>`  | `<medoid>`

## Other input parameters

The input parameter `n_clust` determines the number of clusters,i.e., representative periods.

`n_init` determines the number of random starting points. As a rule of thumb, use:
`n_init` should be chosen 1000 or 10000 if you use k-means or k-medoids
`n_init` should be chosen 1 if you use k-medoids_exact or hierarchical clustering

`iterations` is defaulted to 300, which is a good value for kmeans and kmedoids in our experience. The parameter iterations does not matter when you use k-medoids exact or hierarchical clustering.

=#

# A clustering run with different options chosen as an example
ts_clust_result_2 = run_clust(ts_input_data; method="kmedoids", representation="medoid", n_init=100, n_clust=4, iterations=500)

######
# CEP
######
# Using a Solver called Clp (if not installed run `using Pkg; Pkg.add("Clp")`):
using Clp
optimizer=Clp.Optimizer
# Some extra data for nodes, costs and so on:
cep_data = load_cep_data_provided(ts_clust_data.region)
# Running a simple CEP with a co2-limit of 1000 kg/MWh
co2_result = run_opt(ts_clust_data,cep_data,optimizer;descriptor="co2",limit_emission=Dict{String,Number}("CO2/electricity"=>200))
# co2_result.
gen_var=co2_result.variables["GEN"]
# get all operating decision variables (GEN)
gen_var[:,:,:,:,:]

# get specific operating decision variables (the one for the electricity sector, the coal technology at the node germany) indicated with names
coal_gen_var=gen_var["wind","electricity",:,:,"germany"]
# plot the specific operation decision variables
plot(axes(gen_var,"time_T"),coal_gen_var, xlabel="Time [h]", ylabel="Wind Generation [MW]", label=axes(gen_var,"time_K"), legendtitle="Period K", linewidth=2)