text:
if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open OpMul.jl")
TabLevel=TabLevel*"    "

import Base.*
import Base.^

*(x::MajTypes,y::MajTypes) = OpMul(x,y)
*(x::Scalars,y::MajTypes) = OpScale(x,y)
*(x::MajTypes,y::Scalars) = OpScale(y,x)
^(x::MajTypes,y::Integer) = OpPower(x,y)

function OpMul(A) return OpCopy(A) end ##Of only one factor then the result is trivial

function OpMul(a::MajoranaSp,b::MajoranaSp)
    ###When both terms are superpositoins we add them term by term

    tot=OpScale(0)
    #print("tot=")
    #OpPrintln(tot)
    #println("tot=",tot)
    #print("b=")
    #OpPrintln(b)
    for aOP in a.SpList
        #print("aOP=")
        #OpPrintln(aOP)
        aOPb=aOP*b
        #print("aOPb=")
        #OpPrintln(aOPb)
        tot=aOPb+tot
        #print("tot=")
        #OpPrintln(tot)
    end
    return tot
end

function OpMul(a::MajoranaOp,b_in::MajoranaSp)
    if a == 0
        return 0
    end
    b = OpCopy(b_in)
    for indx in 1:length(b.SpList)
        b.SpList[indx]=a*b.SpList[indx]
    end
    return OpReduce(b)

end

function OpMul(a_in::MajoranaSp,b::MajoranaOp)
    if b == 0
        return 0
    end
    a = OpCopy(a_in)
    ###Loop over the term an multiply the on by one
    for indx in 1:length(a.SpList)
        a.SpList[indx]=a.SpList[indx]*b
    end
    return OpReduce(a)
end

function OpMul(a::MajoranaOp,b::MajoranaOp)
    ##Check that a is not a ket state
    if a.KetState
        error("The state "*OpToStr(a)*" may not be to the left of a matrix multiplication ")
    end
    OpReduce(MajoranaOp(a.scale*b.scale,vcat(a.OpList,b.OpList),b.KetState))
end

function OpMul(a::MajElem,b::MajElem)
    OpReduce(MajoranaOp(1,[a,b]))
end

OpMul(a::MajElem,b::MajoranaOp) = MajoranaOp(1,[a])*b
OpMul(b::MajoranaOp,a::MajElem) = b*MajoranaOp(1,[a])

###  Methods to Make majoranas equal
#### Scaling operators
function OpScale(scale::Scalars)
    return scale
end

function OpScale(scale::Scalars,a_in::MajoranaSp;verbose=false)
    verbose && println("Scale $a_in by  $scale")
    if a_in.Reduced ###Reduce the Majorana first
        a = OpCopy(a_in)
    else
        a = OpReduce(a_in,verbose=verbose)
    end
    for indx in 1:length(a.SpList) ###Updated the terms one by one
        a.SpList[indx]=OpScale(scale,a.SpList[indx],verbose=verbose)
    end
    return a
end

function OpScale(scale::Scalars,a_in::MajoranaOp;verbose=false)
    verbose && println("Scale $a_in by  $scale")
    if scale == 0
        return 0
    end
    a = OpReduce(a_in)
    a.scale*=scale
    return a
end

function OpScale(scale::Scalars,a::MajElem;verbose=false)
    verbose && println("Scale $a_in by  $scale")
    if scale == 0
        return 0
    else
        return MajoranaOp(scale,[a])
    end
end

function OpPower(x::Union{MajoranaOp,MajoranaSp,MajElem},y::Integer)
    if y < 0
        throw(DomainError("Majorana cannot be raised to a negative power $y"))
    elseif y==0
        return 1
    elseif y==1
        return 1*x ###To make sure it copied and reduced
    elseif y==2
        return x*x ###To make sure it copied and reduced
    else
        return x*(x^(y-1)) ###Recursively apply the product
    end
end

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close OpMul.jl")


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    PredicatePartition(predicate)

A method for partitioning spatial objects with a given `predicate`
function. Two locations `i` and `j` are part of the same subset
whenever `predicate(i, j) == true`
"""
struct PredicatePartition <: PredicatePartitionMethod
  pred::Function
end

(p::PredicatePartition)(i, j) = p.pred(i, j)

ext: jl


text:
include("lib_for_robot.jl")
function task13(r)
    putmarker!(r)
    side=(Nord,Ost)
    x,y = move_diog(r,side)
    rtn_diog(r,side,x,y)
    side=(West,Nord)
    x,y = move_diog(r,side)
    rtn_diog(r,side,x,y)
    side=(Sud,West)
    x,y = move_diog(r,side)
    rtn_diog(r,side,x,y)
    side=(Sud,Ost)
    x,y = move_diog(r,side)
    rtn_diog(r,side,x,y)
end


text:
module Day6

TESTINPUT = """3,4,3,1,2"""

struct LanternFish
    timer::Int64
end

LanternFish(s::String) = LanternFish.(parse.(Int64, split(s, ",")))

age(l::LanternFish) = l.timer == 0 ? LanternFish(6) : LanternFish(l.timer - 1)
spawn(l::LanternFish) = l.timer == 0 ? LanternFish(8) : nothing
step_day(ls::Vector{LanternFish}) = vcat(age.(ls), filter(!isnothing, spawn.(ls)) |> x -> convert(Vector{LanternFish}, x))
simulate(ls::Vector{LanternFish}, days::Int64) = days == 0 ? ls : simulate(step_day(ls), days - 1)

run_part1(data_dir::String) = joinpath(data_dir, "day6.txt") |>
    f -> read(f, String) |>
    LanternFish |>
    ls -> simulate(ls, 80) |>
    length

mutable struct LanternFishEcosystem
    state::Dict{Int64, Int64}
end

function LanternFishEcosystem(ls::Vector{LanternFish})
    state = Dict(i => 0 for i in 0:8)
    for l in ls
        state[l.timer] += 1
    end
    LanternFishEcosystem(state)
end

function step_day(eco::LanternFishEcosystem)
    Dict(7 => eco.state[8],
         6 => eco.state[7] + eco.state[0],
         5 => eco.state[6],
         4 => eco.state[5],
         3 => eco.state[4],
         2 => eco.state[3],
         1 => eco.state[2],
         0 => eco.state[1],
         8 => eco.state[0]) |> LanternFishEcosystem
end
simulate(ls::LanternFishEcosystem, days::Int64) = days == 0 ? ls : simulate(step_day(ls), days - 1)

Base.length(eco::LanternFishEcosystem) = values(eco.state) |> sum

run_part2(data_dir::String) = joinpath(data_dir, "day6.txt") |>
    f -> read(f, String) |>
    LanternFish |>
    LanternFishEcosystem |>
    eco -> simulate(eco, 256) |>
    length

end # module


text:
"""
Reproduction of the empirical results concerning an empirical application
in the prediction of transportation mode choice using the best subset
maximum score cross validation approach of Chen and Lee (2017).
The work-trip mode choice dataset of Horowitz (1993) is included in the package.
"""
module transportationModeCV

using DataFrames, Compat
using MathProgBase, Gurobi
using ..BSBP

""" Best subset maximum score approach with cross validation on the work-trip mode choice dataset of Horowitz.
Args:
    warm_start (int,optional):  Set to 1 for warm start strategy. Default = 1.
    tau (float,optional):       Tuning paramater to construct the refined bound used in the warm start approach. Default = 1.5.
    mio (int,optional):         Set to 1 for using Method 1 for the MIO formulation. Set to 2 for Method 2. Default to 1.
    q (int, optional):          Value of the variable selection bound. Default = 1.
    series_exp (int, optional): Set to 1 to use quadratic expansion terms as covariates. Default = 1.
    b (int, optional):          Bound value. Default = 1.
    time_limit (int, optional): MIO solver time limit. Default = 8640000.
"""
function TransportationModeCV(;warm_start=1,tau=1.5,mio=1,series_exp=0,beta0=1,b=10,q=1,time_limit=86400)
    data = readdlm("./data_horowitz.csv",',')
    tr_ind = readdlm("./tr_ind.csv",',') .== "TRUE"
    test_ind = readdlm("./test_ind.csv",',') .== "TRUE"

    # create variables from training and validation samples
    # data columins : [Y DCOST CARS DOVTT DIVTT]

    fold=size(tr_ind)[2]

    if series_exp==1
        bhat=zeros(10,fold)
    else
        bhat=zeros(4,fold)
    end

    score=zeros(fold)
    gap=zeros(fold)
    in_score=zeros(fold)
    rtime=zeros(fold)
    ncount=zeros(fold)
    p_ratio=zeros(fold)

    for i=1:fold
        data_tr=data[tr_ind[:,i],:]
        data_v=data[test_ind[:,i],:]

        println("estimation based on training sample at fold: ",i)
        Y_tr=data_tr[:,1]
        temp=data_tr[:,2:end]

        n_tr=length(Y_tr)
        # [DCOST CARS DOVTT DIVTT]
        x_std=(temp.-repmat(mean(temp,1),n_tr,1))./repmat(std(temp,1),n_tr,1)
        x_foc = [x_std[:,1] ones(n_tr,1)]  # [DCOST Intercept]

        if series_exp == 1
            z2 = x_std[:,2]
            z3 = x_std[:,3]
            z4 = x_std[:,4]
            x_aux1 = [z2 z3 z4] # linear terms
            x_aux2 = [z2.*z3 z3.*z4 z2.*z4]
            x_aux3 = [z2.*z2 z3.*z3 z4.*z4]
            x_aux = [x_aux1 x_aux2 x_aux3]
        else
            x_aux = x_std[:,2:4] # [CARS DOVTT DIVTT]
        end
        k=size(x_foc)[2]
        d=size(x_aux)[2]

        bnd= [-b*ones(k-1+d,1) b*ones(k-1+d,1)] # set the initial parameter bounds

        tol = floor(sqrt(log(n_tr)*n_tr)/2)
        println("tolerance level: ", tol)

        if warm_start == 1 # warm start MIO
            bhat[:,i],score[i],gap[i],rtime[i],ncount[i]  = WarmStartMaxScore(Y_tr,x_foc,x_aux,beta0,q,time_limit,tol,bnd,mio,tau)
        else # cold start MIO
            bhat[:,i],score[i],gap[i],rtime[i],ncount[i]  = MaxScore(Y_tr,x_foc,x_aux,beta0,q,time_limit,tol,bnd,mio)
        end

        println("coefficient values: ")
        println("gurobi score: ",score[i])
        println("gurobi absolute gap: ",gap[i])
        println("gurobi running time: ",rtime[i])
        println("gurobi node count: ",ncount[i])
        in_score[i]=sum(Int.(Y_tr.==(([x_foc x_aux]*[1;bhat[:,i]]).>0)))
        println("in-sample score: ",in_score[i])

        # validation sample

        Y_val=data_v[:,1]
        n_val=length(Y_val)
        temp=data_v[:,2:size(data_v)[2]]
        x_std=(temp.-repmat(mean(temp,1),n_val,1))./repmat(std(temp,1),n_val,1)

        if series_exp == 1
            z2 = x_std[:,2]
            z3 = x_std[:,3]
            z4 = x_std[:,4]
            x_aux1 = [z2 z3 z4] # linear terms
            x_aux2 = [z2.*z3 z3.*z4 z2.*z4]
            x_aux3 = [z2.*z2 z3.*z3 z4.*z4]
            x_aux = [x_aux1 x_aux2 x_aux3]
            x_v = [x_std[:,1] ones(n_val,1) x_aux]
        else
            x_v = [x_std[:,1] ones(n_val,1) x_std[:,2:4]]

        end

        y_hat= Int.((x_v*[1;bhat[:,i]]).>0)
        p_ratio[i]=mean(Int.(Y_val.==y_hat))
        println("out-of-sample performance: ", p_ratio[i])
    end

    println("Average coefficient vector: ",mean(bhat,2))
    println("Average score: ",mean(score))
    println("Average gap: ",mean(gap))
    println("Average running time: ",mean(rtime))
    println("Average node count: ",mean(ncount))
    println("Average in-sample score: ",mean(in_score))
    println("average out-of-sample performance: ",mean(p_ratio))
    return score,gap,rtime,ncount
end

export TransportationModeCV
end


text:
"""
```
function impulse_responses(m::AbstractDSGEVARModel{S}, data::AbstractArray{S}, method::Symbol,
                           n_obs_shock::Int; horizon::Int = 0,
                           flip_shocks::Bool = false, verbose::Symbol = :none) where {S <: Real}

function impulse_responses(m::AbstractDSGEVARModel{S}, method::Symbol,
                           n_obs_shock::Int; horizon::Int = 0, use_intercept::Bool = false,
                           flip_shocks::Bool = false, verbose::Symbol = :none) where {S <: Real}

```
computes the VAR impulse responses identified by the DSGE
```
sₜ = TTT × sₜ₋₁ + RRR × impact[:, i],
yₜ = ZZ × sₜ + DD + MM × impact[:, i],
```
where `impact[:, i]` is a linear combination of
(orthogonal) structural shocks `ϵₜ ∼ 𝒩 (0, I)`, and
`MM × impact[:, i]` are the correlated measurement errors.

We draw a β and Σᵤ from the posterior implied by the DSGE
and data, and we then compute normal VAR impulse responses given those
coefficients and innovations variance-covariance matrix. The weight placed
on the DSGE is encoded by the field `λ` of the DSGEVAR object `m`.

Given β, Σᵤ, we compute impulse responses to the VAR system
```
ŷₜ₊₁ = X̂ₜ₊₁β + uₜ₊₁,
```
where `X̂ₜ₊₁` are the lags of observables in period `t + 1`, i.e. `yₜ, yₜ₋₁, ..., yₜ₋ₚ₊₁`
using one of the available identification methods for VARs

If the second function is used (where `data` is not an input), then we assume
the user wants to compute the VAR approximation of the DSGE,
regardless of the `λ` value in `m`. Note that this function will not
update the value of `λ` in `m` (even though we are computing the DSGE-VAR(∞) approximation).

### Inputs
* `method::Symbol`: The available methods are `:cholesky`, `:maxBC`, and `:choleskyLR`.
    See the docstrings `impulse_responses` for VARs specifically.
* `n_obs_shock::Int`: The index of the observable corresponding to the orthogonalized shock
    causing the impulse response.

### Keyword Arguments
* `horizon::Int`: the desired horizon of the impulse responses.
* `use_intercept::Bool`: use an intercept term for the VAR approximation
* `flip_shocks::Bool`: default is a "negative" impulse response on impact.
    Set to `true` for the positive impulse response.
"""
function impulse_responses(m::AbstractDSGEVARModel{S}, data::AbstractArray{S}, method::Symbol,
                           n_obs_shock::Int; horizon::Int = impulse_response_horizons(m),
                           flip_shocks::Bool = false, frequency_band::Tuple{S, S} = (2*π/32, 2*π/6),
                           verbose::Symbol = :none) where {S <: Real}
    β, Σ = compute_system(m, data; verbose = verbose)
    Σ += Σ'
    Σ ./= 2

    return impulse_responses(β, Σ, n_obs_shock, horizon;
                             method = method, flip_shocks = flip_shocks,
                             use_intercept = true, frequency_band = frequency_band)
end

function impulse_responses(m::AbstractDSGEVARModel{S}, method::Symbol,
                           n_obs_shock::Int; horizon::Int = 0, use_intercept::Bool = false,
                           flip_shocks::Bool = false, frequency_band::Tuple{S, S} = (2*π/32, 2*π/6),
                           verbose::Symbol = :none) where {S <: Real}
    β, Σ = compute_system(m; verbose = verbose, use_intercept = use_intercept)
    Σ += Σ'
    Σ ./= 2

    return impulse_responses(β, Σ, n_obs_shock,
                             horizon > 0 ? horizon : impulse_response_horizons(m);
                             method = method, use_intercept = use_intercept,
                             frequency_band = frequency_band,
                             flip_shocks = flip_shocks)
end

"""
```
function impulse_responses(m::AbstractDSGEVARModel{S}, data::AbstractArray{S},
                           X̂::Matrix{S} = Matrix{S}(undef, 0, 0);
                           horizon::Int = 0, MM::Matrix{S} = Matrix{S}(undef, 0, 0),
                           flip_shocks::Bool = false, draw_shocks::Bool = false,
                           deviations::Bool = false, normalize_rotation::Bool = false,
                           verbose::Symbol = :none) where {S <: Real}
```
computes the VAR impulse responses identified by the DSGE
```
sₜ = TTT × sₜ₋₁ + RRR × impact[:, i],
yₜ = ZZ × sₜ + DD + MM × impact[:, i],
```
where `impact[:, i]` is a linear combination of
(orthogonal) structural shocks `ϵₜ ∼ 𝒩 (0, I)`, and
`MM × impact[:, i]` are the correlated measurement errors.

The VAR impulse responses are computed according to
```
ŷₜ₊₁ = X̂ₜ₊₁β + uₜ₊₁,
```
where `X̂ₜ₊₁` are the lags of observables in period `t + 1`, i.e. `yₜ, yₜ₋₁, ..., yₜ₋ₚ`.
Note these impulse responses are *not* computed in deviations
from the baseline forecast `ŷₜ₊₁ = X̂ₜ₊₁β`. To compute these
impulse responses, use the keyword `deviations`.

The shock `uₜ₊₁` is identified by assuming
```
Σᵤ = 𝔼[u × u'] = chol(Σᵤ) × Ω × ϵₜ,
```
where the rotation matrix `Ω` is the `Q` matrix from a QR decomposition
of the impact response matrix corresponding to the state space system, i.e.
```
Ω, _ = qr(∂yₜ / ∂ϵₜ').
```

For reference, see Del Negro and Schorfheide (2004),
Del Negro and Schorfheide (2006), and Del Negro and Schorfheide (2009).

### Inputs
* `X̂::Matrix{S}`: covariates for the first "forecast" period
    of the impulse response, i.e. if we have a VAR with `p` lags, then
```
X̂ = [1, ŷₜ, ŷₜ₋₁, ..., ŷₜ₋ₚ₊₁]
```
so that, when β is the vector of VAR coefficients, then
```
𝔼[ŷₜ₊₁] = kron(I, X̂') * β.
```
Internally, we do equivalent matrix operations to avoid allocating
the Kronecker product.

****
NOTE: this function generally involves taking random draws from
probability distributions, so seeds need to be set
to achieve reproducibility.
****

### Keywords
* `horizon::Int`: horizon of impulse responses
* `flip_shocks::Bool`: impulse response shocks are negative by default. Set to `true` for
    a positive signed shock.
* `draw_shocks::Bool`: true if you want to draw shocks along the entire horizon
* `deviations::Bool`: set true to compute the impulse response in deviations
    rather than as a forecast. Mechnically, we ignore `X̂` (treated as zeros)
    and the intercept term.
* `normalize_rotation::Bool`: set to true to normalize the rotation
    so that rows have the correct sign. Requires same number of structural shocks
    as observables.
* `verbose::Symbol`: quantity of output desired
"""
function impulse_responses(m::AbstractDSGEVARModel{S}, data::AbstractArray{S},
                           X̂::Vector{S} = Vector{S}(undef, 0);
                           horizon::Int = 0, MM::Matrix{S} = Matrix{S}(undef, 0, 0),
                           flip_shocks::Bool = false, draw_shocks::Bool = false,
                           deviations::Bool = false, normalize_rotation::Bool = false,
                           verbose::Symbol = :none) where {S <: Real}

    # Prepare X̂
    nobs = size(data, 1)
    k = nobs * get_lags(m) + 1
    if deviations
        X̂ = zeros(S, k)
    elseif isempty(X̂)
        XX = lag_data(data, get_lags(m); use_intercept = true)
        X̂ = vcat(1, data[:, end], XX[end, 1+1:k - nobs])
    end
    h = (horizon > 0) ? horizon : impulse_response_horizons(m) # grab horizons

    # Compute underlying state space system
    system = compute_system(m; get_system = true, use_intercept = true)

    # Compute VAR coefficients
    β, Σ = compute_system(m, data; verbose = verbose)
    Σ += Σ'
    Σ ./= 2.

    # Use rotation identification
    if isempty(MM)
        if hasmethod(measurement_error, (typeof(m),))
            _, MM = measurement_error(m)
        else
            MM = zeros(S, n_observables(m), n_shocks(m))
        end
    end
    return impulse_responses(system[:TTT], system[:RRR], system[:ZZ], system[:DD], MM,
                             system[:QQ], k, β, Σ, h, X̂; flip_shocks = flip_shocks,
                             draw_shocks = draw_shocks, deviations = deviations,
                             normalize_rotation = normalize_rotation)
end

"""
```
function impulse_responses(TTT::Matrix{S}, RRR::Matrix{S}, ZZ::Matrix{S},
                           DD::Vector{S}, MM::Matrix{S}, impact::Matrix{S}, horizon::Int;
                           accumulate::Bool = false,
                           cum_inds::Union{Int,UnitRange{Int},Vector{Int}} = 0) where {S<:Real}
```
computes the impulse responses of the linear state space system to linear combinations
of (orthogonal) structural shocks specified by `impact`.
Measurement error that is correlated with
the impact matrix is allowed. We also include the option to accumulate certain observables.

This state space model takes the form

```
sₜ = TTT × sₜ₋₁ + RRR × impact[:, i],
yₜ = ZZ × sₜ + DD + MM × impact[:, i],
```
where `impact[:, i]` is a linear combination of orthogonal structural shocks
with mean zero and identity covariance, and `MM × impact[:, i]` are the
correlated measurement errors.

The `impact` matrix must be `nshock × nirf`, where `nshock` is
 the number of structural shocks and `nirf` is the number
of desired impulse responses. For each row of `impact`,
we compute the corresponding impulse responses.

A standard choice for `impact` is a square diagonal matrix. In this case,
we compute the impulse response of observables to each structural shock,
scaled by the desired size of the shock.

### Keywords
* `accumulate`: set to true if an observable should be accumulated.
* `cum_inds`: specifies the indices of variables to accumulated.

### Outputs
* `irf_results::Matrix{S}`: a `nobs x horizon × nirf` matrix, where
     `nobs` is the number of observables.
"""
function impulse_responses(TTT::Matrix{S}, RRR::Matrix{S}, ZZ::Matrix{S},
                           DD::Vector{S}, MM::Matrix{S}, impact::Matrix{S}, horizon::Int;
                           accumulate::Bool = false,
                           cum_inds::Union{Int,UnitRange{Int},Vector{Int}} = 0) where {S<:Real}
    # Get dimensions
    nshock, nirf = size(impact)
    nstate       = size(TTT, 1)
    nobs         = size(ZZ, 1)

    # Compute impulse response to identified impact matrix
    irf_results = Array{S, 3}(undef, nobs, horizon, nirf)
    for i = 1:nirf
        imp    = impact[:, i]
        states = zeros(S, nstate, horizon)
        obs    = zeros(S, nobs, horizon)

        states[:, 1] = RRR * imp
        obs[:, 1]    = ZZ * states[:, 1] + MM * imp
        for t = 2:horizon
            states[:, t] = TTT * states[:, t - 1]
            obs[:, t]    = ZZ * states[:, t] + DD
        end
        if accumulate
            obs[cum_inds, :] = cumsum(obs[cum_inds, :], dims = length(cum_inds) > 1 ? 2 : 1)
        end
        irf_results[:, :, i] = obs
    end

    return irf_results
end

"""
```
function impulse_responses(TTT::Matrix{S}, RRR::Matrix{S}, ZZ::Matrix{S},
                           DD::Vector{S}, MM::Matrix{S}, QQ::Matrix{S},
                           k::Int, β::Matrix{S}, Σ::Matrix{S},
                           horizon::Int, X̂::Matrix{S} = zeros(S, k);
                           flip_shocks::Bool = false, draw_shocks::Bool = false,
                           deviations::Bool = false, normalize_rotation::Bool = false,
                           test_shocks::Matrix{S} =
                           Matrix{S}(undef, 0, 0)) where {S<:Real}
```
computes the VAR impulse responses identified by the state space system
```
sₜ = TTT × sₜ₋₁ + RRR × ϵₜ
yₜ = ZZ × sₜ + DD + MM × ϵₜ
```
where `ϵₜ ∼ 𝒩 (0, QQ)` and `MM × ϵₜ` are the correlated measurement errors.

The VAR impulse responses are computed according to
```
ŷₜ₊₁ = X̂ₜ₊₁β + uₜ₊₁,
```
where `X̂ₜ₊₁` are the lags of observables in period `t + 1`, i.e. `yₜ, yₜ₋₁, ..., yₜ₋ₚ`.
Note these impulse responses are *not* computed in deviations
from the baseline forecast `ŷₜ₊₁ = X̂ₜ₊₁β`. To compute these
impulse responses, set the keyword `deviations = true`.

The shock `uₜ₊₁` is identified via
```
Σᵤ = 𝔼[u × u'] = chol(Σᵤ) × Ω × ϵₜ,
```
where the rotation matrix `Ω` is the `Q` matrix from a QR decomposition
of the impact response matrix corresponding to the state space system, i.e.
```
Ω, _ = qr(∂yₜ / ∂ϵₜ').
```
To normalize the rotation so that rows have the correct sign,
set `normalize_rotation = true`. This requires that there are
as many structural shocks as observables.

For reference, see Del Negro and Schorfheide (2004),
Del Negro and Schorfheide (2006), and Del Negro and Schorfheide (2009).
"""
function impulse_responses(TTT::Matrix{S}, RRR::Matrix{S}, ZZ::Matrix{S},
                           DD::Vector{S}, MM::Matrix{S}, QQ::Matrix{S},
                           k::Int, β::Matrix{S}, Σ::Matrix{S},
                           horizon::Int, X̂::Vector{S} = zeros(S, k);
                           flip_shocks::Bool = false, draw_shocks::Bool = false,
                           deviations::Bool = false, normalize_rotation::Bool = false,
                           test_shocks::Matrix{S} =
                           Matrix{S}(undef, 0, 0)) where {S <: Real}

    # Compute impulse responses of predicted values for each β, Σ, and rotation
    nobs = size(ZZ, 1)
    a0_m = convert(Matrix{S},
                   dropdims(impulse_responses(TTT, RRR, ZZ, DD, MM,
                                              sqrt.(QQ), 1, # 1 b/c just want impact and sqrt -> 1 sd shock
                                              accumulate = false); dims = 2)')
    rotation, r_a = qr(a0_m)
    if normalize_rotation
        rotation = sign.(diag(r_a)) .* rotation' # normalizes each row (change signs) so that lower diagonal (r_a') has positive diagonal elements
        Σ_chol = cholesky(Σ).L * rotation # mapping from structural shocks to innovations in VAR
    else
        Σ_chol = cholesky(Σ).L * rotation'
    end

    # Ignore intercept if in deviations
    if deviations
        β = β[2:end, :]
        X̂ = zeros(S, size(β, 1))
        k -= 1
    end

    if draw_shocks || !isempty(test_shocks)
        ŷ = Matrix{S}(undef, nobs, horizon)

        if isempty(test_shocks)
            shocks = randn(size(RRR, 2), horizon) # shocks getting drawn are structural shocks
            if flip_shocks
                @warn "Shocks are being drawn, so flipping shocks does not make sense in this context. Ignoring `flip_shocks` . . ."
            end
        else
            shocks = test_shocks
        end

        for t = 1:horizon
            out     = vec(X̂' * β) + Σ_chol * shocks[:, t] # X̂ normally would be [X̂ 0 0; 0 X̂ 0; 0 0 X̂] if nobs = 3, but this way of coding it results in less memory storage
            ŷ[:, t] = out

            X̂       = deviations ? vcat(out, X̂[1:k - nobs]) :
                vcat(1., out, X̂[1 + 1:k - nobs]) # XXl = X̂[1 + 1:k - nobs]
        end
    else
        nshocks = size(RRR, 2)
        ŷ       = Array{S, 3}(undef, nobs, horizon, nshocks)
        old_X̂   = X̂
        shocks  = zeros(S, nshocks)

        for i = 1:nshocks
            X̂ = copy(old_X̂)
            shocks[i] = flip_shocks ? sqrt(QQ[i, i]) :
                -sqrt(QQ[i, i]) # a negative 1 s.d. shock by default
            out        = vec(X̂' * β) + Σ_chol * shocks # do impact separately
            shocks[i]  = 0. # set back to zero
            ŷ[:, 1, i] = out
            X̂          = deviations ? vcat(out, X̂[1:k - nobs]) :
                vcat(1., out, X̂[1 + 1:k - nobs]) # XXl = X̂[1 + 1:k - nobs]
            for t = 2:horizon
                out        = vec(X̂' * β)
                ŷ[:, t, i] = out
                X̂          = deviations ? vcat(out, X̂[1:k - nobs]) :
                    vcat(1., out, X̂[1 + 1:k - nobs]) # XXl = X̂[1 + 1:k - nobs]
            end
        end
    end

    return ŷ
end


text:
module FourierCollocation

using ..NumericalContinuation: Var, Data, Func, Problem
using ..NumericalContinuation: parameter, parameters, parameter_names
using DocStringExtensions

export fourier_collocation

struct FourierColl{F}
    f!::F
    n_dim::Int64
end

# TODO: allow for out-of-place functions (parameterise FourierColl further?)

function (fourier::FourierColl)(res, (u, p, t), data)
    mul!(data.Du, reshape(u, (fourier.n_dim, data.n_mesh)), data.D)
    ii = 1:(fourier.n_dim)
    T = t[2] - t[1]
    for i in 1:(data.n_mesh)
        @views fourier.f!(res[ii], u[ii], p, t[1] + T * (i - 1) / data.n_mesh)
        @views res[ii] .= data.Du[ii] .- T .* res[ii]
        ii = ii .+ fourier.n_dim
    end
    return nothing
end

function fourier_collocation(
    name::String,
    f!,
    trange,
    u0,
    p0;
    pnames=nothing,
    top_level=true,
    phase=true,
    fix_t0=true,
    fix_t1=false,
)
    if length(trange) == 1
        t0 = [0, trange[1]]
    elseif length(trange) == 2
        t0 = trange
    else
        throw(ArgumentError("Expected trange to contain start and end times only"))
    end
    # Generate parameter names
    _pnames = parameter_names(p0, pnames)
    # Create the continuation problem
    problem = Problem(name)
    # Create differentiation matrix and temporary storage to avoid allocations
    D = -fourier_diff(eltype(u0), size(u0, 2)) * 2π
    Du = similar(u0)
    # Create the necessary continuation variables and add the function
    u = Var("u"; initial_u=vec(u0))
    p = Var("p"; initial_u=p0)
    t = Var("t"; initial_u=t0)
    coll = Data("coll", (n_mesh=size(u0, 2), D=D, Du=Du))
    func = Func(
        "f", FourierColl(f!, size(u0, 1)), (u, p, t), (coll,); initial_dim=length(u0)
    )
    push!(problem, func)
    # Continuation parameters
    append!(problem, parameters(_pnames, p; top_level=top_level))
    push!(problem, parameter("t0", t; index=1, active=!fix_t0, top_level=false))
    push!(problem, parameter("t1", t; index=2, active=!fix_t1, top_level=false))
    return problem
end

"""
    $(SIGNATURES)

Create a Fourier differentiation matrix with numerical type T on the domain
`x = range(0, 2π, length=N+1)[1:end-1]`.
"""
function fourier_diff(T::Type{<:Number}, N::Integer; order=1)
    D = zeros(T, N, N)
    n1 = (N - 1) ÷ 2
    n2 = N ÷ 2
    x = LinRange{T}(0, π, N + 1)
    if order == 1
        for i in 2:N
            sgn = (one(T) / 2 - iseven(i))
            D[i, 1] = iseven(N) ? sgn * cot(x[i]) : sgn * csc(x[i])
        end
    elseif order == 2
        D[1, 1] =
            iseven(N) ? -N^2 * one(T) / 12 - one(T) / 6 : -N^2 * one(T) / 12 + one(T) / 12
        for i in 2:N
            sgn = -(one(T) / 2 - iseven(i))
            D[i, 1] = iseven(N) ? sgn * csc(x[i]) .^ 2 : sgn * cot(x[i]) * csc(x[i])
        end
    else
        error("Not implemented")
    end
    for j in 2:N
        D[1, j] = D[N, j - 1]
        D[2:N, j] .= D[1:(N - 1), j - 1]
    end
    return D
end
fourier_diff(N::Integer; kwargs...) = fourier_diff(Float64, N; kwargs...)

end  # module


text:
using Statistics, PaddedMatrices
using AsymptoticPosteriors, Parameters, Rmath
using Base.Cartesian, Random, LinearAlgebra

@inline beta_lpdf(x, alpha, beta) = (alpha-1)*log(x) + (beta-1)*log(1-x)

"""Log Density of the No Gold-Standard Model"""
function NoGoldStandard(priordata, S₁::T, S₂::T, C₁::T, C₂::T, π::Vararg{T,K}) where {T,K}
    @unpack n, απ, βπ, αS1, βS1, αS2, βS2,
        αC1, βC1, αC2, βC2 = priordata
    target = beta_lpdf(S₁, αS1, βS1) + beta_lpdf(S₂, αS2, βS2)
    target += beta_lpdf(C₁, αC1, βC1) + beta_lpdf(C₂, αC2, βC2)
    omS₁, omS₂, omC₁, omC₂ = 1-S₁, 1-S₂, 1-C₁, 1-C₂
    πₛ₁, πc₁ = S₁*S₂, C₁*C₂
    πₛ₂, πₛ₃ = S₁ - πₛ₁, S₂ - πₛ₁
    πc₂, πc₃ = C₁ - πc₁, C₂ - πc₁
    πₛ₄, πc₄ = omS₁ - πₛ₃, omC₁ - πc₃
    @inbounds for k in 1:K
        πₖ = π[k]
        omπₖ = 1 - πₖ
        target += beta_lpdf( πₖ, απ, βπ )

        target += n[1,k] * log( πₖ*πₛ₁ + omπₖ*πc₄ )
        target += n[2,k] * log( πₖ*πₛ₂ + omπₖ*πc₃ )
        target += n[3,k] * log( πₖ*πₛ₃ + omπₖ*πc₂ )
        target += n[4,k] * log( πₖ*πₛ₄ + omπₖ*πc₁ )

        target += n[5,k] * log( πₖ*  S₁ + omπₖ*omC₁ )
        target += n[6,k] * log( πₖ*omS₁ + omπₖ*  C₁ )

        target += n[7,k] * log( πₖ*  S₂ + omπₖ*omC₂ )
        target += n[8,k] * log( πₖ*omS₂ + omπₖ*  C₂ )
    end
    target
end

"""Log Density of the No Gold-Standard Model, with correlations on log-odds scale"""
function NoGoldStandardCorr(priordata, S₁::T, S₂::T, C₁::T, C₂::T, ρₛ::T, ρc::T, π::Vararg{T,K}) where {T,K}
    @unpack n, απ, βπ, αS1, βS1, αS2, βS2, αC1, βC1, αC2, βC2 = priordata
    target = beta_lpdf(S₁, αS1, βS1) + beta_lpdf(S₂, αS2, βS2) + beta_lpdf(C₁, αC1, βC1) + beta_lpdf(C₂, αC2, βC2) - 4ρₛ - 4ρc#exponential priors.

    omS₁, omS₂, omC₁, omC₂ = 1-S₁, 1-S₂, 1-C₁, 1-C₂
    S₁S₂ = S₁*S₂
    Sₐ = S₁ + S₂ - S₁S₂
    Sₒ = Sₐ - S₁S₂
    expnρₛ = exp(-ρₛ)
    Sdenom = 1 / (1 - Sₒ + Sₒ * expnρₛ)
    πₛ₁ = S₁S₂ * Sdenom
    πₛ₂ = (S₁-S₁S₂) * expnρₛ * Sdenom
    πₛ₃ = (S₂-S₁S₂) * expnρₛ * Sdenom
    πₛ₄ = (1 - Sₐ) * Sdenom
    C₁C₂ = C₁*C₂
    Cₐ = C₁ + C₂ - C₁C₂
    Cₒ = Cₐ - C₁C₂
    expnρc = exp(-ρc)
    Cdenom = 1 / (1 - Cₒ + Cₒ * expnρc)
    πc₁ = C₁C₂ * Cdenom
    πc₂ = (C₁-C₁C₂) * expnρc * Cdenom
    πc₃ = (C₂-C₁C₂) * expnρc * Cdenom
    πc₄ = (1 - Cₐ) * Cdenom

    @inbounds for k in 1:K
        πₖ = π[k]
        omπₖ = 1 - πₖ
        target += beta_lpdf( πₖ, απ, βπ )

        target += n[1,k] * log( πₖ*πₛ₁ + omπₖ*πc₄ )
        target += n[2,k] * log( πₖ*πₛ₂ + omπₖ*πc₃ )
        target += n[3,k] * log( πₖ*πₛ₃ + omπₖ*πc₂ )
        target += n[4,k] * log( πₖ*πₛ₄ + omπₖ*πc₁ )

        target += n[5,k] * log( πₖ*  S₁ + omπₖ*omC₁ )
        target += n[6,k] * log( πₖ*omS₁ + omπₖ*  C₁ )

        target += n[7,k] * log( πₖ*  S₂ + omπₖ*omC₂ )
        target += n[8,k] * log( πₖ*omS₂ + omπₖ*  C₂ )
    end

    target
end

@generated function NGS(data, params::AbstractVector{T}, ::Val{N}) where {N,T}
    call = generate_call_transform(:NoGoldStandard, N )
#    Nm4 = N - 4
    quote
        @nexprs $N i -> begin
            @inbounds p_i = inv_logit(params[i])
        end
        $call
        target
    end
end

@inline logit(x) = log(x/(1-x))
@inline inv_logit(x) = 1/(1+exp(-x))

struct Multinomial{N,T}
    n::Int
    p::NTuple{N,T}
    cumulative_p::NTuple{N,T}
end
function Multinomial(n, p::NTuple{N,T}) where {N,T}
    total = Ref{T}(zero(T)) #Cumsum not defined for ntuples; this gets around closure bug.
    cumulative_p = ntuple( i -> begin
        total[] += p[i]
    end, Val{N}())
    Multinomial(n, p, cumulative_p)
end
function Random.rand(m::Multinomial{N,T}, n::Int...) where {N,T}
    rand!(Array{NTuple{N,Int}}(undef, n...), m)
end
@generated function Random.rand!(out::AbstractArray{NTuple{N,I}}, m::Multinomial{N,T}, rng = Compat.Random.GLOBAL_RNG) where {N,T,I<:Integer}
    quote
        for i ∈ eachindex(out)
            @nexprs $N d -> x_d = zero(I)
            for j ∈ 1:m.n
                r = rand(rng)
                @nif $N d -> r < m.cumulative_p[d] d -> x_d += one(I)
            end
            out[i] = rand(m)
        end
        out
    end
end
@generated function Random.rand(m::Multinomial{N,T}, rng = Compat.Random.GLOBAL_RNG) where {N,T}
    quote
        @nexprs $N d -> x_d = 0
        for j ∈ 1:m.n
            r = rand(rng)
            @nif $N d -> r < m.cumulative_p[d] d -> x_d += 1
        end
        @ntuple $N d -> x_d
    end
end
# struct CorrErrors{p, T<: Real}
#     π::SVector{p,T}
#     S::SVector{2,T}
#     C::SVector{2,T}
#     covsc::SVector{2,T}
# end
# function CorrErrors(π::NTuple{p,T}, S::Tuple{T,T}, C::Tuple{T,T}, covsc::Tuple{T,T} = (0.0,0.0)) where {p,T}
#     CorrErrors{p,T}(SVector(π), SVector(S), SVector(C), SVector(covsc))
# end
mutable struct CorrErrors{P,T<: Real,PP}
    π::NTuple{P,T}
    S::NTuple{2,T}
    C::NTuple{2,T}
    covsc::NTuple{2,T}
    p_stick::NTuple{P,NTuple{3,T}}
    p_stick_independent::NTuple{P,NTuple{3,T}}
    p_marginal1::NTuple{P,T}
    p_marginal2::NTuple{P,T}
    Θ::NTuple{PP,T}
end
function CorrErrors(π::NTuple{P,T}, S::Tuple{T,T}, C::Tuple{T,T}, covsc::Tuple{T,T} = (0.0,0.0)) where {P,T}
    p = common_p(π, S, C, covsc)
    p_ind = common_p(π, S, C, (0.0,0.0))
    CorrErrors(π, S, C, covsc, break_sticks.(p), break_sticks.(p_ind), p_i(π, S, C, 1), p_i(π, S, C, 2), unconstrain(π, S, C, covsc))
end
#function Random.rand!(t::CoreErrors{P,T,PP})
#
#end
# function CorrErrors
@generated function unconstrain(π::NTuple{P}, S::NTuple{2}, C::NTuple{2}, covsc::NTuple{2}) where P
    PP = P+6
    quote
        x_1 = logit(2S[1]-1)
        x_2 = logit(2S[2]-1)
        x_3 = logit(2C[1]-1)
        x_4 = logit(2C[2]-1)
        x_5 = log(covsc[1])
        x_6 = log(covsc[2])
        @nexprs $P i -> begin
            x_{i+6} = logit(π[i])
        end
        @ntuple $PP i -> x_i
    end
end
function common_p(Θ::CorrErrors{P,Float64}) where P
    @unpack π, S, C, covsc = Θ
    common_p(π, S, C, covsc)
end
function common_p(π::NTuple{P}, S::NTuple{2}, C::NTuple{2}, covsc::NTuple{2}) where P
    omπ = 1 .- π
    S₁, S₂ = S
    C₁, C₂ = C
    ρₛ, ρc = covsc
    omS₁, omS₂, omC₁, omC₂ = 1-S₁, 1-S₂, 1-C₁, 1-C₂
    S₁S₂ = S₁*S₂
    Sₐ = S₁ + S₂ - S₁S₂
    Sₒ = Sₐ - S₁S₂
    expnρₛ = exp(-ρₛ)
    Sdenom = 1 / (1 - Sₒ + Sₒ * expnρₛ)
    πₛ₁ = S₁S₂ * Sdenom
    πₛ₂ = (S₁-S₁S₂) * expnρₛ * Sdenom
    πₛ₃ = (S₂-S₁S₂) * expnρₛ * Sdenom
    πₛ₄ = (1 - Sₐ) * Sdenom
    C₁C₂ = C₁*C₂
    Cₐ = C₁ + C₂ - C₁C₂
    Cₒ = Cₐ - C₁C₂
    expnρc = exp(-ρc)
    Cdenom = 1 / (1 - Cₒ + Cₒ * expnρc)
    πc₁ = C₁C₂ * Cdenom
    πc₂ = (C₁-C₁C₂) * expnρc * Cdenom
    πc₃ = (C₂-C₁C₂) * expnρc * Cdenom
    πc₄ = (1 - Cₐ) * Cdenom
  ntuple( i -> (
      π[i]*πₛ₁ + omπ[i]*πc₄,
      π[i]*πₛ₂ + omπ[i]*πc₃,
      π[i]*πₛ₃ + omπ[i]*πc₂,
      π[i]*πₛ₄ + omπ[i]*πc₁
  ), Val{P}())
end
function p_i(Θ::CorrErrors, i::Int)
    @unpack π, S, C = Θ
    p_i(π, S, C, i)
end
function p_i(π::NTuple{P}, S::NTuple{2}, C::NTuple{2}, i) where P
    ntuple(j -> π[j]*S[i] + (1-π[j])*(1-C[i]), Val{P}())
end
function p_i2(Θ::CorrErrors, i::Int)
    @unpack π, S, C = Θ
    p_i2(π, S, C, i)
end
function p_i2(π::NTuple{P}, S::NTuple{2}, C::NTuple{2}, i) where P
  [(π[j]*S[i] + (1-π[j])*(1-C[i]), π[j]*(1-S[i]) + (1-π[j])*C[i]) for j ∈ 1:P]
end
@generated ValM1(::Val{N}) where N = Val{N-1}()
@generated function break_sticks(tup::NTuple{N,T}) where {N,T}
    Nm2 = N-2
    Nm1 = N-1
    quote
        p_0 = tup[1]
        cumulative_complete = one(T)
        @nexprs $Nm2 d -> begin
            cumulative_complete -= tup[d]
            p_d = tup[d+1] / cumulative_complete
        end
        @ntuple $Nm1 d -> p_{d-1}
    end
end

# gen_data(Θ::CorrErrors{P}, n_common, n_1_only, n_2_only) where P = gen_data!(Matrix{Int}(undef, 8, P), Θ, n_common, n_1_only, n_2_only)
# function gen_data!(out::Matrix{Int}, Θ::CorrErrors{P,T}, n_common::Int, n_1_only::Int, n_2_only::Int) where {P,T}
#     double_test = common_p(Θ)
#     p_1_only = p_i2(Θ, 1)
#     p_2_only = p_i2(Θ, 2)
#     for i ∈ 1:P
#         out[1:4,i] .= rand(Multinomial(n_common, double_test[i]))
#         out[5:6,i] .= rand(Multinomial(n_1_only, p_1_only[i]))
#         out[7:8,i] .= rand(Multinomial(n_2_only, p_2_only[i]))
#     end
#     out
# end

gen_data(Θ::CorrErrors{P}, n_common, n_1_only, n_2_only) where {P} = gen_data!(Matrix{Int}(undef, 8, P), Θ, n_common, n_1_only, n_2_only)
function gen_data!(out::Matrix{Int}, Θ::CorrErrors{P,T}, n_common::Int, n_1_only::Int, n_2_only::Int) where {P,T}
    for i ∈ 1:P
        cum_p = zero(T)
        x = rbinom(1, n_common, Θ.p_stick[i][1])[1]
        out[1,i] = x
        n = n_common - x
        x = rbinom(1, n, Θ.p_stick[i][2])[1]
        out[2,i] = x
        n -= x
        x = rbinom(1, n, Θ.p_stick[i][3])[1]
        out[3,i] = x
        out[4,i] = n - x
        x = rbinom(1, n_1_only, Θ.p_marginal1[i])[1]
        out[5,i] = x
        out[6,i] = n_1_only - x
        x = rbinom(1, n_2_only, Θ.p_marginal2[i])[1]
        out[7,i] = x
        out[8,i] = n_2_only - x
    end
    out
end
function gen_data!(out::Matrix{Int}, Θ::CorrErrors{P,T}, n_common::NTuple{P,Int}, n_1_only::NTuple{P,Int}, n_2_only::NTuple{P,Int}) where {P,T}
    for p ∈ 1:P
        cum_p = zero(T)
        x = rbinom(1, n_common[p], Θ.p_stick[p][1])[1]
        out[1,p] = x
        n = n_common[p] - x
        x = rbinom(1, n, Θ.p_stick[p][2])[1]
        out[2,p] = x
        n -= x
        x = rbinom(1, n, Θ.p_stick[p][3])[1]
        out[3,p] = x
        out[4,p] = n - x
        x = rbinom(1, n_1_only[p], Θ.p_marginal1[p])[1]
        out[5,p] = x
        out[6,p] = n_1_only[p] - x
        x = rbinom(1, n_2_only[p], Θ.p_marginal2[p])[1]
        out[7,p] = x
        out[8,p] = n_2_only[p] - x
    end
    out
end
# corr_errors_independent2 = CorrErrors((0.1,        0.4), (0.9, 0.95), (0.85,0.97));
# corr_errors_independent3 = CorrErrors((0.1,  0.25, 0.4), (0.9, 0.95), (0.85,0.97));
# corr_errors_independent4 = CorrErrors((0.1,0.2,0.3,0.4), (0.9, 0.95), (0.85,0.97));

# init6 = logit.([0.9*2-1, 0.95*2-1, 0.85*2-1, 0.97*2-1, 0.1, 0.4]);
# init7 = logit.([0.9*2-1, 0.95*2-1, 0.85*2-1, 0.97*2-1, 0.1, 0.25, 0.4]);
# init8 = logit.([0.9*2-1, 0.95*2-1, 0.85*2-1, 0.97*2-1, 0.1, 0.2, 0.3, 0.4]);
# initCorr8 = logit.([0.9*2-1, 0.95*2-1, 0.85*2-1, 0.97*2-1, 0.5/1.5, 0.45/1.45, 0.1, 0.4]);
# initCorr9 = logit.([0.9*2-1, 0.95*2-1, 0.85*2-1, 0.97*2-1, 0.5/1.5, 0.45/1.45, 0.1, 0.25, 0.4]);
# initCorr10 = logit.([0.9*2-1, 0.95*2-1, 0.85*2-1, 0.97*2-1, 0.5/1.5, 0.45/1.45, 0.1, 0.2, 0.3, 0.4]);

# @with_kw struct NoGoldData
#     απ::Float64 = 1.0
#     βπ::Float64 = 1.0
#     αS1::Float64 = 1.0
#     βS1::Float64 = 1.0
#     αS2::Float64 = 1.0
#     βS2,::Float64 = 1.0
#     αC1::Float64 = 1.0
#     βC1::Float64 = 1.0
#     αC2::Float64 = 1.0
#     βC2::Float64 = 1.0
#     n::Matrix{Int}
# end

abstract type NGD{N} end
@with_kw struct NoGoldData{N} <: NGD{N}
    n::Matrix{Int}
    απ::Float64 = 1.0 # απ
    βπ::Float64 = 1.0 # βπ
    αS1::Float64 = 1.0 # αS1
    βS1::Float64 = 1.0 # βS1
    αS2::Float64 = 1.0 # αS2
    βS2::Float64 = 1.0 # βS2
    αC1::Float64 = 1.0 # αC1
    βC1::Float64 = 1.0 # βC1
    αC2::Float64 = 1.0 # αC2
    βC2::Float64 = 1.0 # βC2
end
@with_kw struct NoGoldDataCorr{N} <: NGD{N}
    n::Matrix{Float64}
    απ::Float64 = 1.0 # απ
    βπ::Float64 = 1.0 # βπ
    αS1::Float64 = 38.0 # αS1
    βS1::Float64 = 2.0 # βS1
    αS2::Float64 = 38.0 # αS2
    βS2::Float64 = 2.0 # βS2
    αC1::Float64 = 38.0 # αC1
    βC1::Float64 = 2.0 # βC1
    αC2::Float64 = 38.0 # αC2
    βC2::Float64 = 2.0 # βC2
end

function generate_call_transform(fname, n)
    # verind = VERSION > v"0.6.9" ? 3 : 2
    out = quote
        target = $fname( data, (p_1+1)/2 )
        target += log(p_1) + log(1 - p_1)
        target += log(p_5) + log(1 - p_5)
        target
    end
    for i ∈ 2:4
        p_i = Symbol(:p_, i)
        push!(out.args[2].args[2].args, :( ($p_i+1)/2 ))
        push!(out.args[4].args[2].args, :( log($p_i) ))
        push!(out.args[4].args[2].args, :( log(1-$p_i)  ))
    end
    push!(out.args[2].args[2].args, :(p_5))
    for i ∈ 6:n
        p_i = Symbol(:p_, i)
        push!(out.args[2].args[2].args, :($p_i))
        push!(out.args[6].args[2].args, :( log($p_i) ))
        push!(out.args[6].args[2].args, :( log(1-$p_i)  ))
    end
    out
end

@generated function NGS(data::NoGoldData{N}, params::AbstractVector{T}) where {N,T}
    call = generate_call_transform(:NoGoldStandard, N )
#    Nm4 = N - 4
    quote
        @nexprs $N i -> begin
            @inbounds p_i = inv_logit(params[i])
        end
        $call
        target
    end
end

function generate_call_transformc(fname, n)
    out = quote
        target = 2p_5 + 2p_6 + log(p_1) + log(1 - p_1)  #2p_5 and 2p_6 because Gamma(2,4) prior + Jacobian (+p_5+p_6).
        p_1, p_2, p_3, p_4 = (p_1+1)/2, (p_2+1)/2, (p_3+1)/2, (p_4+1)/2
        target += $fname( data, p_1 )
        target += log(p_7) + log(1 - p_7)
        #target += log(p_5) + log(1 - p_5) + log(p_6) + log(1 - p_6)
        target
    end
    for i ∈ 2:4
        p_i = Symbol(:p_, i)
        push!(out.args[6].args[2].args, :( $p_i ))
        push!(out.args[2].args[2].args, :( log($p_i)  ))
        push!(out.args[2].args[2].args, :( log(1-$p_i)  ))
    end
    push!(out.args[6].args[2].args, :(exp(p_5)))
    push!(out.args[6].args[2].args, :(exp(p_6)))
    push!(out.args[6].args[2].args, :p_7)
    for i ∈ 8:n
        p_i = Symbol(:p_, i)
        push!(out.args[6].args[2].args, :($p_i))
        push!(out.args[8].args[2].args, :( log($p_i) ))
        push!(out.args[8].args[2].args, :( log(1-$p_i)  ))
    end
    out
end
@generated function NGSC(data::NoGoldDataCorr{N}, params::AbstractVector{T}) where {N,T}
    call = generate_call_transformc(:NoGoldStandardCorr, N )
    Nm6 = N - 6
    quote
        @inbounds  begin
            @nexprs 4 i -> begin
                p_i = inv_logit(params[i])
            end
            p_5 = params[5]
            p_6 = params[6]
            @nexprs $Nm6 i -> begin
                p_{i+6} = inv_logit(params[i+6])
            end
        end
        $call
        target
    end
end
(data::NoGoldData{N})(x) where N = NGS(data, x)
(data::NoGoldDataCorr{N})(x) where N = NGSC(data, x)

@generated ValP6(::Val{N}) where N = Val{N+6}()
@generated ValP4(::Val{N}) where N = Val{N+4}()
@generated function NoGoldData(Θ::CorrErrors{P,T}, n_common::Int, n_1_only::Int, n_2_only::Int; kwargs...) where {P,T}
    :(NoGoldData{$(P+4)}(n = gen_data(Θ,  n_common,  n_1_only,  n_2_only), kwargs...))
end
@generated function NoGoldData(Θ::CorrErrors{P,T}; kwargs...) where {P,T}
    :(NoGoldData{$(P+4)}(n = Matrix{Int}(undef,8,$P), kwargs...))
end
function NoGoldData(Θ::CorrErrors{P}, n_common::Int, n_1_only::Int, n_2_only::Int; kwargs...) where P
    NoGoldData(ValP4(Val{P}()); n = gen_data(Θ,  n_common,  n_1_only,  n_2_only), kwargs...)
end
function NoGoldData(::Val{N}; kwargs...) where {N}
    NoGoldData{N}(; kwargs...)
end
function NoGoldDataCorr(Θ::CorrErrors{P}, n_common::Int, n_1_only::Int, n_2_only::Int; kwargs...) where P
    NoGoldDataCorr(ValP6(Val{P}()); n = gen_data(Θ,  n_common,  n_1_only,  n_2_only), kwargs...)
end
function NoGoldDataCorr(::Val{N}; kwargs...) where {N}
    NoGoldDataCorr{N}(; kwargs...)
end
function resample!(data::NGD{N}, Θ, n_common, n_1_only, n_2_only) where N
    gen_data!(data.n, Θ, n_common, n_1_only, n_2_only)
end
#function resample!(data::NGDC{N}, Θ, n_common, n_1_only, n_2_only) where N
#    gen_data!(data.n, Θ, n_common, n_1_only, n_2_only)
#end

# n_small_2  = NoGoldData(corr_errors_independent2,  100,  200,  16);
# n_medium_2 = NoGoldData(corr_errors_independent2,  400,  800,  64);
# n_big_2    = NoGoldData(corr_errors_independent2, 1600, 3200, 252);

# n_small_3  = NoGoldData(corr_errors_independent3,  100,  200,  16);
# n_medium_3 = NoGoldData(corr_errors_independent3,  400,  800,  64);
# n_big_3    = NoGoldData(corr_errors_independent3, 1600, 3200, 252);

# n_small_4  = NoGoldData(corr_errors_independent4,  100,  200,  16);
# n_medium_4 = NoGoldData(corr_errors_independent4,  400,  800,  64);
# n_big_4    = NoGoldData(corr_errors_independent4, 1600, 3200, 252);

# zeros6 = fill(0.0,6);
# zeros7 = fill(0.0,7);
# zeros8 = fill(0.0,8);
# zeros9 = fill(0.0,9);
# zeros10 = fill(0.0,10);

function inverse_transform(x)
    out = similar(x)
    @inbounds for i ∈ 1:4
        out[i] = (inv_logit(x[i]) + 1) / 2
    end
    @inbounds for i ∈ 5:length(x)
        out[i] = inv_logit(x[i])
    end
    out
end
function inverse_transformc(x)
    out = similar(x)
    @inbounds for i ∈ 1:4
        out[i] = (inv_logit(x[i]) + 1) / 2
    end
    @inbounds for i ∈ 5:6
        out[i] = exp(x[i])
    end
    @inbounds for i ∈ 7:length(x)
        out[i] = inv_logit(x[i])
    end
    out
end

interval(ap, alpha, ind) = (quantile(ap,alpha/2,ind), quantile(ap,1-alpha/2,ind))
four_num_sum(ap, ind) = (quantile(ap,0.025,ind), quantile(ap,0.25,ind), quantile(ap,0.75,ind), quantile(ap,0.975,ind))
function Base.summary(ap::AsymptoticPosteriorFD{N,T}; sigfigs=4) where {N,T}
    sig(x, s) = round(x, sigdigits = s)
    println("S1:\t\t\t\t", sig.((inv_logit.(four_num_sum(ap, 1)) .+1) ./ 2, sigfigs))
    println("S2:\t\t\t\t", sig.((inv_logit.(four_num_sum(ap, 2)) .+1) ./ 2, sigfigs))
    println("C1:\t\t\t\t", sig.((inv_logit.(four_num_sum(ap, 3)) .+1) ./ 2, sigfigs))
    println("C2:\t\t\t\t", sig.((inv_logit.(four_num_sum(ap, 4)) .+1) ./ 2, sigfigs))
    for i in 1:N-4
        println("Population $i proportion:\t", sig.(inv_logit.(four_num_sum(ap, i+4)), sigfigs))
    end
end
function corrsummary(ap::AsymptoticPosteriorFD{N,T}; sigfigs=4) where {N,T}
    sig(x, s) = round(x, sigdigits = s)
    println("S1:\t\t\t\t", sig.((inv_logit.(four_num_sum(ap, 1)) .+1) ./ 2, sigfigs))
    println("S2:\t\t\t\t", sig.((inv_logit.(four_num_sum(ap, 2)) .+1) ./ 2, sigfigs))
    println("C1:\t\t\t\t", sig.((inv_logit.(four_num_sum(ap, 3)) .+1) ./ 2, sigfigs))
    println("C2:\t\t\t\t", sig.((inv_logit.(four_num_sum(ap, 4)) .+1) ./ 2, sigfigs))
    println("Cov1:\t\t\t\t", sig.(exp.(four_num_sum(ap, 5)), sigfigs))
    println("Cov2:\t\t\t\t", sig.(exp.(four_num_sum(ap, 6)), sigfigs))
    for i in 1:N-6
        println("Population $i proportion:\t", sig.(inv_logit.(four_num_sum(ap, i+6)), sigfigs))
    end
end
function corrsummarytup(ap::AsymptoticPosteriorFD{8,T}) where {T}
    (inv_logit.(four_num_sum(ap, 1)) .+1) ./ 2, (inv_logit.(four_num_sum(ap, 2)) .+1) ./ 2,
    (inv_logit.(four_num_sum(ap, 3)) .+1) ./ 2, (inv_logit.(four_num_sum(ap, 4)) .+1) ./ 2,
    exp.(four_num_sum(ap, 5)), exp.(four_num_sum(ap, 6)),
    inv_logit.(four_num_sum(ap, 7)), inv_logit.(four_num_sum(ap, 8))
end

# corr_errors_dependent2 = CorrErrors((0.1,        0.4), (0.9, 0.95), (0.85,0.97), (0.5,0.45));
# corr_errors_dependent3 = CorrErrors((0.1,  0.25, 0.4), (0.9, 0.95), (0.85,0.97), (0.5,0.45));
# corr_errors_dependent4 = CorrErrors((0.1,0.2,0.3,0.4), (0.9, 0.95), (0.85,0.97), (0.5,0.45));

# αS1 = 36
# βS1 = 4
# αS2 = 38
# βS2 = 2
# αC1 = 32
# βC1 = 8
# αC2 = 38
# βC2 = 2

# n_small_2c  = NoGoldDataCorr(corr_errors_dependent2,100,200,16,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);
# n_medium_2c = NoGoldDataCorr(corr_errors_dependent2,400,800,64,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);
# n_big_2c    = NoGoldDataCorr(corr_errors_dependent2,1600,3200,252,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);

# n_small_3c  = NoGoldDataCorr(corr_errors_dependent3,100,200,16,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);
# n_medium_3c = NoGoldDataCorr(corr_errors_dependent3,400,800,64,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);
# n_big_3c    = NoGoldDataCorr(corr_errors_dependent3,1600,3200,252,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);

# n_small_4c  = NoGoldDataCorr(corr_errors_dependent4,100,200,16,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);
# n_medium_4c = NoGoldDataCorr(corr_errors_dependent4,400,800,64,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);
# n_big_4c    = NoGoldDataCorr(corr_errors_dependent4,1600,3200,252,αS1=αS1,βS1=βS1,αS2=αS2,βS2=βS2,αC1=αC1,βC1=βC1,αC2=αC2,βC2=βC2);


text:
"""
sanitize_precompile()

Prepares and sanitizes a julia file for precompilation. This removes any non-concrete
methods and anything non-Base or StdLib.
"""
function sanitize_precompile(precompile_lines::Vector{String})
    lines = String[]
    for line in precompile_lines
        # Generally any line with where is non-concrete, so we can skip.
        # Symbol is also runtime dependent so skip as well
        if isempty(line) || contains(line, '#') || contains(line, "where") || contains(line, "Symbol(")
            continue
        else
            try
                if can_precompile(Meta.parse(line))
                    push!(lines, line)
                end
            catch err
                show(err)
            end
        end
    end
    lines
end

"""
can_precompile

Determine if this expr is something we can precompile.
"""
function can_precompile(ex::Expr)
    if ex.head !== :call && ex.args[1] !== :precompile
        return false
    else
        return is_bakeable(ex)
    end
end

# TODO: Some Base are marked nospecialize, so we should filter these out also
"""
recurse through the call and make sure everything is in Base, Core, or a StdLib
"""
function is_bakeable(ex::Expr)

    # handle submodule (this might not be robust)
    if ex.head === :. && is_bakeable(ex.args[1])
        return true
    end

    for arg in ex.args
        #@show arg, typeof(arg)
        if is_bakeable(arg)
            continue
        else
            return false
        end
    end
    return true
end

function is_bakeable(ex::Symbol)
    for lib in bakeable_libs
        if isdefined(lib, ex)
            return true
        else
            continue
        end
    end
    return false
end

function is_bakeable(ex::QuoteNode)
    return is_bakeable(ex.value)
end

function is_bakeable(n::T) where T <: Number
    return true
end


text:
## REGRESSOR METRICS (FOR DETERMINISTIC PREDICTIONS)

struct MAE <: Measure end

"""
    mae(ŷ, y)
    mae(ŷ, y, w)

Mean absolute error.

``\\text{MAE} =  n^{-1}∑ᵢ|yᵢ-ŷᵢ|`` or ``\\text{MAE} = n^{-1}∑ᵢwᵢ|yᵢ-ŷᵢ|``

For more information, run `info(mae)`.
"""
const mae = MAE()
const mav = MAE()

metadata_measure(MAE;
    name                     = "mae",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = true,
    docstring                = "mean absolute error; aliases: `mae`, `mav`")

function (::MAE)(ŷ::Vec{<:Real}, y::Vec{<:Real})
    check_dimensions(ŷ, y)
    ret = zero(eltype(y))
    for i in eachindex(y)
        dev = abs(y[i] - ŷ[i])
        ret += dev
    end
    return ret / length(y)
end

function (::MAE)(ŷ::Vec{<:Real}, y::Vec{<:Real},
                 w::Vec{<:Real})
    check_dimensions(ŷ, y)
    check_dimensions(y, w)
    ret = zero(eltype(y))
    for i in eachindex(y)
        dev = abs(y[i] - ŷ[i])
        ret += w[i]*dev
    end
    return ret / length(y)
end

struct RMS <: Measure end
"""
    rms(ŷ, y)
    rms(ŷ, y, w)

Root mean squared error:

``\\text{RMS} = \\sqrt{n^{-1}∑ᵢ|yᵢ-ŷᵢ|^2}`` or ``\\text{RMS} = \\sqrt{\\frac{∑ᵢwᵢ|yᵢ-ŷᵢ|^2}{∑ᵢwᵢ}}``

For more information, run `info(rms)`.
"""
const rms = RMS()

metadata_measure(RMS;
    name                     = "rms",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    aggregation              = RootMeanSquare(),
    is_feature_dependent     = false,
    supports_weights         = true,
    docstring                = "root mean squared; aliases: `rms`.")

function (::RMS)(ŷ::Vec{<:Real}, y::Vec{<:Real})
    check_dimensions(ŷ, y)
    ret = zero(eltype(y))
    for i in eachindex(y)
        dev = (y[i] - ŷ[i])^2
        ret += dev
    end
    return sqrt(ret / length(y))
end

function (::RMS)(ŷ::Vec{<:Real}, y::Vec{<:Real},
                 w::Vec{<:Real})
    check_dimensions(ŷ, y)
    ret = zero(eltype(y))
    for i in eachindex(y)
        dev = (y[i] - ŷ[i])^2
        ret += w[i]*dev
    end
    return sqrt(ret / length(y))
end

struct L2 <: Measure end

"""
    l2(ŷ, y)
    l2(ŷ, y, w)

L2 per-observation loss.

For more information, run `info(l2)`.
"""
const l2 = L2()

metadata_measure(L2;
    name                     = "l2",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = true,
    is_feature_dependent     = false,
    supports_weights         = true,
    docstring                = "squared deviations; aliases: `l2`.")

function (::L2)(ŷ::Vec{<:Real}, y::Vec{<:Real})
    check_dimensions(ŷ, y)
    return (y - ŷ).^2
end

function (::L2)(ŷ::Vec{<:Real}, y::Vec{<:Real},
                w::Vec{<:Real})
    check_dimensions(ŷ, y)
    check_dimensions(w, y)
    return w .* (y - ŷ).^2
end

struct L1 <: Measure end

"""
    l1(ŷ, y)
    l1(ŷ, y, w)

L1 per-observation loss.

For more information, run `info(l1)`.
"""
const l1 = L1()

metadata_measure(L1;
    name                     = "l1",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = true,
    is_feature_dependent     = false,
    supports_weights         = true,
    docstring                = "absolute deviations; aliases: `l1`.")

function (::L1)(ŷ::Vec{<:Real}, y::Vec{<:Real})
    check_dimensions(ŷ, y)
    return abs.(y - ŷ)
end

function (::L1)(ŷ::Vec{<:Real}, y::Vec{<:Real},
                w::Vec{<:Real})
    check_dimensions(ŷ, y)
    check_dimensions(w, y)
    return w .* abs.(y - ŷ)
end

struct RMSL <: Measure end

"""
    rmsl(ŷ, y)

Root mean squared logarithmic error:

``\\text{RMSL} = n^{-1}∑ᵢ\\log\\left({yᵢ \\over ŷᵢ}\\right)``

For more information, run `info(rmsl)`.

See also [`rmslp1`](@ref).
"""
const rmsl = RMSL()

metadata_measure(RMSL;
    name                     = "rmsl",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    aggregation              = RootMeanSquare(),
    is_feature_dependent     = false,
    supports_weights         = false,
    docstring                = "root mean square logarithm; aliases: `rmsl`.")

function (::RMSL)(ŷ::Vec{<:Real}, y::Vec{<:Real})
    check_dimensions(ŷ, y)
    ret = zero(eltype(y))
    for i in eachindex(y)
        dev = (log(y[i]) - log(ŷ[i]))^2
        ret += dev
    end
    return sqrt(ret / length(y))
end

struct RMSLP1 <: Measure end
"""
    rmslp1(ŷ, y)

Root mean squared logarithmic error with an offset of 1:

``\\text{RMSLP1} = n^{-1}∑ᵢ\\log\\left({yᵢ + 1 \\over ŷᵢ + 1}\\right)``

For more information, run `info(rmslp1)`.

See also [`rmsl`](@ref).
"""
const rmslp1 = RMSLP1()

metadata_measure(RMSLP1;
    name                     = "rmslp1",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    aggregation              = RootMeanSquare(),
    is_feature_dependent     = false,
    supports_weights         = false,
    docstring                = "root mean squared logarithm plus one; " *
                               "aliases: `rmslp1`.")

function (::RMSLP1)(ŷ::Vec{<:Real}, y::Vec{<:Real})
    check_dimensions(ŷ, y)
    ret = zero(eltype(y))
    for i in eachindex(y)
        dev = (log(y[i] + 1) - log(ŷ[i] + 1))^2
        ret += dev
    end
    return sqrt(ret / length(y))
end

struct RMSP <: Measure end

"""
    rmsp(ŷ, y)

Root mean squared proportional loss:

``\\text{RMSP} = m^{-1}∑ᵢ \\left({yᵢ-ŷᵢ \\over yᵢ}\\right)^2``

where the sum is over indices such that `yᵢ≂̸0` and `m` is the number
of such indices.

For more information, run `info(rmsp)`.
"""
const rmsp = RMSP()

metadata_measure(RMSP;
    name                     = "rmsp",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    aggregation              = RootMeanSquare(),
    is_feature_dependent     = false,
    supports_weights         = false,
    docstring                = "root mean square proportions; aliases: `rmsp`.")

function (::RMSP)(ŷ::Vec{<:Real}, y::Vec{<:Real}, tol=eps())
    check_dimensions(ŷ, y)
    ret = zero(eltype(y))
    count = 0
    @inbounds for i in eachindex(y)
        ayi = abs(y[i])
        if ayi > tol
            dev = ((y[i] - ŷ[i]) / ayi)^2
            ret += dev
            count += 1
        end
    end
    return sqrt(ret / count)
end

struct MAPE <: Measure
    tol::Real
end

MAPE(; tol=eps()) = MAPE(tol)

"""
     MAPE(; tol=esp())

Mean Absolute Proportional Error:

``\\text{MAPE} =  m^{-1}∑ᵢ|{(yᵢ-ŷᵢ) \\over yᵢ}|``
where the sum is over indices such that `yᵢ > tol` and `m` is the number
of such indices.

For more information, run `info(mape)`.
"""
const mape = MAPE()

metadata_measure(MAPE;
    name                     = "mape",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = false,
    docstring                = "Mean Absolute Proportional Error; "*
                 "aliases: `mape=MAPE()`.")

function (m::MAPE)(ŷ::Vec{<:Real}, y::Vec{<:Real})
    check_dimensions(ŷ, y)
    ret = zero(eltype(y))
    count = 0
    @inbounds for i in eachindex(y)
        ayi = abs(y[i])
        if ayi > m.tol
        #if y[i] != zero(eltype(y))
            dev = abs((y[i] - ŷ[i]) / ayi)
            #dev = abs((y[i] - ŷ[i]) / y[i])
            ret += dev
            count += 1
        end
    end
    return ret / count
end

struct LogCosh <: Measure end

"""
    log_cosh(ŷ, y)

Log-Cosh per-observation loss:

``\\text{Log-Cosh Loss} = log(cosh(ŷᵢ-yᵢ))``

where `yᵢ` is the target and `ŷᵢ` is the output.

For more information, run `info(log_cosh)`.
"""
const log_cosh = LogCosh()

metadata_measure(LogCosh;
    name                     = "log_cosh",
    target_scitype           = Union{Vec{Continuous},Vec{Count}},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = true,
    is_feature_dependent     = false,
    supports_weights         = false,
    docstring                = "log cosh loss; aliases: `log_cosh`.")

_softplus(x::T) where T<:Real = x > zero(T) ? x + log1p(exp(-x)) : log1p(exp(x))
_log_cosh(x::T) where T<:Real = x + _softplus(-2x) - log(convert(T, 2))

function (log_cosh::LogCosh)(ŷ::Vec{<:T}, y::Vec{<:T}) where T <:Real
    check_dimensions(ŷ, y)
    return _log_cosh.(ŷ - y)
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using StatsBase
using Plots
using DataFrames
using Underscores
using CSV
using HTTP, JSON
using Dates
## ============================================================================
#                    apple
#  ============================================================================
old2new_names = Dict(
    "Capital Region of Denmark" => "capital",
    "Central Denmark Region" => "central",
    "North Denmark Region" => "north",
    "Region Zealand" => "zealand",
    "Region of Southern Denmark" => "south"
)
regions = values(old2new_names)

function replacemissing(xs)
    for i in eachindex(xs)
        if ismissing(xs[i])
            xs[i] = @_ vcat( xs[i-2:i-1], xs[i+1:i+2] ) |>
                skipmissing |>
                mean
        end
    end
    convert(Vector{Float64}, xs)
end

function custom_mean(xs)
    n = 0
    s = 0.
    for x in xs
        if !ismissing(x)
            n += 1
            s += x
        end
    end
    if n == 0
        @warn "all missing"
        return missing
    else
        m = s / n
        @assert isfinite(m)
        return m
    end
end

apple = @_ projectdir("data/apple/applemobilitytrends.csv") |>
    CSV.read(__, DataFrame) |>
    filter(:country => x->!ismissing(x) , __) |>
    # filter(:geo_type => ==("country/region"), __) |>
    filter(:country => ==("Denmark") , __) |>
    filter(:geo_type => ==("sub-region") , __) |>
    # filter(:region => x-> x ∈ keys(old2new_names) , __) |>
    transform(__, :region => ByRow(x -> old2new_names[x]) => :region) |>
    select(__, Not(["geo_type", "alternative_name", "sub-region", "country"])) |>
    DataFrames.stack(__, Not([:region, :transportation_type]), variable_name="date", value_name="mobility") |>
    unstack(__, :transportation_type, :mobility) |>
    select(__, :date => ByRow(Date) => :date, Not(:date), AsTable([:driving,:walking,:transit]) => ByRow(x->custom_mean(x)) => :apple) |>
    rename(__, [:driving,:walking,:transit] .=> [:apple_driving,:apple_walking,:apple_transit])

## ============================================================================
#                    google
#  ============================================================================
const mobility_names = [
    "retail_and_recreation_percent_change_from_baseline",
    "grocery_and_pharmacy_percent_change_from_baseline",
    "workplaces_percent_change_from_baseline",
    "residential_percent_change_from_baseline"
]

const new_mobility_names = [
    "google_retail",
    "google_grocery",
    "google_workplaces",
    "google_residential"
]

google = @_ "/home/and/data/google-mobility/Global_Mobility_Report.csv" |>
    CSV.File(__) |>
    DataFrame |>
    filter(:country_region => ==("Denmark"), __) |>
    filter( :sub_region_1 => x->!ismissing(x), __) |>
    transform(__, :sub_region_1 => disallowmissing => :region) |>
    filter(:sub_region_1 => x-> x ∈ keys(old2new_names) , __) |>
    filter(:sub_region_2 => ismissing, __) |>
    select(__, :date, :region, mobility_names .=> new_mobility_names) |>
    select(__, All(), :region => ByRow(x -> old2new_names[x]) => :region) |>
    select(__, All(), AsTable([:google_retail,:google_grocery,:google_workplaces]) => ByRow(mean) => :google)

## ============================================================================
#                    telco
#  ============================================================================
eng2dan = Dict(
        "Brondby" => "Brøndby",
        "Copenhagen" => "København",
        "Lyngby Taarbæk" => "Lyngby Taarbæk",
        "Nordfyn" => "Nordfyns",
        "Aarhus" => "Århus"
    )

kom2reg = @_ "/home/and/data/google-mobility/Global_Mobility_Report.csv" |>
    CSV.File(__) |>
    DataFrame |>
    filter(:country_region => ==("Denmark"), __) |>
    filter( :sub_region_1 => x->!ismissing(x), __) |>
    filter( :sub_region_2 => x->!ismissing(x), __) |>
    select(__, :sub_region_1 => :region, :sub_region_2 => :kommune) |>
    disallowmissing |>
    groupby(__, :kommune) |>
    combine(__, unique) |>
    select(__,
        :region => ByRow(x -> old2new_names[x]) => :region,
        :kommune => ByRow( x->string(first(split(x, " Municipality"))) ) => :kommune,
    ) |>
    transform(__, :kommune => ByRow( x-> x ∈ keys(eng2dan) ? eng2dan[x] : x ) => :kommune ) |>
    transform(__, :kommune => ByRow( x->replace(x, "-"=>" ") ) => :kommune) |>
    append!(__,
        DataFrame( Dict(
        "kommune" => ["Ærø", "Samsø", "Fanø", "Læsø"],
        "region"  => ["south", "central", "south", "north"]
    ))) |>
    Dict([row.kommune => row.region for row in eachrow(__)])

function json2long(raw)
    num_kommunes = length(raw["locations"])
    data = Array(hcat(raw["data"]...)')
    source = string.(repeat(raw["locations"], inner=num_kommunes))
    target = string.(repeat(raw["locations"], outer=num_kommunes))
    dates = Date.(raw["dates"], DT)

    df = DataFrame(
        "date" => Vector{Date}(),
        "source" => Vector{String}(),
        "target" => Vector{String}(),
        "flow" =>  Vector{Float64}()
        )

    for (d, flow) in zip(dates, eachcol(data))
        append!(df, DataFrame( (; date = fill(d, num_kommunes^2), source, target, flow) ) )
    end

    return df
end

const DT = dateformat"yyyy-mm-dd HH:MM:SS"
telco = @_ HTTP.get("https://covid19.compute.dtu.dk/data/telco_data.json") |>
    String(__.body) |>
    JSON.Parser.parse |>
    json2long |>
    transform(__, :source => ByRow(x->kom2reg[x]) => :source) |>
    transform(__, :target => ByRow(x->kom2reg[x]) => :target) |>
    filter(_.source == _.target, __) |>
    select(__, :date, :source => :region, :flow => :telco) |>
    groupby(__, [:date, :region]) |>
    combine(__, :telco => sum => :telco)

## ============================================================================
#                    combine google, apple and telco
#  ============================================================================
mobil = @_ outerjoin(apple, google; on=[:date,:region]) |>
    outerjoin(__, telco; on=[:date,:region]) |>
    sort(__, :date) |>
    filter( Date("2020-02-15") <= _.date <= Date("2021-06-30"), __)
CSV.write(projectdir("data/mobility/", "mobility_regional.csv"), mobil)

## regional data sets

function rollingmean(xs, window)
    n  = length(xs)-window+1
    ys = Vector{Float64}(undef, n)
    for i in 1:n
        x = xs[i:i+window-1]
        @assert sum( ismissing.(x) ) < 4 "nr of missing: $(sum( ismissing.(x) )) at i=$i"
        ys[i] = mean(skipmissing(x))
        @assert isfinite(ys[i]) "infinite at i=$i"
    end
    return ys
end

numericcols(df::DataFrame) =
    filter(x->nonmissingtype( eltype( df[!,x] ) ) <: Number, names(df))

function postprocessing(df, window)
    @assert isodd(window)
    Δ = window ÷ 2
    smoothed = DataFrame( :date => df.date[1+Δ:end-Δ] )
    cols = filter(x->nonmissingtype( eltype( df[!,x] ) ) <: Number, names(df)) #[eltype(col) <: Number for col in eachcol(mobil)]
    for col in cols
        @info col
        smoothed[!,col] = rollingmean(df[:,col], window)
        smoothed[!,col] ./= abs(smoothed[1,col])
        smoothed[!,col] .-= 1
        smoothed[!,col] ./= StatsBase.std(smoothed[!,col])
        if !all( isfinite.(smoothed[!,col]) )
            @error "is not finite at $col:"
            println("$(findall(x-> !isfinite(x), smoothed[!,col]))")
        end
    end
    smoothed
end

function assertconsistency(df)
    @assert isequal( df.date |> length,
                     df.date[1]:Day(1):df.date[end] |> length)
    @assert all( [!any( @. isnan(col) | ismissing(col) | ! isfinite(col) ) for col in eachcol(df) if eltype(col) <: Number] )
    return df
end

function drop_empy_columns(df)
    drop = []
    for n in names(df)
        col = df[!,n]
        all(ismissing.(col)) && push!(drop, n)
    end
    return select(df, Not(drop))
end

regionalmobil = let
    regionalmobil = Dict()
    for region in regions
        @info region
        regional = @_ mobil |> filter(:region => ==(region), __) |>
            sort(__, :date) |>
            drop_empy_columns |>
            postprocessing(__, 7) |>
            select(__, All(), AsTable([:google,:apple]) => ByRow(mean) => :applegoogle) |>
            disallowmissing |>
            assertconsistency

        CSV.write(projectdir("data/mobility/", "mobility_region=$region.csv"), regional)
        regionalmobil[region] = regional
    end
    regionalmobil
end

## plot
using Plots
plotlyjs()
let
    ps = []
    nplots = length(Regional.regions)
    for region in regions
        mob = regionalmobil[region]
        p = plot(title="$region", xticks=:native)
        plot!(mob.date, mob.apple, lab="apple")
        plot!(mob.date, mob.google, lab="google")
        plot!(mob.date, mob.applegoogle, lab="apple-google")
        plot!(mob.date, mob.telco, lab="telco")
        push!(ps,p)
    end
    plot(ps..., layout=(nplots,1), sharex=true, link=:x, size=(1000, nplots*250))
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
int_rules_4_1_0_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.0.3*(a*csc)^m*(b*sec)^n =#
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * ~b * (~a * csc(~e + ~f * ~x)) ^ (~m - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (~n - 1)) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (EqQ((~m + ~n) - 2, 0) && NeQ(~n, 1))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) ^ ~(m') * sec(~(e') + ~(f') * ~x) ^ ~(n'), ~x) => (1 / ~f) * Subst(Antiderivative((1 + (~x) ^ 2) ^ ((~m + ~n) / 2 - 1) / (~x) ^ ~m, ~x), ~x, tan(~e + ~f * ~x)) <-- FreeQ([~e, ~f], ~x) && IntegersQ(~m, ~n, (~m + ~n) / 2)
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * sec(~(e') + ~(f') * ~x) ^ ~(n'), ~x) => (-1 / (~f * (~a) ^ ~n)) * Subst(Antiderivative((~x) ^ ((~m + ~n) - 1) / (-1 + (~x) ^ 2 / (~a) ^ 2) ^ ((~n + 1) / 2), ~x), ~x, ~a * csc(~e + ~f * ~x)) <-- FreeQ([~a, ~e, ~f, ~m], ~x) && (IntegerQ((~n + 1) / 2) && Not(IntegerQ((~m + 1) / 2) && LtQ(0, ~m, ~n)))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~m * csc(~(e') + ~(f') * ~x) ^ ~(n'), ~x) => (1 / (~f * (~a) ^ ~n)) * Subst(Antiderivative((~x) ^ ((~m + ~n) - 1) / (-1 + (~x) ^ 2 / (~a) ^ 2) ^ ((~n + 1) / 2), ~x), ~x, ~a * sec(~e + ~f * ~x)) <-- FreeQ([~a, ~e, ~f, ~m], ~x) && (IntegerQ((~n + 1) / 2) && Not(IntegerQ((~m + 1) / 2) && LtQ(0, ~m, ~n)))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-(~a) * (~a * csc(~e + ~f * ~x)) ^ (~m - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n + 1)) / (~f * ~b * (~m - 1)) + (((~a) ^ 2 * (~n + 1)) / ((~b) ^ 2 * (~m - 1))) * Antiderivative((~a * csc(~e + ~f * ~x)) ^ (~m - 2) * (~b * sec(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (GtQ(~m, 1) && (LtQ(~n, -1) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b * (~a * csc(~e + ~f * ~x)) ^ (~m + 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~f * ~a * (~n - 1)) + (((~b) ^ 2 * (~m + 1)) / ((~a) ^ 2 * (~n - 1))) * Antiderivative((~a * csc(~e + ~f * ~x)) ^ (~m + 2) * (~b * sec(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (GtQ(~n, 1) && (LtQ(~m, -1) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-(~a) * ~b * (~a * csc(~e + ~f * ~x)) ^ (~m - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (~m - 1)) + (((~a) ^ 2 * ((~m + ~n) - 2)) / (~m - 1)) * Antiderivative((~a * csc(~e + ~f * ~x)) ^ (~m - 2) * (~b * sec(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (GtQ(~m, 1) && (IntegersQ(2 * ~m, 2 * ~n) && Not(GtQ(~n, ~m))))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * ~b * (~a * csc(~e + ~f * ~x)) ^ (~m - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (~n - 1)) + (((~b) ^ 2 * ((~m + ~n) - 2)) / (~n - 1)) * Antiderivative((~a * csc(~e + ~f * ~x)) ^ ~m * (~b * sec(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (GtQ(~n, 1) && IntegersQ(2 * ~m, 2 * ~n))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~b * (~a * csc(~e + ~f * ~x)) ^ (~m + 1) * (~b * sec(~e + ~f * ~x)) ^ (~n - 1)) / (~a * ~f * (~m + ~n)) + ((~m + 1) / ((~a) ^ 2 * (~m + ~n))) * Antiderivative((~a * csc(~e + ~f * ~x)) ^ (~m + 2) * (~b * sec(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (LtQ(~m, -1) && (NeQ(~m + ~n, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-(~a) * (~a * csc(~e + ~f * ~x)) ^ (~m - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~m + ~n)) + ((~n + 1) / ((~b) ^ 2 * (~m + ~n))) * Antiderivative((~a * csc(~e + ~f * ~x)) ^ ~m * (~b * sec(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (LtQ(~n, -1) && (NeQ(~m + ~n, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~a * csc(~e + ~f * ~x)) ^ ~m * (~b * sec(~e + ~f * ~x)) ^ ~n) / tan(~e + ~f * ~x) ^ ~n) * Antiderivative(tan(~e + ~f * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ(~n)) && EqQ(~m + ~n, 0))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * csc(~e + ~f * ~x)) ^ ~m * (~b * sec(~e + ~f * ~x)) ^ ~n * (~a * sin(~e + ~f * ~x)) ^ ~m * (~b * cos(~e + ~f * ~x)) ^ ~n * Antiderivative((~a * sin(~e + ~f * ~x)) ^ -(~m) * (~b * cos(~e + ~f * ~x)) ^ -(~n), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (IntegerQ(~m - 1 / 2) && IntegerQ(~n - 1 / 2))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a) ^ 2 / (~b) ^ 2) * (~a * csc(~e + ~f * ~x)) ^ (~m - 1) * (~b * sec(~e + ~f * ~x)) ^ (~n + 1) * (~a * sin(~e + ~f * ~x)) ^ (~m - 1) * (~b * cos(~e + ~f * ~x)) ^ (~n + 1) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ -(~m) * (~b * cos(~e + ~f * ~x)) ^ -(~n), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && Not(SimplerQ(-(~m), -(~n)))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a) ^ 2 / (~b) ^ 2) * (~a * sec(~e + ~f * ~x)) ^ (~m - 1) * (~b * csc(~e + ~f * ~x)) ^ (~n + 1) * (~a * cos(~e + ~f * ~x)) ^ (~m - 1) * (~b * sin(~e + ~f * ~x)) ^ (~n + 1) * Antiderivative((~a * cos(~e + ~f * ~x)) ^ -(~m) * (~b * sin(~e + ~f * ~x)) ^ -(~n), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x)
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.0.3 (a csc)^m (b sec)^n.jl


text:
"""
Smets and Wouters 2007
Model available at: https://www.macromodelbase.com/ under "US_SW07"
Also at: https://github.com/IMFS-MMB/mmb-rep/tree/master/US_SW07
Paper: Smets, F., Wouters, R., 2007. Shocks and frictions in US business cycles: A bayesian DSGE approach. The American Economic Review 97(3), 586–606.
"""
module SW07

using ModelBaseEcon

model = Model()
model.flags.ssZeroSlope = true

@parameters model begin
    # fixed parameters
    ctou = .025;      # depreciation rate
    clandaw = 1.5;    # SS markup labor market
    cg = 0.18;        # exogenous spending GDP-ratio
    curvp = 10;       # curvature Kimball aggregator goods market
    curvw = 10;       # curvature Kimball aggregator labor market

    # estimated parameters initialisation
    ctrend = 0.4312;                  # quarterly trend growth rate to GDP
    cgamma = @link ctrend / 100 + 1;
    constebeta = 0.1657;
    cbeta = @link 100 / (constebeta + 100);     # discount factor
    constepinf = 0.7869;              # quarterly SS inflation rate
    cpie = @link constepinf / 100 + 1;
    constelab = 0.5509;

    calfa = 0.1901;           # labor share in production

    csigma = 1.3808;          # intertemporal elasticity of substitution
    cfc = 1.6064;
    cgy = 0.5187;

    csadjcost = 5.7606;      # investment adjustment cost
    chabb =     0.7133;       # habit persistence
    cprobw =    0.7061;       # calvo parameter labor market
    csigl =     1.8383;
    cprobp =    0.6523;       # calvo parameter goods market
    cindw =     0.5845;       # indexation labor market
    cindp =     0.2432;       # indexation goods market
    czcap =     0.5462;       # capital utilization
    crpi =      2.0443;       # Taylor rule reaction to inflation
    crr =       0.8103;       # Taylor rule interest rate smoothing
    cry =       0.0882;       # Taylor rule long run reaction to output gap
    crdy =      0.2247;       # Taylor rule short run reaction to output gap

    crhoa =     0.9577;
    crhob =     0.2194;
    crhog =     0.9767;
    crhoqs =    0.7113;
    crhoms =    0.1479;
    crhopinf =  0.8895;
    crhow =     0.9688;
    cmap =      0.7010;
    cmaw  =     0.8503;

    # derived from steady state
    clandap = @link cfc;
    cbetabar = @link cbeta * cgamma^(-csigma);
    cr = @link cpie / (cbeta * cgamma^(-csigma));
    crk = @link (cbeta^(-1)) * (cgamma^csigma) - (1 - ctou);
    cw = @link (calfa^calfa * (1 - calfa)^(1 - calfa) / (clandap * crk^calfa))^(1 / (1 - calfa));
    cikbar = @link (1 - (1 - ctou) / cgamma);
    cik = @link (1 - (1 - ctou) / cgamma) * cgamma;
    clk = @link ((1 - calfa) / calfa) * (crk / cw);
    cky = @link cfc * (clk)^(calfa - 1);
    ciy = @link cik * cky;
    ccy = @link 1 - cg - cik * cky;
    crkky = @link crk * cky;
    cwhlc = @link (1 / clandaw) * (1 - calfa) / calfa * crk * cky / ccy;
    cwly = @link 1 - crk * cky;
end

@variables model begin
    a; b; c; cf; dc;
    dinve; dw; dy; epinfma; ewma;
    g; inve; invef; k; kf;
    kp; kpf; lab; labf; labobs;
    mc; ms; pinf; pinf4; pinfobs;
    pk; pkf; qs; r; rk;
    rkf; robs; rrf; spinf; sw;
    w; wf; y; yf; zcap;
    zcapf;
end

@shocks model begin
    ea; eb; eg; em; epinf; eqs; ew;
end

@autoexogenize model begin
    labobs = eg
    robs = em
    pinfobs = epinf
    dy = ea
    dc = eb
    dinve = eqs
    dw = ew
end

@equations model begin

    # flexible economy
    0 * (1 - calfa) * a[t] + 1 * a[t] = calfa * rkf[t] + (1 - calfa) * (wf[t])
    zcapf[t] = (1 / (czcap / (1 - czcap))) * rkf[t]
    rkf[t] = (wf[t]) + labf[t] - kf[t]
    kf[t] = kpf[t - 1] + zcapf[t]
    "investment Euler equation"
    invef[t] = (1 / (1 + cbetabar * cgamma)) * (invef[t - 1] + cbetabar * cgamma * invef[t + 1] + (1 / (cgamma^2 * csadjcost)) * pkf[t]) + qs[t]
    pkf[t] = -rrf[t] - 0 * b[t] + (1 / ((1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)))) * b[t] + (crk / (crk + (1 - ctou))) * rkf[t + 1] + ((1 - ctou) / (crk + (1 - ctou))) * pkf[t + 1]
    "consumption Euler equation"
    cf[t] = (chabb / cgamma) / (1 + chabb / cgamma) * cf[t - 1] + (1 / (1 + chabb / cgamma)) * cf[t + 1] + ((csigma - 1) * cwhlc / (csigma * (1 + chabb / cgamma))) * (labf[t] - labf[t + 1]) - (1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)) * (rrf[t] + 0 * b[t]) + b[t]
    "aggregate resource constraint"
    yf[t] = ccy * cf[t] + ciy * invef[t] + g[t] + crkky * zcapf[t]
    "aggregate production function"
    yf[t] = cfc * (calfa * kf[t] + (1 - calfa) * labf[t] + a[t])
    wf[t] = csigl * labf[t] + (1 / (1 - chabb / cgamma)) * cf[t] - (chabb / cgamma) / (1 - chabb / cgamma) * cf[t - 1]
    "accumulation of installed capital"
    kpf[t] = (1 - cikbar) * kpf[t - 1] + (cikbar) * invef[t] + (cikbar) * (cgamma^2 * csadjcost) * qs[t]

    # sticky price - wage economy
    "marginal cost"
    mc[t] = calfa * rk[t] + (1 - calfa) * (w[t]) - 1 * a[t] - 0 * (1 - calfa) * a[t]
    "capital utilization"
    zcap[t] = (1 / (czcap / (1 - czcap))) * rk[t]
    "rental rate of capital"
    rk[t] = w[t] + lab[t] - k[t]
    "Capital installed used one period later in production"
    k[t] = kp[t - 1] + zcap[t]
    "investment Euler equation"
    inve[t] = (1 / (1 + cbetabar * cgamma)) * (inve[t - 1] + cbetabar * cgamma * inve[t + 1] + (1 / (cgamma^2 * csadjcost)) * pk[t]) + qs[t]
    "arbitrage equation for value of capital"
    pk[t] = -r[t] + pinf[t + 1] - 0 * b[t] + (1 / ((1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)))) * b[t] + (crk / (crk + (1 - ctou))) * rk[t + 1] + ((1 - ctou) / (crk + (1 - ctou))) * pk[t + 1]
    "consumption Euler equation"
    c[t] = (chabb / cgamma) / (1 + chabb / cgamma) * c[t - 1] + (1 / (1 + chabb / cgamma)) * c[t + 1] + ((csigma - 1) * cwhlc / (csigma * (1 + chabb / cgamma))) * (lab[t] - lab[t + 1]) - (1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)) * (r[t] - pinf[t + 1] + 0 * b[t]) + b[t]
    "aggregate resource constraint"
    y[t] = ccy * c[t] + ciy * inve[t] + g[t] + 1 * crkky * zcap[t]
    "aggregate production function"
    y[t] = cfc * (calfa * k[t] + (1 - calfa) * lab[t] + a[t])
    "Phillips Curve"
    pinf[t] = (1 / (1 + cbetabar * cgamma * cindp)) * (cbetabar * cgamma * pinf[t + 1] + cindp * pinf[t - 1] + ((1 - cprobp) * (1 - cbetabar * cgamma * cprobp) / cprobp) / ((cfc - 1) * curvp + 1) * (mc[t])) + spinf[t]
    w[t] = (1 / (1 + cbetabar * cgamma)) * w[t - 1] + (cbetabar * cgamma / (1 + cbetabar * cgamma)) * w[t + 1] + (cindw / (1 + cbetabar * cgamma)) * pinf[t - 1] - (1 + cbetabar * cgamma * cindw) / (1 + cbetabar * cgamma) * pinf[t] + (cbetabar * cgamma) / (1 + cbetabar * cgamma) * pinf[t + 1] + (1 - cprobw) * (1 - cbetabar * cgamma * cprobw) / ((1 + cbetabar * cgamma) * cprobw) * (1 / ((clandaw - 1) * curvw + 1)) * (csigl * lab[t] + (1 / (1 - chabb / cgamma)) * c[t] - ((chabb / cgamma) / (1 - chabb / cgamma)) * c[t - 1] - w[t]) + 1 * sw[t]

    "Monetary Policy Rule"
    r[t] = crpi * (1 - crr) * pinf[t] + cry * (1 - crr) * (y[t] - yf[t]) + crdy * (y[t] - yf[t] - y[t - 1] + yf[t - 1]) + crr * r[t - 1] + ms[t]
    a[t] = crhoa * a[t - 1] + ea[t]
    b[t] = crhob * b[t - 1] + eb[t]
    "exogenous spending (also including net exports)"
    g[t] = crhog * (g[t - 1]) + eg[t] + cgy * ea[t]
    qs[t] = crhoqs * qs[t - 1] + eqs[t]
    ms[t] = crhoms * ms[t - 1] + em[t]
    "cost push shock"
    spinf[t] = crhopinf * spinf[t - 1] + epinfma[t] - cmap * epinfma[t - 1]
    epinfma[t] = epinf[t]
    sw[t] = crhow * sw[t - 1] + ewma[t] - cmaw * ewma[t - 1]
    ewma[t] = ew[t]
    "accumulation of installed capital"
    kp[t] = (1 - cikbar) * kp[t - 1] + cikbar * inve[t] + cikbar * cgamma^2 * csadjcost * qs[t]

    # measurement equations
    dy[t] = y[t] - y[t - 1] + ctrend
    dc[t] = c[t] - c[t - 1] + ctrend
    dinve[t] = inve[t] - inve[t - 1] + ctrend
    dw[t] = w[t] - w[t - 1] + ctrend
    pinfobs[t] = 1 * (pinf[t]) + constepinf
    pinf4[t] = pinf[t] + pinf[t - 1] + pinf[t - 2] + pinf[t - 3]
    robs[t] = 1 * (r[t]) + constebeta
    labobs[t] = lab[t] + constelab
end

@initialize(model)

end


text:
using CAGD
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

V = [
	0.0 2.0 4.0 1.0 2.0 3.0 1.0 3.0 0.0 4.0
	4.0 4.0 4.0 3.0 3.0 3.0 1.0 1.0 0.0 0.0
];
EV = [
	[ 1, 2], [ 2, 3], [ 1, 9], [ 2, 5],
	[ 3,10], [ 4, 5], [ 5, 6], [ 4, 7],
	[ 6, 8], [ 7, 8], [ 9,10]
];
FV = [
    [1, 2, 2, 3, 4, 5, 5, 6, 7, 8, 9, 10],
    [            4, 5,    6, 7, 8       ],
    [1, 2,    3,                   9, 10]
]
FE = [
	[1, 2, 3, 4,  4, 5, 6, 7,  8, 9, 10, 11],
	[                   6, 7,  8, 9, 10    ],
	[1, 2, 3,        5,                  11]
]

copEV = SparseArrays.sparse(Int8.([
	1  1  0  0  0  0  0  0  0  0
	0  1  1  0  0  0  0  0  0  0
	1  0  0  0  0  0  0  0  1  0
	0  1  0  0  1  0  0  0  0  0
	0  0  1  0  0  0  0  0  0  1
	0  0  0  1  1  0  0  0  0  0
	0  0  0  0  1  1  0  0  0  0
	0  0  0  1  0  0  1  0  0  0
	0  0  0  0  0  1  0  1  0  0
	0  0  0  0  0  0  1  1  0  0
	0  0  0  0  0  0  0  0  1  1
]))

copFE = SparseArrays.sparse(Int8.([
	1  1  1  2  1  1  1  1  1  1  1
	0  0  0  0  0  1  1  1  1  1  0
	1  1  1  0  1  0  0  0  0  0  1
]))

copFV = SparseArrays.sparse(Int8.([
    1  2  1  1  2  1  1  1  1  1
    0  0  0  1  1  1  1  1  0  0
    1  1  1  0  0  0  0  0  1  1
]))

#--- Using CAGD
model = CAGD.Model(V)
CAGD.addModelCells!(model, 1, EV, signed = false)
CAGD.addModelCells!(model, 2, FE, signed = false)
copEV = model.T[1]
copFE = model.T[2]
#--- Otherwise

copFV = map(x -> floor(Int8, x / 2), copFE * copEV)

copFE = map(x -> floor(Int8, x / 2), copFV * copEV')

# EXAMPLE 2
#

V = [
    1.0 3.0  0.0 1.0 3.0 4.0  2.0  2.0  1.0 3.0
    4.0 4.0  3.0 3.0 3.0 3.0  2.0  1.0  0.0 0.0
]

copEV = SparseArrays.sparse(Int8.([
    1 0 1 0 0  0 0 0 0 0
    1 0 0 1 0  0 0 0 0 0
    0 1 0 0 1  0 0 0 0 0
    0 1 0 0 0  1 0 0 0 0
    0 0 1 1 0  0 0 0 0 0
    0 0 0 0 1  1 0 0 0 0
    0 0 0 1 0  0 1 0 0 0
    0 0 0 0 1  0 1 0 0 0
    0 0 0 0 0  0 1 1 0 0
    0 0 0 0 0  0 0 1 1 0
    0 0 0 0 0  0 0 1 0 1
    0 0 0 0 0  0 0 0 1 1
]))

copFE = SparseArrays.sparse(Int8.([
    1 1 0 0 1  0 0 0 0 0  0 0
    0 0 1 1 0  1 0 0 0 0  0 0
    0 0 0 0 0  0 0 0 0 1  1 1
    1 1 1 1 1  1 2 2 2 1  1 1
]))


text:
create_ksp() = KSPCreate()

function create_ksp(A::PetscMat)
    ksp = KSPCreate()
    KSPSetOperators(ksp, A, A)
    return ksp
end

function create_ksp(Amat::PetscMat, Pmat::PetscMat)
    ksp = KSPCreate()
    KSPSetOperators(ksp, Amat, Pmat)
    return ksp
end

set_operators!(ksp::PetscKSP, Amat::PetscMat) = KSPSetOperators(ksp, Amat, Amat)
set_operators!(ksp::PetscKSP, Amat::PetscMat, Pmat::PetscMat) = KSPSetOperators(ksp, Amat, Pmat)

solve!(ksp::PetscKSP, b::PetscVec, x::PetscVec) = KSPSolve(ksp, b, x)

function solve(ksp::PetscKSP, b::PetscVec)
    x = VecDuplicate(b)
    KSPSolve(ksp, b, x)
    return x
end

set_up!(ksp::PetscKSP) = KSPSetUp(ksp)

set_from_options!(ksp::PetscKSP) = KSPSetFromOptions(ksp)

Base.show(::IO, ksp::PetscKSP) = KSPView(ksp)

destroy!(ksp::PetscKSP) = KSPDestroy(ksp)


text:
using StochasticSemiDiscretizationMethod

function createHayesProblem(a,β)
    AMx =  ProportionalMX(a*ones(1,1));
    τ1=1.
    BMx1 = DelayMX(τ1,zeros(1,1));
    cVec = Additive(1)
    noiseID = 1
    αMx1 = stCoeffMX(noiseID,ProportionalMX(zeros(1,1)))
    βMx11 = stCoeffMX(noiseID,DelayMX(τ1,β*ones(1,1)))
    σ = stAdditive(1,Additive(ones(1)))
    LDDEProblem(AMx,[BMx1],[αMx1],[βMx11],cVec,[σ])
end

hayes_lddep=createHayesProblem(-6.,2.); # LDDE problem for Hayes equation
method=SemiDiscretization(0,0.1) # 0th order semi discretization with Δt=0.1
τmax=1. # the largest τ of the system
# Second Moment mapping
mapping=DiscreteMapping_M2(hayes_lddep,method,τmax,n_steps=10,calculate_additive=true); #The discrete mapping of the system

@show spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
statM2=VecToCovMx(fixPointOfMapping(mapping), length(mapping.M1_Vs[1])); # stationary second moment matrix of the hayes equation (equilibrium position)
@show statM2[1,1]

# Moment stability chart for the Hayes equation
using MDBM
using Plots
gr();
using LaTeXStrings

method=SemiDiscretization(0,0.1);
τmax=1.

foo(a,b) = log(spectralRadiusOfMapping(DiscreteMapping_M2(createHayesProblem(a,b),method,τmax,
    n_steps=10))); # No additive term calculated

axis=[Axis(-9.0:1.0,:a),
    Axis(-5.0:5.0,:β)]

iteration=4;
stab_border_points=getinterpolatedsolution(solve!(MDBM_Problem(foo,axis),iteration));

scatter(stab_border_points...,xlim=(-9.,1.),ylim=(-5.,5.),
    label="",title="2nd moment stability border of the Hayes equation",xlabel=L"A",ylabel=L"$\beta$",
    guidefontsize=14,tickfont = font(10),markersize=2,markerstrokewidth=0)


text:
## CD3 γϵ and δϵ cooperative binding to TCRαβ ecotodomain
##

# bimolecular dissociations
#       γϵ + TCR ← γϵ-TCR
#       δϵ + TCR ← δϵ-TCR

##
using Revise
# include("model_fit_global.jl")
include("models.jl")
include("models_dissociation.jl")
include("helper_functions.jl")
include("binning_helper_fs.jl")
include("fitting.jl")
include("helper_plotting.jl")
using DataFrames, Plots

## load thermal data

CD3_data = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/CD3_data/dissociations.csv",DataFrame)

t_γϵ_thrml = delete_missing(CD3_data[:,5:5])
t_δϵ_thrml = delete_missing(CD3_data[:,11:11])
t_mix_thrml = delete_missing(CD3_data[:,17:17])
p_γϵ_thrml = survival_p.(t_γϵ_thrml)
p_δϵ_thrml = survival_p.(t_δϵ_thrml)
p_mix_thrml = survival_p.(t_mix_thrml)

# γϵ_thrml = DissData(p_γϵ_thrml,t_γϵ_thrml,[[]])
# δϵ_thrml = DissData(p_δϵ_thrml,t_δϵ_thrml,[[]])
# mix_thrml = DissData(p_mix_thrml,t_mix_thrml,[[]])

γϵ_thrml = DissData(p_γϵ_thrml,rescale_t.(t_γϵ_thrml),[[]])
δϵ_thrml = DissData(p_δϵ_thrml,rescale_t.(t_δϵ_thrml),[[]])
mix_thrml = DissData(p_mix_thrml,rescale_t.(t_mix_thrml),[[]])

## Thermal fluctuation bimolecular dissociation

gr()

u₀ = [1.0]          # Initial condition
tspan = (0.0,10.0)  # Simulation interval
dens = fill([],length(t_γϵ_thrml))
model = bimolecular_diss!
cons_thrml = [[0.0]]

#pack parameters
γϵM_thrml = M(model,u₀,tspan,dens,cons_thrml)

opts = (0.,100.)
pinit = [0.0]

# CD3γϵ thermal fitting---------------------------------------------------------
rates_γϵ_thrml,l_γϵ_thrml = optimization(γϵM_thrml, γϵ_thrml, opts)
k₋₁ᵒ = rates_γϵ_thrml

# CD3δϵ thermal fitting---------------------------------------------------------
δϵM_thrml = M(model,u₀,tspan,dens,cons_thrml)
rates_δϵ_thrml,l_δϵ_thrml = optimization(δϵM_thrml, δϵ_thrml, opts)
k₋₂ᵒ = rates_δϵ_thrml

# CD3 mix thermal fitting-------------------------------------------------------
u₀ = [0.4,0.4,0.2]
cons_mix_thrml = [[k₋₁ᵒ[1],k₋₂ᵒ[1]]]
cons_mix_thrml = [[14.1,11.1]]
model = cd3_dissociation!

opts = (0.,100.)
pinit = [0.0,0.0]
mixM_thrml = M(model,u₀,tspan,dens,cons_mix_thrml)

rates_mix_thrml,l_mix_thrml = optimization(mixM_thrml, mix_thrml, opts)

pyplot()
plotly()
# thermal fluctuation plotting
# gr()
p1 = plot()
show_diss_fit!(p1,γϵM_thrml,rates_γϵ_thrml, γϵ_thrml; cmap=:blue)
show_diss_fit!(p1,δϵM_thrml,rates_δϵ_thrml, δϵ_thrml; cmap=:red)
show_diss_fit!(p1,mixM_thrml,rates_mix_thrml, mix_thrml; cmap=:green)

show_diss_fit!(p1,γϵM_thrml,[14.1], γϵ_thrml; cmap=:blue)
show_diss_fit!(p1,δϵM_thrml,[11.1], δϵ_thrml; cmap=:red)
u3 = rates_mix_thrml[1]
mixM_thrml.u₀ = [(1-u3)*0.610, (1-u3)*0.390, u3]
show_diss_fit!(p1,mixM_thrml,rates_mix_thrml[2], mix_thrml; cmap=:green)

## load force dependent data

γϵ_data = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/CD3_data/cd3ge-tcr-unbinned.txtbinned bond lifetime.csv",DataFrame)
δϵ_data = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/CD3_data/cd3de-tcr-unbinned.txtbinned bond lifetime.csv",DataFrame)
mix_data = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/CD3_data/cd3mix-tcr-unbinned.txtbinned bond lifetime.csv",DataFrame)

γϵ = clean_dissociation_data(γϵ_data)
δϵ = clean_dissociation_data(δϵ_data)
mix = clean_dissociation_data(mix_data)

t_γϵ, lnp_γϵ = separate_diss_data(γϵ)
t_δϵ, lnp_δϵ = separate_diss_data(δϵ)
t_mix, lnp_mix = separate_diss_data(mix)

p_γϵ = lnp2p.(lnp_γϵ)
p_δϵ = lnp2p.(lnp_δϵ)
p_mix = lnp2p.(lnp_mix)

using StatsBase
f_γϵ_avg = extract_avg_forces(γϵ)
f_δϵ_avg = extract_avg_forces(δϵ)
f_mix_avg = extract_avg_forces(mix)

γϵ = DissData(p_γϵ,t_γϵ,f_γϵ_avg)
δϵ = DissData(p_δϵ,t_δϵ,f_δϵ_avg)
mix = DissData(p_mix,t_mix,f_mix_avg)
# mix_thrml = DissData(p_mix_thrml,t_mix_thrml,[[]])

## force dependent dissociation

u₀ = [1.0, 0.0]          # Initial condition
tspan = (0.0,20.0)       # Simulation interval
dens = fill([],length(t_γϵ))
model = bimolecular_diss_bell_activated_state_w_cons!

# CD3γϵ force-dependent dissociations-------------------------------------------
cons_γϵ = [[k₋₁ᵒ[1], f] for f in f_γϵ_avg]
cons_γϵ = [[14.1, f] for f in f_γϵ_avg]
#pack parameters
γϵM = M(model,u₀,tspan,dens,cons_γϵ)

opts = (0.,100.)
pinit = zeros(7)

# # CD3γϵ force fitting
rates_γϵ,l_γϵ = optimization(γϵM, γϵ, opts)

pyplot()
p2 = plot()
show_diss_fit!(p2,γϵM,rates_γϵ, γϵ)

# CD3δϵ force-dependent dissociations-------------------------------------------

cons_δϵ = [[11.1, f] for f in f_δϵ_avg]
#pack parameters
δϵM = M(model,u₀,tspan,dens,cons_δϵ)

opts = (0.,100.)
pinit = zeros(7)

# # CD3γϵ force fitting
rates_δϵ,l_δϵ = optimization(δϵM, δϵ, opts)

p3 = plot()
show_diss_fit!(p3,δϵM,rates_δϵ, δϵ)

# plot()
# plot_F_vs_t(clean_dissociation_data(δϵ_data))

# CD3mix force-dependent dissociations-------------------------------------------

u₀ = [0.4,0.0, 0.4,0.0, 0.2,0.0]          # Initial condition
tspan = (0.0,20.0)       # Simulation interval
dens = fill([],length(t_mix))
model = cd3_tri_dissociation_bell!

k₋₃ᵒ = rates_mix_thrml[1]
cons_mix = [[vcat(k₋₁ᵒ,rates_γϵ), vcat(k₋₂ᵒ,rates_δϵ), k₋₃ᵒ,  f] for f in f_mix_avg]

#pack parameters
mixM = M(model,u₀,tspan,dens,cons_mix)

opts = (0.,100.)
opts = [(0.,0.1),(0.,100.),(0.,0.1),(0.,100.),(0.,0.1),(0.,100.),(0.,0.1)]
pinit = zeros(7)

# # CD3γϵ force fitting
rates_mix,l_mix = optimization(mixM, mix, opts)
p4 = plot()
show_diss_fit!(p4,mixM,rates_mix, mix)
show_diss_fit!(p4,mixM, zeros(7), mix)

ft1=plot()
plot_F_vs_t(clean_dissociation_data(γϵ_data))
title!("γϵ")
ft2=plot()
plot_F_vs_t(clean_dissociation_data(δϵ_data))
title!("δϵ")
ft3=plot()
plot_F_vs_t(clean_dissociation_data(mix_data))
title!("mix")
plot(ft1,ft2,ft3)


text:
ENV["JULIA_MARGO_LOAD_PYPLOT"] = "no thank you"
using JuMP
using ClimateMARGO
using ClimateMARGO.Models, ClimateMARGO.Optimization, ClimateMARGO.Diagnostics
using Test
using BenchmarkTools

@testset "Temperature optimization" begin
    @testset "Temp goal: $(temp_goal)" for temp_goal in 1.5:0.5:4.0
        model = ClimateModel(deepcopy(ClimateMARGO.IO.included_configurations["default"]))
        status = optimize_controls!(model, temp_goal=temp_goal)
        @test raw_status(status) == "Solve_Succeeded"
        @test isapprox(
            maximum(T_adapt(model, M=true, R=true, G=true, A=true)),
            temp_goal,
            rtol=1.e-5
        )
    end
end

function default_parameters(dt=12)::ClimateModelParameters
    p = deepcopy(ClimateMARGO.IO.included_configurations["default"])

    p.domain = Domain(Float64(dt), 2020.0, 2200.0)
    p.economics.baseline_emissions = ramp_emissions(p.domain)
    p.economics.extra_CO₂ = zeros(size(p.economics.baseline_emissions))

    return p
end

@testset "Pinned results" begin
    # We generate results using ClimateMARGO.jl once, and then add a test to check that these same results are still generated. This allows us to "pin" a correct state of the model, and make internal changes with the guarantee that the outcome does not change.

    # In the future, to "pin" new results as the correct ones, run the tests, and copy the expected result from the test failure report.

    @testset "Forward" begin
        model = ClimateModel(default_parameters(20))

        bump() = [0.0:0.2:0.9;  1.0:-0.2:0.1] # == [0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 0.8, 0.6, 0.4, 0.2]
        slope() = 0.0:0.1:0.99
        constant() = fill(1.0, size(slope())...)

        model.controls.mitigate = bump()
        model.controls.remove = 0.2 .* bump()
        model.controls.geoeng = 0.1 .* slope()
        model.controls.adapt = 0.3 .* constant()

        result_T = T_adapt(model; M=true, R=true, G=true, A=true)

        @test result_T ≈ [1.2767663110173484, 1.608852123606971, 1.8978258561722958, 2.108526090432568, 2.211931788367322, 2.1713210884072347, 2.158586858753497, 2.1224302958426446, 2.0894535810555555, 2.0607413902621032] rtol=1e-5
    end

    @testset "Optimization" begin
        model = ClimateModel(default_parameters(20))

        optimize_controls!(model, temp_goal=2.2)

        @test model.controls.mitigate ≈ [0.0, 0.3670584821282429, 0.4393272792260638, 0.5257940783265341, 0.6245380250211872, 0.7074912042403808, 0.7676867396677578, 0.0, 0.0, 0.0] rtol=1e-3
        @test model.controls.remove ≈ [0.0, 0.10882651060773696, 0.13025296274983308, 0.15588887768405002, 0.18516475517861308, 0.20975894242613452, 0.22760588089905764, 0.2367501073337253, 0.23175978166812158, 0.19659824324032557] rtol=1e-3
        @test model.controls.geoeng ≈ [0.0, 0.05049106345373785, 0.05264838036542512, 0.06460796652594906, 0.10819144740515446, 0.12584271884062354, 0.12767445817227918, 0.12318301615805025, 0.11794556458720538, 0.1132343245663971] rtol=1e-3
        @test model.controls.adapt ≈ [0.0, 0.0001537202313707738, 0.00015374376329529388, 0.03816957588616602, 0.10942294994448329, 0.13786594819465353, 0.14363875559967432, 0.14044599067883523, 0.1365010682189578, 0.13366455981096753] rtol=1e-3
    end

end

const RUN_BENCHMARKS = true

if RUN_BENCHMARKS
    @info "Running benchmark..."

    function go()
        model = ClimateModel(default_parameters(20))

        optimize_controls!(model; temp_goal=2.2, print_raw_status=false)
    end

    go()

    display(@benchmark go())
end


text:
#=
@author Alvaro Lucas Castillo Calabacero
contact alvarolucascc96@gmail.com
=#

#=
Recursive algorithm for resolve MergeSort
=#
function mergeSort!(array, start = 1, finish = length(array))
    if start < finish
        middle = div(start+finish, 2)
        mergeSort!(array, start, middle)
        mergeSort!(array, middle+1, finish)
        merge(array, start, middle, finish)
    end
    return array
end

#=
algorithm for merge
 =#
function merge(array, start, half, finish)
    ward = typemax(eltype(array))
    leftSubArray = array[start:half]
    rightSubArray = array[(half+1):finish]
    push!(leftSubArray, ward)
    push!(rightSubArray, ward)
    auxLeft=1
    auxRight = 1
    for i in start:finish
      if leftSubArray[auxLeft] <= rightSubArray[auxRight]
          array[i] = leftSubArray[auxLeft]
          auxLeft += 1
      else
          array[i] = rightSubArray[auxRight]
          auxRight += 1
      end
    end
end
#can change this total numbers for process of sort
totalNumbers = 10000000;
list = [rem(rand(Int128),100000) for i =1:totalNumbers]

print("Sorted List in :")
#get time to solve whit secuencial algorithm
@time mergeSort!(list)
print(list)