text:
# Some utility functions to keep code clean in statespace_functions.jl

function compute_gensys_gensys2_regimes(m::AbstractDSGEModel)

    regime_switching = haskey(get_settings(m), :regime_switching) && get_setting(m, :regime_switching)
    n_regimes        = (regime_switching && haskey(get_settings(m), :n_regimes)) ? get_setting(m, :n_regimes) : 1

    if haskey(get_settings(m), :gensys2) && get_setting(m, :gensys2)
        if haskey(get_settings(m), :regime_eqcond_info)
            sorted_eqcond = sort!(collect(get_setting(m, :regime_eqcond_info)), by = x -> x[1]) # x -> x[1] sorts by keys
            first_gensys2_ind = findfirst(x -> x[2].alternative_policy.key in # finds the first value (EqcondEntry) whose key
                                          get_setting(m, :temporary_altpolicy_names), sorted_eqcond) # is a temporary altpolicy
            first_gensys2_regime = !isnothing(first_gensys2_ind) ? sorted_eqcond[first_gensys2_ind][1] : nothing
        elseif haskey(get_settings(m), :gensys_regimes) && haskey(get_settings(m), :gensys2_regimes)
            return get_setting(m, :gensys_regimes), get_setting(m, :gensys2_regimes)
        else
            first_gensys2_regime = nothing
        end

        if isnothing(first_gensys2_regime)
            # throw(GensysError("No equilibrium conditions in any regime are being temporarily replaced, " *
            #                   "but the setting :gensys2 is true."))
            gensys2_regimes = Vector{UnitRange{Int}}(undef, 0)
            gensys_regimes  = UnitRange{Int}[1:n_regimes]
        else
            last_gensys2_regime = haskey(get_settings(m), :temporary_altpolicy_length) ?
                min(first_gensys2_regime + get_setting(m, :temporary_altpolicy_length), n_regimes) :
                n_regimes # NOTE removed a +1 here--if tests start failing, check here first

            gensys_regimes = UnitRange{Int}[1:(first_gensys2_regime - 1)]
            if last_gensys2_regime != n_regimes
                append!(gensys_regimes, [(last_gensys2_regime + 1):n_regimes])
            end
            gensys2_regimes = [first_gensys2_regime-1:last_gensys2_regime]
        end
    else
        gensys2_regimes = Vector{UnitRange{Int}}(undef, 0)
        gensys_regimes  = UnitRange{Int}[1:n_regimes]
    end

    return gensys_regimes, gensys2_regimes
end

"""
```
compute_system_function(system::System{S}) where S<:AbstractFloat
```

### Inputs

- `system::System`: The output of compute_system(m), i.e. the matrix outputs from solving a given model, m.

### Outputs

- `Φ::Function`: transition equation
- `Ψ::Function`: measurement equation
- `F_ϵ::Distributions.MvNormal`: shock distribution
- `F_u::Distributions.MvNormal`: measurement error distribution
"""
function compute_system_function(system::System{S}) where S<:AbstractFloat
    # Unpack system
    TTT    = system[:TTT]
    RRR    = system[:RRR]
    CCC    = system[:CCC]
    QQ     = system[:QQ]
    ZZ     = system[:ZZ]
    DD     = system[:DD]
    EE     = system[:EE]

    # Define transition and measurement functions
    @inline Φ(s_t1::Vector{S}, ϵ_t::Vector{S}) = TTT*s_t1 + RRR*ϵ_t + CCC
    @inline Ψ(s_t::Vector{S}) = ZZ*s_t + DD

    # Define shock and measurement error distributions
    nshocks = size(QQ, 1)
    nobs    = size(EE, 1)
    F_ϵ = Distributions.MvNormal(zeros(nshocks), QQ)
    F_u = Distributions.MvNormal(zeros(nobs),    EE)

    return Φ, Ψ, F_ϵ, F_u
end

function zero_system_constants(system::System{S}) where S<:AbstractFloat
    system = copy(system)

    system.transition.CCC = zeros(size(system[:CCC]))
    system.measurement.DD = zeros(size(system[:DD]))
    system.pseudo_measurement.DD_pseudo = zeros(size(system[:DD_pseudo]))

    return system
end

function zero_system_constants(system::RegimeSwitchingSystem{S}) where S<:AbstractFloat
    system = copy(system)

    for i in 1:n_regimes(system)
        system.transitions[i].CCC = zeros(size(system[i, :CCC]))
        system.measurements[i].DD = zeros(size(system[i, :DD]))
        system.pseudo_measurements[i].DD_pseudo = zeros(size(system[i, :DD_pseudo]))
    end

    return system
end

mutable struct ForwardExpectationsMemo{T}
    time_varying_memo::Dict{Int64, T} # for products of time-varying TTT matrices
    permanent_memo::Dict{Int64, T} # for products of non time-varying TTT matrices
end

"""
```
ForwardExpectationsMemo(TTTs::Vector{<: AbstractMatrix{S}},
                        current_regime::Int64, last_tv_period::Int64,
                        first_perm_period::Int64, min_perm_power::Int64 = 0,
                        max_perm_power::Int64 = 0) where {S <: Real}
```
computes the memo dictionaries of the necessary products/powers of TTTs
for computing forward expectations of states and sums of states.

### Inputs
- `TTTs`: complete sequence of time-varying transition matrices from regime 1 to the final regime
- `current_regime`: the current period's regime
- `last_tv_period`: the last period in which the transition matrix is believed to have changed relative to the previous period.
- `first_perm_period`: the first period in which the transition matrix in `TTTs` is permanently imposed. This period
    should be from the perspective of an omniscient econometrician rather than the agents' perspective.
- `min_perm_power`: minimum power of the permanent matrix to be calculated, i.e. we calculate at least `TTTs[first_perm_period] ^ min_perm_power`
- `max_perm_power`: maximum power of the permanent matrix to be calculated, i.e. we calculate at most `TTTs[first_perm_period] ^ max_perm_power`

### Notes
- To clarify what `last_tv_period` should be, note that
if every matrix in `TTTs` is time-varying, AND agents believe every matrix is time-varying, then `last_tv_period = length(TTTs)`.
However, if agents are myopic and believe that the last time-varying matrix is the second to last one, then
`last_tv_period = length(TTTs) - 1`. This flexibility allows the user to specify different degrees of awareness about `TTTs`.

- Note that when `last_tv_period == first_perm_period`, we do not calculate

```
time_varying_memo[last_tv_period] = TTTs[first_perm_period] * TTTs[first_perm_period - 1] * ... TTTs[t + 1]
```

Instead, we calculate

```
time_varying_memo[last_tv_period] = TTTs[first_perm_period - 1] * ... TTTs[t + 1]
```

The reason is that it is more efficient to include that first `TTTs[first_perm_period]` in the powers of `TTT[first_perm_period]`
computed for `permanent_memo`. So if `length(TTTs) == 7`, to get the correct `k`-periods ahead forward
expectations from some period `t < 7`, you need to run

```
memo = ForwardExpectationsMemo(TTTs, t, 7, 7, 1, t + k + 1 - 7)
# or equivalently . . .
# memo = ForwardExpectationsMemo(TTTs, t, length(TTTs), length(TTTs), 1, t + k + 1 - length(TTTs))
```
"""
function ForwardExpectationsMemo(TTTs::Vector{<: AbstractMatrix{S}},
                                 current_regime::Int64, last_tv_period::Int64,
                                 first_perm_period::Int64, min_perm_power::Int64 = 0,
                                 max_perm_power::Int64 = 0) where {S <: Real}

    @assert current_regime <= last_tv_period
    @assert min_perm_power <= max_perm_power
    @assert last_tv_period <= first_perm_period

    tv_memo   = Dict{Int64, eltype(TTTs)}()
    perm_memo = Dict{Int64, eltype(TTTs)}()

    # When the last time-varying period equals the first period in which the
    # transition matrix is permanently imposed, we only calculate the products
    # of the time-varying matrices and ignore the permanent matrix
    if last_tv_period == first_perm_period
        last_tv_period -= 1
    end
    tv_memo[last_tv_period] = TTTs[last_tv_period]
    for k in (last_tv_period - 1):-1:(current_regime + 1)
        tv_memo[k] = tv_memo[k + 1] * TTTs[k]
    end

    if !(min_perm_power == max_perm_power == 0)
        TTT_perm                  = TTTs[first_perm_period]
        perm_memo[min_perm_power] = TTT_perm ^ min_perm_power
        for k in (min_perm_power + 1):max_perm_power
            perm_memo[k] = perm_memo[k - 1] * TTT_perm
        end
    end

    return ForwardExpectationsMemo(tv_memo, perm_memo)
end

mutable struct ForwardMultipleExpectationsMemo{T}
    time_varying_memo::Dict{Int64, Dict{Int64, T}} # for products of time-varying TTT matrices
    permanent_memo::Dict{Int64, T} # for products of non time-varying TTT matrices
end

"""
```
ForwardMultipleExpectationsMemo(TTTs::Vector{<: AbstractMatrix{S}},
                                current_regime::Int64, min_last_tv_period::Int64, max_last_tv_period::Int64,
                                first_perm_period::Int64, min_perm_power::Int64 = 0,
                                max_perm_power::Int64 = 0) where {S <: Real}
```
computes the memo dictionaries of the necessary products/powers of TTTs
for computing forward expectations of states and sums of states for
multiple different horizon lengths for the "forward-looking" window, i.e. variation
in how many periods ahead for which expectations are taken. This approach
avoids repeating computations that would be incurred if the user
instead repeatedly created `ForwardExpectationsMemo` for each horizon length.

### Inputs
- `TTTs`: complete sequence of time-varying transition matrices from regime 1 to the final regime
- `current_regime`: the current period's regime
- `min_last_tv_period`: the minimum last period in which the transition matrix is believed to have changed relative to the previous period.
- `max_last_tv_period`: the maximum last period in which the transition matrix is believed to have changed relative to the previous period.
- `first_perm_period`: the first period in which the transition matrix in `TTTs` is permanently imposed. This period
    should be from the perspective of an omniscient econometrician rather than the agents' perspective.
- `min_perm_power`: minimum power of the permanent matrix to be calculated, i.e. we calculate at least `TTTs[first_perm_period] ^ min_perm_power`
- `max_perm_power`: maximum power of the permanent matrix to be calculated, i.e. we calculate at most `TTTs[first_perm_period] ^ max_perm_power`

### Notes
- To clarify what `min_last_tv_period` and `max_last_tv_period` should be, first consider the meaning of input argument
`last_tv_period` for `ForwardExpectationsMemo`.
If every matrix in `TTTs` is time-varying, AND agents believe every matrix is time-varying, then `last_tv_period = length(TTTs)`.
However, if agents are myopic and believe that the last time-varying matrix is the second to last one, then
`last_tv_period = length(TTTs) - 1`. This flexibility allows the user to specify different degrees of awareness about `TTTs`.

The argument `min_last_tv_period` and `max_last_tv_period` allow the user to construct one memo object for
computing varying lengths of forward expectations. For example, suppose we want expectations of the nominal rate
from 1 to 6 periods ahead in the measurement equation in the regime `t`. Then you should run

```
ForwardExpectationsMemo(TTTs, t, t + 1, t + 6, ...) # ellipsis omits the remaining input args
```

- Note that when `max_last_tv_period == first_perm_period`,

```
time_varying_memo[max_last_tv_period] = TTTs[first_perm_period] * TTTs[first_perm_period - 1] * ... TTTs[t + 1]
```

Instead, we calculate

```
time_varying_memo[max_last_tv_period] = TTTs[first_perm_period - 1] * ... TTTs[t + 1]
```

The reason is that it is more efficient to include that first `TTTs[first_perm_period]` in the powers of `TTT[first_perm_period]`
computed for `permanent_memo`.
"""
function ForwardMultipleExpectationsMemo(TTTs::Vector{<: AbstractMatrix{S}},
                                current_regime::Int64, min_last_tv_period::Int64, max_last_tv_period::Int64,
                                first_perm_period::Int64, min_perm_power::Int64 = 0,
                                max_perm_power::Int64 = 0) where {S <: Real}

    @assert current_regime <= min_last_tv_period <= max_last_tv_period
    @assert min_perm_power <= max_perm_power
    @assert max_last_tv_period <= first_perm_period

    tv_memo   = Dict{Int64, Dict{Int64, eltype(TTTs)}}()
    perm_memo = Dict{Int64, eltype(TTTs)}()

    # Iterate over the length of the window of the forward-looking period
    for window_len in min_last_tv_period:max_last_tv_period
        tv_memo[window_len] = Dict{Int64, eltype(TTTs)}()

        # When the last time-varying period equals the first period in which the
        # transition matrix is permanently imposed, we only calculate the products
        # of the time-varying matrices and ignore the permanent matrix
        last_tv_period = window_len == first_perm_period ? window_len - 1 : window_len
        tv_memo[window_len][last_tv_period] = TTTs[last_tv_period]
        for k in (last_tv_period - 1):-1:(current_regime + 1)
            tv_memo[window_len][k] = tv_memo[window_len][k + 1] * TTTs[k]
        end
    end
    if !(min_perm_power == max_perm_power == 0)
        TTT_perm                  = TTTs[first_perm_period]
        perm_memo[min_perm_power] = TTT_perm ^ min_perm_power
        for k in (min_perm_power + 1):max_perm_power
            perm_memo[k] = perm_memo[k - 1] * TTT_perm
        end
    end

    return ForwardMultipleExpectationsMemo(tv_memo, perm_memo)
end

"""
```
k_periods_ahead_expectations(TTT, CCC, TTTs, CCCs, t, k; permanent_t = length(TTTs),
                             integ_series = false, memo = nothing)
```

calculates the matrices associated with the expected state `k` periods ahead from `t`.
This function should NOT be used with linear state space system matrices with any unit roots.

The `TTT` and `CCC` inputs are the transition matrix and constant vector associated with
the current period `t`, while the `TTTs` and `CCCs` are vectors containing the time-varying
transition matrices and constant vectors, such that `TTTs[t]` retrieves the time-varying
transition matrix associated with period `t` and `TTTs[t + k]` retrieves the time-varying
transition matrix associated with period `t + k`. The optional argument `permanent_t`
indicates the period for which the matrices/vectors are no longer time-varying, i.e.
if `t >= permanent_t`, then `TTTs[permanent_t]` is the transition matrix.

The formula implemented by this function is
```
𝔼ₜ[sₜ₊ₖ] = (∏ⱼ=₁ᵏ Tₜ₊ⱼ) sₜ + (∑ₘ₌₁ᵏ⁻¹ (∏ⱼ₌ₘ₊₁ᵏ Tₜ₊ⱼ) Cₜ₊ₘ) + Cₜ₊ₖ.
```
Additional simplifications are made if it is known that `t + k > permanent_t`
since this implies some matrices are the same. This recognition reduces
unnecessary computations.

### Keyword Arguments
- `integ_series::Bool`: set to true if there are some transition matries in `TTT`
    that result in integrated series, in which case we cannot speed up
    computations by using left-divides.

- `memo::Union{ForwardExpectationsMemo, Nothing}`: pass a properly formed
    `ForwardExpectationsMemo` to avoid calculating unnecessary products and powers of
    the matrices in `TTTs`. Typically, the memo you want to compute is

```
# min_t is minimum t you will use, maximum_t is maximum t you will use, and
# max_k is the maximum window for forward expectations.
memo = ForwardExpectationsMemo(TTTs, min_t, length(TTTs), length(TTTs), min_t + max_k - length(TTTs),
                               max_t + max_k + 1 - length(TTTs))
```
"""
function k_periods_ahead_expectations(TTT::AbstractMatrix, CCC::AbstractVector,
                                      TTTs::Vector{<: AbstractMatrix}, CCCs::Vector{<: AbstractVector},
                                      t::Int, k::Int, permanent_t::Int = length(TTTs);
                                      integ_series::Bool = false,
                                      memo::Union{ForwardExpectationsMemo, Nothing} = nothing)
# TODO: add a version that infers the case of a memo with kwarg type ForwardMultipleExpectationsMemo
# and translate it to the correct ForwardExpectationsMemo based on t, k, permanent_t

    if isempty(TTTs) || isempty(CCCs)
        Tᵏ = TTT^k
        if all(CCC .≈ 0.)
            return Tᵏ, CCC
        else
            if integ_series # are there integrated series? If so, we have to accumulate instead of doing a left-division
                T_memo = Dict{Int, typeof(TTT)}()
                T_memo[1] = TTT
                for i in 2:k
                    T_memo[i] = T_memo[i - 1] * TTT
                end
                Tᵏsum = sum([T_memo[j] for j in 1:k])
            else
                Tᵏsum = (I - TTT) \ (I - Tᵏ)
            end
            return Tᵏ, Tᵏsum * CCC
        end
    else
        if t + k <= permanent_t
            # Cannot save computation speed by not calculating further times b/c always time-varying
            t_plus_k = t + k
            if isnothing(memo)
                T_memo = Dict{Int, eltype(TTTs)}()
                T_memo[k] = TTTs[t_plus_k]
                for i in (k-1):-1:1
                    T_memo[i] = T_memo[i + 1] * TTTs[t + i]
                end
            end

            C_accum = deepcopy(CCCs[t_plus_k])
            for i in 1:(k - 1)
                if isnothing(memo)
                    C_accum .+= T_memo[i + 1] * CCCs[t + i]
                elseif t_plus_k == permanent_t
                    # This is an edge case b/c the memo computation will not add
                    # the permanent matrix to the product (and instead stores it
                    # in the memo.permanent_memo dict)
                    if i == k - 1
                        C_accum .+= TTTs[permanent_t] * CCCs[t + i]
                    else
                        C_accum .+= TTTs[permanent_t] * memo.time_varying_memo[t + (i + 1)] * CCCs[t + i]
                    end
                else
                    C_accum .+= memo.time_varying_memo[t + (i + 1)] * CCCs[t + i]
                end
            end

            if isnothing(memo)
                return T_memo[1], C_accum
            elseif t_plus_k == permanent_t
                if k == 1
                    return TTTs[permanent_t], C_accum
                else
                    return TTTs[permanent_t] * memo.time_varying_memo[t + 1], C_accum
                end
            else
                return memo.time_varying_memo[t + 1], C_accum
            end
        elseif t == permanent_t
            # None of the matrices are time-varying anymore
            if all(CCCs[permanent_t] .≈ 0.)
                # TODO: we may want to add a way to figure out if the power requested here
                # could be pre-computed with the memo.
                return TTTs[permanent_t]^k, CCCs[permanent_t]
            else
                if isnothing(memo)
                    Tᵏₜ₊₁ = TTTs[permanent_t]^k
                else
                    Tᵏₜ₊₁ = memo.permanent_memo[k]
                end
                if integ_series
                    T_memo = Dict{Int, eltype(TTTs)}()
                    T_memo[1] = TTTs[permanent_t]
                    for i in 2:k
                        T_memo[i] = T_memo[i - 1] * TTTs[permanent_t]
                    end
                    Tᵏsum = sum([T_memo[j] for j in 1:k])
                else
                    Tᵏsum = (I - TTTs[permanent_t]) \ (I - Tᵏₜ₊₁)
                end
                return Tᵏₜ₊₁, Tᵏsum * CCCs[permanent_t]
            end
        else
            # Computation time can be saved by realizing some matrices are not time-varying
            h = (permanent_t - 1) - t # last time of time-variation is permanent_t - 1
            if isnothing(memo)
                Tᵏ⁻ʰₜ₊ₕ₊₁ = (k == h) ? Diagonal(ones(length(CCCs[permanent_t]))) : TTTs[permanent_t]^(k - h) # why using Diagonal instead of I?
            else
                # Tᵏ⁻ʰₜ₊ₕ₊₁ = (k == h) ? I : memo.permanent_memo[k - h] # why using Diagonal instead of I?
                Tᵏ⁻ʰₜ₊ₕ₊₁ = (k == h) ? Diagonal(ones(length(CCCs[permanent_t]))) : memo.permanent_memo[k - h] # why using Diagonal instead of I?
            end

            if isnothing(memo)
                T_memo = Dict{Int, eltype(TTTs)}()
            end

            if h > 0
                if isnothing(memo)
                    T_memo[h] = TTTs[t + h] # maps i to ∏ⱼ₌ᵢʰ Tₜ₊ⱼ, so T_memo[h] = Tₜ₊ₕ, T_memo[h-1] = Tₜ₊ₕ * Tₜ₊ₕ₋₁, ...
                    for i in (h-1):-1:1
                        T_memo[i] = T_memo[i + 1] * TTTs[t + i]
                    end
                    T_accum = Tᵏ⁻ʰₜ₊ₕ₊₁ * T_memo[1]
                else
                    T_accum = Tᵏ⁻ʰₜ₊ₕ₊₁ * memo.time_varying_memo[t + 1]
                end
            else
                T_accum = Tᵏ⁻ʰₜ₊ₕ₊₁
            end

            if h == 0 # Nothing to accumulate from the past
                C_accum = zeros(size(T_accum, 1))
            else
                C_accum = deepcopy(CCCs[t + h])
                for i in 1:(h - 1)
                    if isnothing(memo)
                        C_accum .+= T_memo[i + 1] * CCCs[t + i]
                    else
                        C_accum .+= memo.time_varying_memo[t + i + 1] * CCCs[t + i]
                    end
                end
                C_accum .= Tᵏ⁻ʰₜ₊ₕ₊₁ * C_accum
            end

            if all(CCCs[permanent_t] .≈ 0.)
                return T_accum, C_accum
            else
                if integ_series
                    T_memo = Dict{Int, eltype(TTTs)}()
                    T_memo[1] = TTTs[permanent_t]
                    for i in 2:(k - h - 1)
                        T_memo[i] = T_memo[i - 1] * TTTs[permanent_t]
                    end
                    Tᵏ⁻ᵐsum = sum([T_memo[j] for j in 1:(k - h - 1)]) + I
                else
                    Tᵏ⁻ᵐsum = (I - TTTs[permanent_t]) \ (I - Tᵏ⁻ʰₜ₊ₕ₊₁)
                end
                C_accum .+= Tᵏ⁻ᵐsum * CCCs[permanent_t]
                return T_accum, C_accum
            end
        end
    end
end

"""
```
k_periods_ahead_expected_sums(TTT, CCC, TTTs, CCCs, t, k; permanent_t = length(TTTs),
                              integ_series = false, memo = nothing)
```

calculates the matrices associated with the sum of the expected states between periods
`t + 1` and  `t + k`. This function should NOT be used with
linear state space system matrices with any unit roots.

The `TTT` and `CCC` inputs are the transition matrix and constant vector associated with
the current period `t`, while the `TTTs` and `CCCs` are vectors containing the time-varying
transition matrices and constant vectors, such that `TTTs[t]` retrieves the time-varying
transition matrix associated with period `t` and `TTTs[t + k]` retrieves the time-varying
transition matrix associated with period `t + k`. The optional argument `permanent_t`
indicates the period for which the matrices/vectors are no longer time-varying, i.e.
if `t >= permanent_t`, then `TTTs[permanent_t]` is the transition matrix.

The formula implemented by this function is
```
∑ⱼ₌₁ᵏ 𝔼ₜ[sₜ₊ⱼ] = ∑ⱼ₌₁ᵏ(∏ⱼ=₁ᵏ Tₜ₊ⱼ) sₜ + ∑ᵣ₌₁ᵏ⁻¹(I + ∑ⱼ₌ᵣ₊₁ᵏ (∏ₘ₌ᵣ₊₁ʲ Tₜ₊ₘ))Cₜ₊ᵣ + Cₜ₊ₖ.
```
Additional simplifications are made if it is known that `t + k > permanent_t`
since this implies some matrices are the same. This recognition reduces
unnecessary computations.

### Keyword Arguments
- `integ_series::Bool`: set to true if there are some transition matries in `TTT`
    that result in integrated series, in which case we cannot speed up
    computations by using left-divides.

- `memo::Union{ForwardMultipleExpectationsMemo, Nothing}`: pass a properly formed
    `ForwardExpectationsMemo` to avoid calculating unnecessary products and powers of
    the matrices in `TTTs`. Typically, the memo you want to compute is

```
# min_t is minimum t you will use, maximum_t is maximum t you will use, and
# max_k is the maximum window for forward expectations.
memo = ForwardMultipleExpectationsMemo(TTTs, min_t, length(TTTs), length(TTTs), min_t + max_k - length(TTTs),
                               max_t + max_k + 1 - length(TTTs))
```
"""
function k_periods_ahead_expected_sums(TTT::AbstractMatrix, CCC::AbstractVector,
                                       TTTs::Vector{<: AbstractMatrix}, CCCs::Vector{<: AbstractVector},
                                       t::Int, k::Int, permanent_t::Int = length(TTTs);
                                       integ_series::Bool = false,
                                       memo::Union{ForwardMultipleExpectationsMemo, Nothing} = nothing)

    if integ_series # Do this by directly summing the k-periods ahead expectations. Not fully efficient but also not the typical case
        T_accum = Vector{eltype(TTTs)}(undef, k)
        C_accum = Vector{eltype(CCCs)}(undef, k)
        for i in 1:k
            # Create ForwardExpectationsMemo corresponding for computing the states
            # i-periods ahead from time t. We then take min(t + i, permanent_t)
            # b/c each key of memo.time_varying_memo specifies the last time-varying period,
            # with a minimum of permanent_t.
            #
            # Note that this construction is non-allocating b/c we're just wrapping
            # precomputed dictionaries, so data shouldn't be copied.
            ith_memo = isnothing(memo) ? nothing : ForwardExpectationsMemo(memo.time_varying_memo[min(t + i, permanent_t)],
                                                                           memo.permanent_memo)
            T_accum[i], C_accum[i] = k_periods_ahead_expectations(TTT, CCC, TTTs, CCCs, t, i, permanent_t,
                                                                  integ_series = integ_series,
                                                                  memo = ith_memo)
        end
        return sum(T_accum), sum(C_accum)
    end

    if isempty(TTTs) || isempty(CCCs)
        Tᵏsum = (I - TTT) \ (TTT - TTT^(k + 1))
        if all(CCC .≈ 0.)
            return Tᵏsum, CCC
        else
            TTTʲmemo = Vector{typeof(TTT)}(undef, k)
            TTTʲmemo[1] = TTT
            for q in 2:k
                TTTʲmemo[q] = TTTʲmemo[q - 1] * TTT
            end
            return Tᵏsum, (I + sum([(I - TTT) \ (I - TTTʲmemo[k - q + 1]) for q in 1:(k - 1)])) * CCC
        end
    else
        if t + k <= permanent_t
            # Cannot save computation speed by not calculating further times b/c always time-varying
            total_Tsum = zeros(eltype(TTTs[t]), size(TTTs[t]))
            total_Csum = zeros(eltype(CCCs[t]), size(CCCs[t]))

            # Time-varying matrices so need to calculate in loop
            if isnothing(memo)
                # T_accum and C_accum store expected values in period i for some i
                T_accum = TTTs[t+1]
                C_accum = CCCs[t+1]
                total_Tsum .+= T_accum
                total_Csum .+= C_accum

                for i in (t+2):(t+k)
                    T_accum = TTTs[i] * T_accum
                    total_Tsum .+= T_accum

                    C_accum = TTTs[i] * C_accum + CCCs[i]
                    total_Csum .+= C_accum
                end
            else
                for i in 1:k
                    ith_memo = ForwardExpectationsMemo(memo.time_varying_memo[min(t + i, permanent_t)],
                                                       memo.permanent_memo) # see comments above for why this is correct
                    tmp1, tmp2 = k_periods_ahead_expectations(TTT, CCC, TTTs, CCCs, t, i, permanent_t,
                                                              integ_series = integ_series, memo = ith_memo)
                    total_Tsum .+= tmp1
                    total_Csum .+= tmp2
                    # total_Tsum .+= i == k ? memo.time_varying_memo[t+1] : memo.time_varying_memo[1] / memo.time_varying_memo[t+i]
                    # total_Csum .+= (I + memo.time_varying_memo[t+i]) * CCCs[t+i]
                end
            end

            return total_Tsum, total_Csum

        elseif t == permanent_t
            # None of the matrices are time-varying anymore
            Tᵏsum = (I - TTTs[permanent_t]) \ (TTTs[permanent_t] - TTTs[permanent_t]^(k + 1))
            if all(CCC .≈ 0.)
                return Tᵏsum, CCCs[permanent_t]
            else
                TTTʲmemo = Vector{eltype(TTTs)}(undef, k)
                TTTʲmemo[1] = TTTs[permanent_t]
                for q in 2:k
                    TTTʲmemo[q] = TTTʲmemo[q - 1] * TTTs[permanent_t]
                end
                return Tᵏsum, (I + sum([(I - TTTs[permanent_t]) \ (I - TTTʲmemo[k - q + 1]) for q in 1:(k - 1)])) * CCCs[permanent_t]
            end
        else
            # Computation time can be saved by realizing some matrices are not time-varying
            total_Tsum = zeros(eltype(TTTs[t]), size(TTTs[t]))
            total_Csum = zeros(eltype(CCCs[t]), size(CCCs[t]))

            # Do we need to adjust for C after permanent_t?
            do_C = !(all(CCCs[permanent_t] .≈ 0.0))

            # T_accum and C_accum store the expectation of T and C at time i
            ## So, T_accum at time i will store ∏_{j=t+1}^i T_j
            ## C_accum at time i stores ∑_{q=t+1}^i C_q + ∑_{j=t+1}^q T_j C_q
            # total_sum will stores sum over the expectations

            if isnothing(memo)
                T_accum = TTTs[t+1]
                C_accum = CCCs[t+1]
                total_Tsum .+= T_accum
                total_Csum .+= C_accum

                # Time-varying matrices so need to calculate in loop
                for i in (t+2):permanent_t
                    T_accum = TTTs[i] * T_accum
                    total_Tsum .+= T_accum

                    C_accum = TTTs[i] * C_accum + CCCs[i]
                    total_Csum .+= C_accum
                end

                # Post permanent_t, fixed matrices, so can compute here
                T_sums = (I - TTTs[permanent_t]) \ (TTTs[permanent_t] - TTTs[permanent_t]^(t+k-permanent_t+1))
                total_Tsum .+= T_sums * T_accum

                # If C unchanging after permanent_t, can compute with sum expression
                ## else do for loop
                if !do_C
                    total_Csum .+= T_sums * C_accum
                else
                    for i in (permanent_t+1):(t+k)
                        C_accum = TTTs[permanent_t] * C_accum + CCCs[permanent_t]
                        total_Csum .+= C_accum
                    end
                end
            else
                T_accum = Vector{eltype(TTTs)}(undef, k)
                C_accum = Vector{eltype(CCCs)}(undef, k)
                for i in 1:k
                    ith_memo = ForwardExpectationsMemo(memo.time_varying_memo[min(t + i, permanent_t)],
                                                                                memo.permanent_memo) # see comments above for why this is correct
                    T_accum[i], C_accum[i] = k_periods_ahead_expectations(TTT, CCC, TTTs, CCCs, t, i, permanent_t,
                                                                          integ_series = integ_series, memo = ith_memo)
                end
                total_Tsum = sum(T_accum)
                total_Csum = sum(C_accum)
                #=for i in 1:permanent_t # Using memo
                    #=ith_memo = isnothing(memo) ? nothing : ForwardExpectationsMemo(memo.time_varying_memo[min(t + i, permanent_t)],
                                                                                memo.permanent_memo) # see comments above for why this is correct
                    T_accum[i], C_accum[i] = k_periods_ahead_expectations(TTT, CCC, TTTs, CCCs, t, i, permanent_t,
                                                                          integ_series = integ_series, memo = ith_memo)=#
                    total_Tsum .+= i == permanent_t ? memo.time_varying_memo[t+1] : memo.time_varying_memo[1] / memo.time_varying_memo[t+i]
                    total_Csum .+= (I + memo.time_varying_memo[t+i]) * CCCs[t+i]
                end
                for i in (permanent_t+1):k
                    total_Tsum .+= memo.permanent_memo[i-permanent_t] * memo.time_varying_memo[1]
                    total_Csum .+= i == k ? CCCs[t+k] : (I + memo.permanent_memo[k-i]) * CCCs[t+i]
                end=#
            end

            return total_Tsum, total_Csum

            # This code seems to work (match directly adding each k_periods_ahead_expectations)
            # but sometimes causes mild numerical differences,
            # so we calculate the expected sum by adding each k_periods_ahead_expectations.
            #=            h = (permanent_t - 1) - t # last time of time-variation is permanent_t - 1

            Tₜ₊ₕ₊₁_memo = Vector{eltype(TTTs)}(undef, k - h)
            Tₜ₊ₕ₊₁_memo[1] = TTTs[permanent_t] # maps j to Tₜ₊ₕ₊₁ʲ⁻ʰ for j in (h + 1):k or Tₜ₊ₕ₊₁ʲ for j in 1:(k-h)
            for j in 2:(k - h)
            Tₜ₊ₕ₊₁_memo[j] = TTTs[permanent_t] * Tₜ₊ₕ₊₁_memo[j - 1]
            end

            Tₜ₊ₕ₊₁ʲsum = (I - TTTs[permanent_t]) \ (TTTs[permanent_t] - TTTs[permanent_t] ^ (k - h + 1)) # ∑ⱼ₌₁ᵏ⁻ʰ Tₜ₊ₕ₊₁ʲ

            TC_memo = Vector{eltype(TTTs)}(undef, h) # matrices to be used for both accumulated TTT and CCC
            TC_memo[h] = TTTs[t + h] # maps i to ∏ⱼ₌ᵢʰ Tₜ₊ⱼ, so T_memo[h] = Tₜ₊ₕ, T_memo[h-1] = Tₜ₊ₕ * Tₜ₊ₕ₋₁, ...
            for i in (h-1):-1:1
            TC_memo[i] = TC_memo[i + 1] * TTTs[t + i]
            end

            T1_term = Vector{eltype(TTTs)}(undef, h) # first part of the accumulated TTT matrix
            T1_term[1] = TTTs[t + 1] # maps i to ∏ⱼ₌₁ⁱ Tₜ₊ⱼ, so T_memo[h] = Tₜ₊ₕ * ⋯  * Tₜ₊₁
            for i in 2:(h - 1)
            T1_term[i] = TTTs[t + i] * T1_term[i - 1]
            end
            T1_term[h] = TC_memo[1] # This one was calculated already

            # second part of the accumulated TTT matrix
            # maps j to Tₜ₊ₕ₊₁ʲ ∏ᵣ₌₁ʰ Tₜ₊ᵣ, so T_memo[h] = Tₜ₊ₕ * ⋯  * Tₜ₊₁
            T2_term = [Tₜ₊ₕ₊₁_memo[j] * T1_term[h] for j in 1:(k - h)]

            T_accum = sum(T1_term) + sum(T2_term) # calculated accumulated matrix

            # Calculate final portions of accumulated CCC vector first
            I₊Tʲ = (I + Tₜ₊ₕ₊₁ʲsum)
            C_accum   = I₊Tʲ * CCCs[t + h]
            if any(.!(CCCs[permanent_t] .≈ 0.))
            C_accum .+= (I + (k - h - 1) .* I₊Tʲ) * CCCs[permanent_t]
            end
            if h > 1
            C1_term = Vector{eltype(TTTs)}(undef, h - 1)
            C2_term = Vector{eltype(TTTs)}(undef, h - 1)

            for q in 1:(h - 1)
            C1_term[q] = sum([prod([TTTs[t + m] for m in (q + 1):j]) for j in (q + 1):h])
            C2_term[q] = sum([Tₜ₊ₕ₊₁_memo[j] * TC_memo[q + 1] for j in 1:(k - h)])
            end
            C_accum .+= sum([(I + C1_term[q] + C2_term[q]) * CCCs[t + q] for q in 1:(h - 1)])
            end=#
        end
    end
end


text:
module rbc_2
using SparseArrays, LinearAlgebra, DifferentiableStateSpaceModels, LaTeXStrings, ModelingToolkit
const n_y = 3
const n_x = 2
const n = 5
const n_p = 3
const n_ϵ = 1
const n_z = 2
const functions_type() = DenseFunctions()
const η = reshape([0; -1], 2, 1)
const Q = [1.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0]
const Γ! = (var"##MTIIPVar#557", var"##MTKArg#553", var"##MTKArg#554", var"##MTKArg#555")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#553"[1], var"##MTKArg#553"[2], var"##MTKArg#553"[3], var"##MTKArg#554"[1], var"##MTKArg#554"[2], var"##MTKArg#554"[3])
                    var"##MTIIPVar#557"[1] = σ
                end
            end
        nothing
    end
const Γ_p! = (var"##MTIIPVar#564", var"##MTKArg#560", var"##MTKArg#561", var"##MTKArg#562")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#560"[1], var"##MTKArg#560"[2], var"##MTKArg#560"[3], var"##MTKArg#561"[1], var"##MTKArg#561"[2], var"##MTKArg#561"[3])
                end
            end
        nothing
    end
const Ω! = (var"##MTIIPVar#571", var"##MTKArg#567", var"##MTKArg#568", var"##MTKArg#569")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#567"[1], var"##MTKArg#567"[2], var"##MTKArg#567"[3], var"##MTKArg#568"[1], var"##MTKArg#568"[2], var"##MTKArg#568"[3])
                    var"##MTIIPVar#571"[1] = Ω_1
                    var"##MTIIPVar#571"[2] = Ω_1
                end
            end
        nothing
    end
const Ω_p! = (var"##MTIIPVar#578", var"##MTKArg#574", var"##MTKArg#575", var"##MTKArg#576")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#574"[1], var"##MTKArg#574"[2], var"##MTKArg#574"[3], var"##MTKArg#575"[1], var"##MTKArg#575"[2], var"##MTKArg#575"[3])
                end
            end
        nothing
    end
const H! = (var"##MTIIPVar#591", var"##MTKArg#581", var"##MTKArg#582", var"##MTKArg#583", var"##MTKArg#584", var"##MTKArg#585", var"##MTKArg#586", var"##MTKArg#587", var"##MTKArg#588", var"##MTKArg#589")->begin
        @inbounds begin
                let (c_p, q_p, i_p, c, q, i, c_ss, q_ss, i_ss, k_p, z_p, k, z, k_ss, z_ss, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#581"[1], var"##MTKArg#581"[2], var"##MTKArg#581"[3], var"##MTKArg#582"[1], var"##MTKArg#582"[2], var"##MTKArg#582"[3], var"##MTKArg#583"[1], var"##MTKArg#583"[2], var"##MTKArg#583"[3], var"##MTKArg#584"[1], var"##MTKArg#584"[2], var"##MTKArg#585"[1], var"##MTKArg#585"[2], var"##MTKArg#586"[1], var"##MTKArg#586"[2], var"##MTKArg#587"[1], var"##MTKArg#587"[2], var"##MTKArg#587"[3], var"##MTKArg#588"[1], var"##MTKArg#588"[2], var"##MTKArg#588"[3])
                    var"##MTIIPVar#591"[1] = (-)((/)(1, c), (*)((/)(β, c_p), (+)((*)((*)(α, (exp)(z_p)), (^)(k_p, (-)(α, 1))), (-)(1, δ))))
                    var"##MTIIPVar#591"[2] = (-)((-)((+)(c, k_p), (*)((-)(1, δ), k)), q)
                    var"##MTIIPVar#591"[3] = (-)(q, (*)((exp)(z), (^)(k, α)))
                    var"##MTIIPVar#591"[4] = (-)(z_p, (*)(ρ, z))
                    var"##MTIIPVar#591"[5] = (-)(i, (-)(k_p, (*)((-)(1, δ), k)))
                end
            end
        nothing
    end
const H_yp! = (var"##MTIIPVar#600", var"##MTKArg#594", var"##MTKArg#595", var"##MTKArg#596", var"##MTKArg#597", var"##MTKArg#598")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#594"[1], var"##MTKArg#594"[2], var"##MTKArg#594"[3], var"##MTKArg#595"[1], var"##MTKArg#595"[2], var"##MTKArg#596"[1], var"##MTKArg#596"[2], var"##MTKArg#596"[3], var"##MTKArg#597"[1], var"##MTKArg#597"[2], var"##MTKArg#597"[3])
                    var"##MTIIPVar#600"[1] = (*)((+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 2), β)
                end
            end
        nothing
    end
const H_y! = (var"##MTIIPVar#609", var"##MTKArg#603", var"##MTKArg#604", var"##MTKArg#605", var"##MTKArg#606", var"##MTKArg#607")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#603"[1], var"##MTKArg#603"[2], var"##MTKArg#603"[3], var"##MTKArg#604"[1], var"##MTKArg#604"[2], var"##MTKArg#605"[1], var"##MTKArg#605"[2], var"##MTKArg#605"[3], var"##MTKArg#606"[1], var"##MTKArg#606"[2], var"##MTKArg#606"[3])
                    var"##MTIIPVar#609"[1] = (*)(-1, (^)((inv)(c), 2))
                    var"##MTIIPVar#609"[2] = 1
                    var"##MTIIPVar#609"[7] = -1
                    var"##MTIIPVar#609"[8] = 1
                    var"##MTIIPVar#609"[15] = 1
                end
            end
        nothing
    end
const H_xp! = (var"##MTIIPVar#618", var"##MTKArg#612", var"##MTKArg#613", var"##MTKArg#614", var"##MTKArg#615", var"##MTKArg#616")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#612"[1], var"##MTKArg#612"[2], var"##MTKArg#612"[3], var"##MTKArg#613"[1], var"##MTKArg#613"[2], var"##MTKArg#614"[1], var"##MTKArg#614"[2], var"##MTKArg#614"[3], var"##MTKArg#615"[1], var"##MTKArg#615"[2], var"##MTKArg#615"[3])
                    var"##MTIIPVar#618"[1] = (*)((*)((*)((*)(-1, (exp)(z)), (^)((inv)(c), 1)), α), (*)((*)(β, (^)(k, (+)(-2, α))), (+)(-1, α)))
                    var"##MTIIPVar#618"[2] = 1
                    var"##MTIIPVar#618"[5] = -1
                    var"##MTIIPVar#618"[6] = (*)((*)((*)((*)(-1, (exp)(z)), (^)((inv)(c), 1)), α), (*)(β, (^)(k, (+)(-1, α))))
                    var"##MTIIPVar#618"[9] = 1
                end
            end
        nothing
    end
const H_x! = (var"##MTIIPVar#627", var"##MTKArg#621", var"##MTKArg#622", var"##MTKArg#623", var"##MTKArg#624", var"##MTKArg#625")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#621"[1], var"##MTKArg#621"[2], var"##MTKArg#621"[3], var"##MTKArg#622"[1], var"##MTKArg#622"[2], var"##MTKArg#623"[1], var"##MTKArg#623"[2], var"##MTKArg#623"[3], var"##MTKArg#624"[1], var"##MTKArg#624"[2], var"##MTKArg#624"[3])
                    var"##MTIIPVar#627"[2] = (*)(-1, (+)(1, (*)(-1, δ)))
                    var"##MTIIPVar#627"[3] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    var"##MTIIPVar#627"[5] = (+)(1, (*)(-1, δ))
                    var"##MTIIPVar#627"[8] = (*)(-1, (exp)(z), (^)(k, α))
                    var"##MTIIPVar#627"[9] = (*)(-1, ρ)
                end
            end
        nothing
    end
const H_yp_p! = (var"##MTIIPVar#636", var"##MTKArg#630", var"##MTKArg#631", var"##MTKArg#632", var"##MTKArg#633", var"##MTKArg#634")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#630"[1], var"##MTKArg#630"[2], var"##MTKArg#630"[3], var"##MTKArg#631"[1], var"##MTKArg#631"[2], var"##MTKArg#632"[1], var"##MTKArg#632"[2], var"##MTKArg#632"[3], var"##MTKArg#633"[1], var"##MTKArg#633"[2], var"##MTKArg#633"[3])
                    (var"##MTIIPVar#636"[1])[1] = (*)((^)((inv)(c), 2), β, (+)((*)((exp)(z), (^)(k, (+)(-1, α))), (*)((log)(k), (exp)(z), α, (^)(k, (+)(-1, α)))))
                    (var"##MTIIPVar#636"[2])[1] = (*)((+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 2))
                end
            end
        nothing
    end
const H_y_p! = (var"##MTIIPVar#645", var"##MTKArg#639", var"##MTKArg#640", var"##MTKArg#641", var"##MTKArg#642", var"##MTKArg#643")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#639"[1], var"##MTKArg#639"[2], var"##MTKArg#639"[3], var"##MTKArg#640"[1], var"##MTKArg#640"[2], var"##MTKArg#641"[1], var"##MTKArg#641"[2], var"##MTKArg#641"[3], var"##MTKArg#642"[1], var"##MTKArg#642"[2], var"##MTKArg#642"[3])
                end
            end
        nothing
    end
const H_xp_p! = (var"##MTIIPVar#654", var"##MTKArg#648", var"##MTKArg#649", var"##MTKArg#650", var"##MTKArg#651", var"##MTKArg#652")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#648"[1], var"##MTKArg#648"[2], var"##MTKArg#648"[3], var"##MTKArg#649"[1], var"##MTKArg#649"[2], var"##MTKArg#650"[1], var"##MTKArg#650"[2], var"##MTKArg#650"[3], var"##MTKArg#651"[1], var"##MTKArg#651"[2], var"##MTKArg#651"[3])
                    (var"##MTIIPVar#654"[1])[1] = (+)((*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α))), (*)(-1, (exp)(z), (^)((inv)(c), 1), β, (^)(k, (+)(-2, α)), (+)(-1, α)), (*)(-1, (log)(k), (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-2, α))))
                    (var"##MTIIPVar#654"[1])[6] = (+)((*)(-1, (exp)(z), (^)((inv)(c), 1), β, (^)(k, (+)(-1, α))), (*)(-1, (log)(k), (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-1, α))))
                    (var"##MTIIPVar#654"[2])[1] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#654"[2])[6] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (^)(k, (+)(-1, α)))
                end
            end
        nothing
    end
const H_x_p! = (var"##MTIIPVar#663", var"##MTKArg#657", var"##MTKArg#658", var"##MTKArg#659", var"##MTKArg#660", var"##MTKArg#661")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#657"[1], var"##MTKArg#657"[2], var"##MTKArg#657"[3], var"##MTKArg#658"[1], var"##MTKArg#658"[2], var"##MTKArg#659"[1], var"##MTKArg#659"[2], var"##MTKArg#659"[3], var"##MTKArg#660"[1], var"##MTKArg#660"[2], var"##MTKArg#660"[3])
                    (var"##MTIIPVar#663"[1])[3] = (+)((*)(-1, (exp)(z), (^)(k, (+)(-1, α))), (*)(-1, (log)(k), (exp)(z), α, (^)(k, (+)(-1, α))))
                    (var"##MTIIPVar#663"[1])[8] = (*)(-1, (log)(k), (exp)(z), (^)(k, α))
                    (var"##MTIIPVar#663"[3])[9] = -1
                end
            end
        nothing
    end
const H_p! = (var"##MTIIPVar#672", var"##MTKArg#666", var"##MTKArg#667", var"##MTKArg#668", var"##MTKArg#669", var"##MTKArg#670")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#666"[1], var"##MTKArg#666"[2], var"##MTKArg#666"[3], var"##MTKArg#667"[1], var"##MTKArg#667"[2], var"##MTKArg#668"[1], var"##MTKArg#668"[2], var"##MTKArg#668"[3], var"##MTKArg#669"[1], var"##MTKArg#669"[2], var"##MTKArg#669"[3])
                    var"##MTIIPVar#672"[1] = (*)(-1, (^)((inv)(c), 1), β, (+)((*)((exp)(z), (^)(k, (+)(-1, α))), (*)((log)(k), (exp)(z), α, (^)(k, (+)(-1, α)))))
                    var"##MTIIPVar#672"[3] = (*)(-1, (log)(k), (exp)(z), (^)(k, α))
                    var"##MTIIPVar#672"[6] = (*)(-1, (+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 1))
                    var"##MTIIPVar#672"[14] = (*)(-1, z)
                end
            end
        nothing
    end
const Ψ! = (var"##MTIIPVar#681", var"##MTKArg#675", var"##MTKArg#676", var"##MTKArg#677", var"##MTKArg#678", var"##MTKArg#679")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#675"[1], var"##MTKArg#675"[2], var"##MTKArg#675"[3], var"##MTKArg#676"[1], var"##MTKArg#676"[2], var"##MTKArg#677"[1], var"##MTKArg#677"[2], var"##MTKArg#677"[3], var"##MTKArg#678"[1], var"##MTKArg#678"[2], var"##MTKArg#678"[3])
                    (var"##MTIIPVar#681"[1])[1] = (*)(-2, (+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 3), β)
                    (var"##MTIIPVar#681"[1])[7] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#681"[1])[8] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#681"[1])[34] = (*)(2, (^)((inv)(c), 3))
                    (var"##MTIIPVar#681"[1])[61] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#681"[1])[67] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α))
                    (var"##MTIIPVar#681"[1])[68] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#681"[1])[71] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#681"[1])[77] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#681"[1])[78] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#681"[3])[89] = (*)(-1, (exp)(z), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    (var"##MTIIPVar#681"[3])[90] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#681"[3])[99] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    (var"##MTIIPVar#681"[3])[100] = (*)(-1, (exp)(z), (^)(k, α))
                end
            end
        nothing
    end
const Ψ_p! = (var"##MTIIPVar#690", var"##MTKArg#684", var"##MTKArg#685", var"##MTKArg#686", var"##MTKArg#687", var"##MTKArg#688")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#684"[1], var"##MTKArg#684"[2], var"##MTKArg#684"[3], var"##MTKArg#685"[1], var"##MTKArg#685"[2], var"##MTKArg#686"[1], var"##MTKArg#686"[2], var"##MTKArg#686"[3], var"##MTKArg#687"[1], var"##MTKArg#687"[2], var"##MTKArg#687"[3])
                    ((var"##MTIIPVar#690"[1])[1])[1] = (*)(-2, (^)((inv)(c), 3), β, (+)((*)((exp)(z), (^)(k, (+)(-1, α))), (*)((log)(k), (exp)(z), α, (^)(k, (+)(-1, α)))))
                    ((var"##MTIIPVar#690"[1])[1])[7] = (+)((*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α))), (*)((exp)(z), (^)((inv)(c), 2), β, (^)(k, (+)(-2, α)), (+)(-1, α)), (*)((log)(k), (exp)(z), (^)((inv)(c), 2), α, (+)(-1, α), β, (^)(k, (+)(-2, α))))
                    ((var"##MTIIPVar#690"[1])[1])[8] = (+)((*)((exp)(z), (^)((inv)(c), 2), β, (^)(k, (+)(-1, α))), (*)((log)(k), (exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α))))
                    ((var"##MTIIPVar#690"[1])[1])[61] = (+)((*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α))), (*)((exp)(z), (^)((inv)(c), 2), β, (^)(k, (+)(-2, α)), (+)(-1, α)), (*)((log)(k), (exp)(z), (^)((inv)(c), 2), α, (+)(-1, α), β, (^)(k, (+)(-2, α))))
                    ((var"##MTIIPVar#690"[1])[1])[67] = (+)((*)(-1, (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-3, α))), (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-3, α)), (+)(-2, α)), (*)(-1, (exp)(z), (^)((inv)(c), 1), (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α)), (*)(-1, (log)(k), (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), (+)(-2, α), β, (^)(k, (+)(-3, α))))
                    ((var"##MTIIPVar#690"[1])[1])[68] = (+)((*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α))), (*)(-1, (exp)(z), (^)((inv)(c), 1), β, (^)(k, (+)(-2, α)), (+)(-1, α)), (*)(-1, (log)(k), (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-2, α))))
                    ((var"##MTIIPVar#690"[1])[1])[71] = (+)((*)((exp)(z), (^)((inv)(c), 2), β, (^)(k, (+)(-1, α))), (*)((log)(k), (exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α))))
                    ((var"##MTIIPVar#690"[1])[1])[77] = (+)((*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α))), (*)(-1, (exp)(z), (^)((inv)(c), 1), β, (^)(k, (+)(-2, α)), (+)(-1, α)), (*)(-1, (log)(k), (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-2, α))))
                    ((var"##MTIIPVar#690"[1])[1])[78] = (+)((*)(-1, (exp)(z), (^)((inv)(c), 1), β, (^)(k, (+)(-1, α))), (*)(-1, (log)(k), (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-1, α))))
                    ((var"##MTIIPVar#690"[1])[3])[89] = (+)((*)(-1, (exp)(z), α, (^)(k, (+)(-2, α))), (*)(-1, (exp)(z), (^)(k, (+)(-2, α)), (+)(-1, α)), (*)(-1, (log)(k), (exp)(z), α, (^)(k, (+)(-2, α)), (+)(-1, α)))
                    ((var"##MTIIPVar#690"[1])[3])[90] = (+)((*)(-1, (exp)(z), (^)(k, (+)(-1, α))), (*)(-1, (log)(k), (exp)(z), α, (^)(k, (+)(-1, α))))
                    ((var"##MTIIPVar#690"[1])[3])[99] = (+)((*)(-1, (exp)(z), (^)(k, (+)(-1, α))), (*)(-1, (log)(k), (exp)(z), α, (^)(k, (+)(-1, α))))
                    ((var"##MTIIPVar#690"[1])[3])[100] = (*)(-1, (log)(k), (exp)(z), (^)(k, α))
                    ((var"##MTIIPVar#690"[2])[1])[1] = (*)(-2, (+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 3))
                    ((var"##MTIIPVar#690"[2])[1])[7] = (*)((exp)(z), (^)((inv)(c), 2), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#690"[2])[1])[8] = (*)((exp)(z), (^)((inv)(c), 2), α, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#690"[2])[1])[61] = (*)((exp)(z), (^)((inv)(c), 2), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#690"[2])[1])[67] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (^)(k, (+)(-3, α)), (+)(-1, α), (+)(-2, α))
                    ((var"##MTIIPVar#690"[2])[1])[68] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#690"[2])[1])[71] = (*)((exp)(z), (^)((inv)(c), 2), α, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#690"[2])[1])[77] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#690"[2])[1])[78] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (^)(k, (+)(-1, α)))
                end
            end
        nothing
    end
const Ψ_yp! = (var"##MTIIPVar#699", var"##MTKArg#693", var"##MTKArg#694", var"##MTKArg#695", var"##MTKArg#696", var"##MTKArg#697")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#693"[1], var"##MTKArg#693"[2], var"##MTKArg#693"[3], var"##MTKArg#694"[1], var"##MTKArg#694"[2], var"##MTKArg#695"[1], var"##MTKArg#695"[2], var"##MTKArg#695"[3], var"##MTKArg#696"[1], var"##MTKArg#696"[2], var"##MTKArg#696"[3])
                    ((var"##MTIIPVar#699"[1])[1])[1] = (*)(6, (+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 4), β)
                    ((var"##MTIIPVar#699"[1])[1])[7] = (*)(-2, (exp)(z), (^)((inv)(c), 3), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#699"[1])[1])[8] = (*)(-2, (exp)(z), (^)((inv)(c), 3), α, β, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#699"[1])[1])[61] = (*)(-2, (exp)(z), (^)((inv)(c), 3), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#699"[1])[1])[67] = (*)((exp)(z), (^)((inv)(c), 2), α, (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α))
                    ((var"##MTIIPVar#699"[1])[1])[68] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#699"[1])[1])[71] = (*)(-2, (exp)(z), (^)((inv)(c), 3), α, β, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#699"[1])[1])[77] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#699"[1])[1])[78] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α)))
                end
            end
        nothing
    end
const Ψ_y! = (var"##MTIIPVar#708", var"##MTKArg#702", var"##MTKArg#703", var"##MTKArg#704", var"##MTKArg#705", var"##MTKArg#706")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#702"[1], var"##MTKArg#702"[2], var"##MTKArg#702"[3], var"##MTKArg#703"[1], var"##MTKArg#703"[2], var"##MTKArg#704"[1], var"##MTKArg#704"[2], var"##MTKArg#704"[3], var"##MTKArg#705"[1], var"##MTKArg#705"[2], var"##MTKArg#705"[3])
                    ((var"##MTIIPVar#708"[1])[1])[34] = (*)(-6, (^)((inv)(c), 4))
                end
            end
        nothing
    end
const Ψ_xp! = (var"##MTIIPVar#717", var"##MTKArg#711", var"##MTKArg#712", var"##MTKArg#713", var"##MTKArg#714", var"##MTKArg#715")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#711"[1], var"##MTKArg#711"[2], var"##MTKArg#711"[3], var"##MTKArg#712"[1], var"##MTKArg#712"[2], var"##MTKArg#713"[1], var"##MTKArg#713"[2], var"##MTKArg#713"[3], var"##MTKArg#714"[1], var"##MTKArg#714"[2], var"##MTKArg#714"[3])
                    ((var"##MTIIPVar#717"[1])[1])[1] = (*)(-2, (exp)(z), (^)((inv)(c), 3), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#717"[1])[1])[7] = (*)((exp)(z), (^)((inv)(c), 2), α, (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α))
                    ((var"##MTIIPVar#717"[1])[1])[8] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#717"[1])[1])[61] = (*)((exp)(z), (^)((inv)(c), 2), α, (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α))
                    ((var"##MTIIPVar#717"[1])[1])[67] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), (+)(-2, α), β, (^)(k, (+)(-4, α)), (+)(-3, α))
                    ((var"##MTIIPVar#717"[1])[1])[68] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α))
                    ((var"##MTIIPVar#717"[1])[1])[71] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#717"[1])[1])[77] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α))
                    ((var"##MTIIPVar#717"[1])[1])[78] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#717"[2])[1])[1] = (*)(-2, (exp)(z), (^)((inv)(c), 3), α, β, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#717"[2])[1])[7] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#717"[2])[1])[8] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#717"[2])[1])[61] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#717"[2])[1])[67] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, (+)(-1, α), β, (^)(k, (+)(-3, α)), (+)(-2, α))
                    ((var"##MTIIPVar#717"[2])[1])[68] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#717"[2])[1])[71] = (*)((exp)(z), (^)((inv)(c), 2), α, β, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#717"[2])[1])[77] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#717"[2])[1])[78] = (*)(-1, (exp)(z), (^)((inv)(c), 1), α, β, (^)(k, (+)(-1, α)))
                end
            end
        nothing
    end
const Ψ_x! = (var"##MTIIPVar#726", var"##MTKArg#720", var"##MTKArg#721", var"##MTKArg#722", var"##MTKArg#723", var"##MTKArg#724")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#720"[1], var"##MTKArg#720"[2], var"##MTKArg#720"[3], var"##MTKArg#721"[1], var"##MTKArg#721"[2], var"##MTKArg#722"[1], var"##MTKArg#722"[2], var"##MTKArg#722"[3], var"##MTKArg#723"[1], var"##MTKArg#723"[2], var"##MTKArg#723"[3])
                    ((var"##MTIIPVar#726"[1])[3])[89] = (*)(-1, (exp)(z), α, (^)(k, (+)(-3, α)), (+)(-1, α), (+)(-2, α))
                    ((var"##MTIIPVar#726"[1])[3])[90] = (*)(-1, (exp)(z), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#726"[1])[3])[99] = (*)(-1, (exp)(z), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#726"[1])[3])[100] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#726"[2])[3])[89] = (*)(-1, (exp)(z), α, (^)(k, (+)(-2, α)), (+)(-1, α))
                    ((var"##MTIIPVar#726"[2])[3])[90] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#726"[2])[3])[99] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    ((var"##MTIIPVar#726"[2])[3])[100] = (*)(-1, (exp)(z), (^)(k, α))
                end
            end
        nothing
    end
const H̄! = (var"##MTIIPVar#734", var"##MTKArg#729", var"##MTKArg#730", var"##MTKArg#731", var"##MTKArg#732")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#729"[1], var"##MTKArg#729"[2], var"##MTKArg#729"[3], var"##MTKArg#729"[4], var"##MTKArg#729"[5], var"##MTKArg#730"[1], var"##MTKArg#730"[2], var"##MTKArg#730"[3], var"##MTKArg#731"[1], var"##MTKArg#731"[2], var"##MTKArg#731"[3])
                    var"##MTIIPVar#734"[1] = (+)((^)((inv)(c), 1), (*)(-1, (+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 1), β))
                    var"##MTIIPVar#734"[2] = (+)(c, k, (*)(-1, (+)(q, (*)(k, (+)(1, (*)(-1, δ))))))
                    var"##MTIIPVar#734"[3] = (+)((*)(-1, (exp)(z), (^)(k, α)), q)
                    var"##MTIIPVar#734"[4] = (+)((*)(-1, z, ρ), z)
                    var"##MTIIPVar#734"[5] = (+)(i, (*)(-1, (+)((*)(-1, k, (+)(1, (*)(-1, δ))), k)))
                end
            end
        nothing
    end
const H̄_w! = (var"##MTIIPVar#742", var"##MTKArg#737", var"##MTKArg#738", var"##MTKArg#739", var"##MTKArg#740")->begin
        @inbounds begin
                let (c, q, i, k, z, α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#737"[1], var"##MTKArg#737"[2], var"##MTKArg#737"[3], var"##MTKArg#737"[4], var"##MTKArg#737"[5], var"##MTKArg#738"[1], var"##MTKArg#738"[2], var"##MTKArg#738"[3], var"##MTKArg#739"[1], var"##MTKArg#739"[2], var"##MTKArg#739"[3])
                    var"##MTIIPVar#742"[1] = (+)((*)(-1, (^)((inv)(c), 2)), (*)((+)(1, (*)(-1, δ), (*)((exp)(z), α, (^)(k, (+)(-1, α)))), (^)((inv)(c), 2), β))
                    var"##MTIIPVar#742"[2] = 1
                    var"##MTIIPVar#742"[7] = -1
                    var"##MTIIPVar#742"[8] = 1
                    var"##MTIIPVar#742"[15] = 1
                    var"##MTIIPVar#742"[16] = (*)((*)((*)((*)(-1, (exp)(z)), (^)((inv)(c), 1)), α), (*)((*)(β, (^)(k, (+)(-2, α))), (+)(-1, α)))
                    var"##MTIIPVar#742"[17] = (+)(1, (*)(-1, (+)(1, (*)(-1, δ))))
                    var"##MTIIPVar#742"[18] = (*)(-1, (exp)(z), α, (^)(k, (+)(-1, α)))
                    var"##MTIIPVar#742"[20] = (*)(-1, (+)(1, (*)(-1, (+)(1, (*)(-1, δ)))))
                    var"##MTIIPVar#742"[21] = (*)((*)((*)((*)(-1, (exp)(z)), (^)((inv)(c), 1)), α), (*)(β, (^)(k, (+)(-1, α))))
                    var"##MTIIPVar#742"[23] = (*)(-1, (exp)(z), (^)(k, α))
                    var"##MTIIPVar#742"[24] = (+)(1, (*)(-1, ρ))
                end
            end
        nothing
    end
const ȳ_iv! = (var"##MTIIPVar#749", var"##MTKArg#745", var"##MTKArg#746", var"##MTKArg#747")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#745"[1], var"##MTKArg#745"[2], var"##MTKArg#745"[3], var"##MTKArg#746"[1], var"##MTKArg#746"[2], var"##MTKArg#746"[3])
                    var"##MTIIPVar#749"[1] = (-)((^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(α, (-)(α, 1))), (*)(δ, (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1)))))
                    var"##MTIIPVar#749"[2] = (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(α, (-)(α, 1)))
                    var"##MTIIPVar#749"[3] = (*)(δ, (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1))))
                end
            end
        nothing
    end
const x̄_iv! = (var"##MTIIPVar#756", var"##MTKArg#752", var"##MTKArg#753", var"##MTKArg#754")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#752"[1], var"##MTKArg#752"[2], var"##MTKArg#752"[3], var"##MTKArg#753"[1], var"##MTKArg#753"[2], var"##MTKArg#753"[3])
                    var"##MTIIPVar#756"[1] = (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1)))
                end
            end
        nothing
    end
const ȳ! = (var"##MTIIPVar#763", var"##MTKArg#759", var"##MTKArg#760", var"##MTKArg#761")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#759"[1], var"##MTKArg#759"[2], var"##MTKArg#759"[3], var"##MTKArg#760"[1], var"##MTKArg#760"[2], var"##MTKArg#760"[3])
                    var"##MTIIPVar#763"[1] = (-)((^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(α, (-)(α, 1))), (*)(δ, (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1)))))
                    var"##MTIIPVar#763"[2] = (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(α, (-)(α, 1)))
                    var"##MTIIPVar#763"[3] = (*)(δ, (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1))))
                end
            end
        nothing
    end
const x̄! = (var"##MTIIPVar#770", var"##MTKArg#766", var"##MTKArg#767", var"##MTKArg#768")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#766"[1], var"##MTKArg#766"[2], var"##MTKArg#766"[3], var"##MTKArg#767"[1], var"##MTKArg#767"[2], var"##MTKArg#767"[3])
                    var"##MTIIPVar#770"[1] = (^)((/)((+)((-)((/)(1, β), 1), δ), α), (/)(1, (-)(α, 1)))
                end
            end
        nothing
    end
const ȳ_p! = (var"##MTIIPVar#777", var"##MTKArg#773", var"##MTKArg#774", var"##MTKArg#775")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#773"[1], var"##MTKArg#773"[2], var"##MTKArg#773"[3], var"##MTKArg#774"[1], var"##MTKArg#774"[2], var"##MTKArg#774"[3])
                    var"##MTIIPVar#777"[1] = (+)((*)(-1, δ, (+)((*)(-1, (log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (^)((inv)((+)(-1, α)), 1)), (^)((inv)((+)(-1, α)), 2)), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 2)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1)))), (*)((log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (+)((^)((inv)((+)(-1, α)), 1), (*)(-1, α, (^)((inv)((+)(-1, α)), 2))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (*)(α, (^)((inv)((+)(-1, α)), 1)))), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 1)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (*)(α, (^)((inv)((+)(-1, α)), 1))))), (^)((inv)((+)(-1, α)), 1)))
                    var"##MTIIPVar#777"[2] = (+)((*)((log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (+)((^)((inv)((+)(-1, α)), 1), (*)(-1, α, (^)((inv)((+)(-1, α)), 2))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (*)(α, (^)((inv)((+)(-1, α)), 1)))), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 1)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (*)(α, (^)((inv)((+)(-1, α)), 1))))), (^)((inv)((+)(-1, α)), 1)))
                    var"##MTIIPVar#777"[3] = (*)(δ, (+)((*)(-1, (log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (^)((inv)((+)(-1, α)), 1)), (^)((inv)((+)(-1, α)), 2)), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 2)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1))))
                    var"##MTIIPVar#777"[4] = (+)((*)(-1, (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (*)(α, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1), (^)((inv)(β), 2)), (*)((*)((*)((*)((^)((inv)(α), 1), δ), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1)), (^)((inv)(β), 2)))
                    var"##MTIIPVar#777"[5] = (*)(-1, (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (*)(α, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1), (^)((inv)(β), 2))
                    var"##MTIIPVar#777"[6] = (*)((*)((*)((*)(-1, (^)((inv)(α), 1)), δ), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (*)((^)((inv)((+)(-1, α)), 1), (^)((inv)(β), 2)))
                end
            end
        nothing
    end
const x̄_p! = (var"##MTIIPVar#784", var"##MTKArg#780", var"##MTKArg#781", var"##MTKArg#782")->begin
        @inbounds begin
                let (α, β, ρ, δ, σ, Ω_1) = (var"##MTKArg#780"[1], var"##MTKArg#780"[2], var"##MTKArg#780"[3], var"##MTKArg#781"[1], var"##MTKArg#781"[2], var"##MTKArg#781"[3])
                    var"##MTIIPVar#784"[1] = (+)((*)(-1, (log)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1))), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (^)((inv)((+)(-1, α)), 1)), (^)((inv)((+)(-1, α)), 2)), (*)((*)((*)((*)(-1, (+)(-1, (^)((inv)(β), 1), δ)), (^)((inv)(α), 2)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1)))
                    var"##MTIIPVar#784"[3] = (*)((*)((*)((*)(-1, (^)((inv)(α), 1)), (^)((*)((+)(-1, (^)((inv)(β), 1), δ), (^)((inv)(α), 1)), (+)(-1, (^)((inv)((+)(-1, α)), 1)))), (^)((inv)((+)(-1, α)), 1)), (^)((inv)(β), 2))
                end
            end
        nothing
    end
const steady_state! = nothing
const allocate_solver_cache = m->begin
        #= C:\Google Drive\Projects\HMCDSGE\src\mtk_constructor.jl:512 =#
        SecondOrderSolverCache(m)
    end
const select_p_ss_hash = LinearAlgebra.UniformScaling{Bool}(true)
const select_p_f_ss_hash = LinearAlgebra.UniformScaling{Bool}(true)
const select_p_perturbation_hash = LinearAlgebra.UniformScaling{Bool}(true)
const select_p_f_perturbation_hash = LinearAlgebra.UniformScaling{Bool}(true)
end

ext: jl


text:
function F1d(F::Vector{Float64},Q::Float64,A::Float64,beta::Float64,rho::Float64,ts::Float64,
    xs::Float64,Qs::Float64,As::Float64)
    F[1] = ts*Qs./(As*xs)*Q;
    F[2] = Qs*ts/(As*xs)*(Q^2)/A + ts*(As^1.5)./(Qs*xs)*beta/(3*rho)*A^1.5;
    # println("A flux: $(ts*Qs./(As*xs)*Q)")
    # println("Q flux, term 1: $(Qs*ts/(As*xs)*Q^2/A)")
    # println("Q flux, term 2: $(ts*As^1.5./(Qs*xs)*beta/(3*rho)*A^1.5)")
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
#############################
## Degree
#############################
==(d1::Degree, d2::Degree) = (d1.explvar==d2.explvar) && (d1.conjvar==d2.conjvar)
!=(d1::Degree, d2::Degree) = !(d1 == d2)

hash(deg::Degree, h::UInt) = hash(deg.explvar, hash(deg.conjvar, h))

function isless(deg1::Degree, deg2::Degree)
  return (deg1.explvar<deg2.explvar) || (deg1.explvar==deg2.explvar && deg1.conjvar<deg2.conjvar)
end

#############################
## Variable
#############################
==(x::Variable, y::Variable) = (x.name==y.name) && (x.kind==y.kind)
!=(x::Variable, y::Variable) = !(x == y)
hash(var::Variable, h::UInt) = hash(var.name, hash(var.kind, h))

function isless(a::Variable, b::Variable)
  return isless(a.name, b.name)
end

#############################
## Exponent
#############################
function ==(exp1::Exponent, exp2::Exponent)
  (length(exp1.expo) == length(exp2.expo)) || return false

  for (varname, deg) in exp2.expo
    if !haskey(exp1.expo, varname) || exp1.expo[varname] != deg
      return false
    end
  end
  return true
end
!=(exp1::Exponent, exp2::Exponent) = !(exp1 == exp2)

hash(expo::Exponent, h::UInt) = hash(expo.degree, hash(expo.expo, h))

function isless(exp1::Exponent, exp2::Exponent)
  exp1_deg = exp1.degree.explvar + exp1.degree.conjvar
  exp2_deg = exp2.degree.explvar + exp2.degree.conjvar
  if exp1_deg < exp2_deg
    return true
  elseif exp1_deg == exp2_deg
    vars = union(Set(keys(exp1)), Set(keys(exp2)))
    for var in sort(collect(vars))
      if !haskey(exp2, var) # hence haskey(exp1, var)
        return true
      elseif !haskey(exp1, var) # hence haskey(exp2, var)
        return false
      elseif exp1[var] > exp2[var]
        return true
      elseif exp1[var] == exp2[var]
        continue
      else
        return false
      end
    end
  else
    return false
  end
  return false
end

#############################
## Polynomial
#############################
function ==(p1::Polynomial, p2::Polynomial)
  if (length(p1) != length(p2)) || (p1.degree != p2.degree)
    return false
  end
  for (exponent, coeff) in p1
    if !haskey(p2, exponent) || (coeff != p2.poly[exponent])
      return false
    end
  end
  true
end
!=(pol1::Polynomial, pol2::Polynomial) = !(pol1 == pol2)

# NOTE: order on polynomials ?

#############################
## Point
#############################
function ==(pt1::Point, pt2::Point)
  if length(pt1) != length(pt2)
    return false
  end
  for (var1, val1) in pt1
    if !haskey(pt2, var1) || pt2[var1] != val1
      return false
    end
  end
  true
end
!=(pt1::Point, pt2::Point) = !(pt1 == pt2)

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
export skipindex, @symmetrical

using SpecialFunctions
using Rocket

import Base: show, similar
import Base: IteratorSize, HasLength
import Base: IteratorEltype, HasEltype
import Base: eltype, length, size, sum
import Base: IndexStyle, IndexLinear, getindex

import Rocket: similar_typeof

struct SkipIndexIterator{T, I} <: AbstractVector{T}
    iterator :: I
    skip     :: Int
end

skip(iter::SkipIndexIterator) = iter.skip

function skipindex(iterator::I, skip::Int) where I
    @assert skip >= 1
    @assert length(iterator) >= 1
    return SkipIndexIterator{eltype(I), I}(iterator, skip)
end

Base.IteratorSize(::Type{<:SkipIndexIterator})   = HasLength()
Base.IteratorEltype(::Type{<:SkipIndexIterator}) = HasEltype()
Base.IndexStyle(::Type{<:SkipIndexIterator})     = IndexLinear()

Base.eltype(::Type{<:SkipIndexIterator{T}}) where T = T
Base.length(iter::SkipIndexIterator)                = length(iter.iterator) - 1
Base.size(iter::SkipIndexIterator)                  = (length(iter), )

Base.getindex(iter::SkipIndexIterator, i::Int)               = i < skip(iter) ? @inbounds(iter.iterator[i]) : @inbounds(iter.iterator[i + 1])
Base.getindex(iter::SkipIndexIterator, i::CartesianIndex{1}) = Base.getindex(iter, first(i.I))

Rocket.similar_typeof(::SkipIndexIterator, ::Type{L}) where L = Vector{L}

cast_to_message_subscribable(some::T) where T = cast_to_message_subscribable(as_subscribable(T), some)

cast_to_message_subscribable(::InvalidSubscribableTrait, some)   = of(as_message(some))
cast_to_message_subscribable(::SimpleSubscribableTrait, some)    = some |> map(Message, as_message)
cast_to_message_subscribable(::ScheduledSubscribableTrait, some) = some |> map(Message, as_message)

reduce_with_sum(array) = reduce(+, array)

##

"""
    OneDivNVector(N::Int)
    OneDivNVector(::Type{T}, N::Int) where T

Allocation-free version of `fill(one(T) / N, N)` vector.

# Arguments
- `::Type{T}`: type of elements, optional, Float64 by default, should be a subtype of `Number`
- `N::Int`: number of elements in a container, should be greater than zero

# Examples

```jldoctest
julia> iter = ReactiveMP.OneDivNVector(3)
OneDivNVector(Float64, 3)

julia> length(iter)
3

julia> eltype(iter)
Float64

julia> collect(iter)
3-element Vector{Float64}:
 0.3333333333333333
 0.3333333333333333
 0.3333333333333333

julia> iter = ReactiveMP.OneDivNVector(Float32, 3)
OneDivNVector(Float32, 3)

julia> collect(iter)
3-element Vector{Float32}:
 0.33333334
 0.33333334
 0.33333334
```

See also: [`SampleList`](@ref)
"""
struct OneDivNVector{N, T} <: AbstractVector{T} end

Base.show(io::IO, ::OneDivNVector{N, T}) where { N, T } = print(io, "OneDivNVector($T, $N)")

function OneDivNVector(N::Int)
    return OneDivNVector(Float64, N)
end

function OneDivNVector(::Type{T}, N::Int) where T
    @assert N > 0 "OneDivNVector creation error: N should be greater than zero"
    @assert T <: Number "OneDivNVector creation error: T should be a subtype of `Number`"
    return OneDivNVector{N, T}()
end

Base.IteratorSize(::Type{ <: OneDivNVector })   = Base.HasLength()
Base.IteratorEltype(::Type{ <: OneDivNVector }) = Base.HasEltype()

Base.sum(::OneDivNVector{ N, T }) where { N, T } = one(T)
Base.eltype(::Type{ <: OneDivNVector{N, T} }) where { N, T } = T
Base.length(::OneDivNVector{N})               where N        = N
Base.size(::OneDivNVector{N})                 where N        = (N, )

Base.iterate(::OneDivNVector{N, T})        where { N, T } = (one(T) / N, 1)
Base.iterate(::OneDivNVector{N, T}, state) where { N, T } = state >= N ? nothing : (one(T) / N, state + 1)

Base.getindex(v::OneDivNVector{N, T}, index::Int) where { N, T } = 1 <= index <= N ? (one(T) / N) : throw(BoundsError(v, index))

Base.similar(v::OneDivNVector) = v
Base.similar(v::OneDivNVector{N}, ::Type{ T }) where { N, T } = OneDivNVector(T, N)

Base.vec(::OneDivNVector{N, T}) where { N, T } = fill(one(T) / N, N)

##

import Base: +, -, *, /, convert, float, isfinite, isinf, zero, eltype

struct Infinity
    degree :: Int
end

degree(a::Infinity) = a.degree

const ∞ = Infinity(1)

Base.:+(a::Infinity, b::Infinity) = Infinity(degree(a) + degree(b))
Base.:-(a::Infinity, b::Infinity) = Infinity(degree(a) - degree(b))
Base.:+(a::Infinity)              = Infinity(+degree(a))
Base.:-(a::Infinity)              = Infinity(-degree(a))

Base.:*(i::Infinity, d::Int) = Infinity(degree(i) * d)
Base.:*(d::Int, i::Infinity) = Infinity(degree(i) * d)

Base.:*(::Infinity, ::Real)     = error("Infinity multiplication on real numbers is disallowed")
Base.:*(::Real, ::Infinity)     = error("Infinity multiplication on real numbers is disallowed")
Base.:*(::Infinity, ::Infinity) = error("Infinity multiplication is disallowed")
Base.:/(::Infinity, ::Real)     = error("Infinity division on real numbers is disallowed")
Base.:/(::Real, ::Infinity)     = error("Infinity division on real numbers is disallowed")
Base.:/(::Infinity, ::Infinity) = error("Infinity division is disallowed")

Base.zero(::Type{Infinity})       = Infinity(0)

Base.show(io::IO, a::Infinity) = print(io, "Infinity($(degree(a)))")

struct InfCountingReal{ T <: Real }
    value :: T
    infs  :: Int
end

InfCountingReal(value::T)                 where { T <: Real } = InfCountingReal{T}(value, 0)
InfCountingReal(::Type{T}, inf::Infinity) where { T <: Real } = InfCountingReal{T}(zero(T), degree(inf))

value(a::InfCountingReal) = a.value
infs(a::InfCountingReal)  = a.infs

Base.isfinite(a::InfCountingReal) = infs(a) === 0
Base.isinf(a::InfCountingReal)    = !(isfinite(a))

Base.eltype(::Type{ InfCountingReal{T} }) where T = T
Base.eltype(::Type{ InfCountingReal })            = Real

Base.eltype(::T) where { T <: InfCountingReal } = eltype(T)

Base.:+(a::Infinity, b::Real) = InfCountingReal(b, degree(a))
Base.:-(a::Infinity, b::Real) = InfCountingReal(-b, degree(a))
Base.:+(b::Real, a::Infinity) = InfCountingReal(b, +degree(a))
Base.:-(b::Real, a::Infinity) = InfCountingReal(b, -degree(a))

Base.:+(a::InfCountingReal) = InfCountingReal(+value(a), +infs(a))
Base.:-(a::InfCountingReal) = InfCountingReal(-value(a), -infs(a))

Base.:+(a::InfCountingReal, b::Infinity) = InfCountingReal(value(a), infs(a) + degree(b))
Base.:-(a::InfCountingReal, b::Infinity) = InfCountingReal(value(a), infs(a) - degree(b))

Base.:+(b::Infinity, a::InfCountingReal) = InfCountingReal(+value(a), degree(b) + infs(a))
Base.:-(b::Infinity, a::InfCountingReal) = InfCountingReal(-value(a), degree(b) - infs(a))

Base.:+(a::InfCountingReal, b::Real) = InfCountingReal(value(a) + b, infs(a))
Base.:-(a::InfCountingReal, b::Real) = InfCountingReal(value(a) - b, infs(a))
Base.:+(b::Real, a::InfCountingReal) = InfCountingReal(b + value(a), +infs(a))
Base.:-(b::Real, a::InfCountingReal) = InfCountingReal(b - value(a), -infs(a))

Base.:*(::InfCountingReal, ::Real) = error("InfCountingReal multiplication is dissalowed")
Base.:/(::InfCountingReal, ::Real) = error("InfCountingReal division is dissalowed")
Base.:*(::Real, ::InfCountingReal) = error("InfCountingReal multiplication is dissalowed")
Base.:/(::Real, ::InfCountingReal) = error("InfCountingReal division is dissalowed")

Base.:+(a::InfCountingReal, b::InfCountingReal) = InfCountingReal(value(a) + value(b), infs(a) + infs(b))
Base.:-(a::InfCountingReal, b::InfCountingReal) = InfCountingReal(value(a) - value(b), infs(a) - infs(b))

Base.convert(::Type{ InfCountingReal },    v::T)                  where { T <: Real }            = InfCountingReal(v)
Base.convert(::Type{ InfCountingReal{T} }, v::T)                  where { T <: Real }            = InfCountingReal(v)
Base.convert(::Type{ InfCountingReal{T} }, v::R)                  where { T <: Real, R <: Real } = InfCountingReal(convert(T, v))
Base.convert(::Type{ InfCountingReal{T} }, inf::Infinity)         where { T <: Real }            = InfCountingReal(T, inf)
Base.convert(::Type{ InfCountingReal{T} }, v::InfCountingReal{R}) where { T <: Real, R <: Real } = InfCountingReal{T}(convert(T, value(v)), infs(v))

Base.float(a::InfCountingReal) = isfinite(a) ? value(a) : Inf

Base.zero(::Type{InfCountingReal{T}}) where { T <: Real } = InfCountingReal(zero(T))

Base.show(io::IO, a::InfCountingReal{T}) where T = print(io, "InfCountingReal($(value(a)), $(infs(a))∞)")

# Union helpers

union_types(x::Union) = (x.a, union_types(x.b)...)
union_types(x::Type)  = (x,)

# Symbol helpers

__extract_val_type(::Type{ Type{ Val{ S } } }) where S = S

@generated function split_underscored_symbol(symbol_val)
    S = __extract_val_type(symbol_val)
    R = tuple(map(Symbol, split(string(S), "_"))...)
    return :(Val{ $R })
end

# NamedTuple helpers

fields(::NamedTuple{ F }) where F  = F
hasfield(field::Symbol, ntuple::NamedTuple) = field ∈ fields(ntuple)

function swapped(tuple::Tuple, i, j)
    @assert j > i
    return (tuple[1:i-1]..., tuple[j], tuple[i+1:j-1]..., tuple[i], tuple[j+1:end]...)
end

function swapped(array::AbstractArray, i, j)
    array = copy(array)
    array[i], array[j] = array[j], array[i]
    return array
end

##

__check_all(fn::Function, iterator)  = all(fn, iterator)
__check_all(fn::Function, ::Nothing) = true

## Meta utils

"""
    @symmetrical `function_definition`
Duplicate a method definition with the order of the first two arguments swapped.
This macro is used to duplicate methods that are symmetrical in their first two input arguments,
but require explicit definitions for the different argument orders.
Example:
    @symmetrical function prod!(x, y, z)
        ...
    end
"""
macro symmetrical(fn::Expr)
    # Check if macro is applied to a function definition
    # Valid function definitions include:
    # 1. foo([ args... ]) [ where ... [ where ... [ ... ] ] ] = :block
    # 2. function foo([ args... ]) [ where ... [ where ... [ ... ] ] ]
    #        :block
    #    end
    if (fn.head === :(=) || fn.head === :function) &&
        (fn.args[1] isa Expr && fn.args[2] isa Expr) &&
        (fn.args[2].head === :block)
        return esc(quote
            $fn
            $(swap_arguments(fn))
        end)
    else
        error("@symmetrical macro can be applied only to function definitions")
    end
end

function swap_arguments(fn::Expr)
    swapped = copy(fn)

    if swapped.args[1].head === :where
        swapped.args[1] = swap_arguments(swapped.args[1])
    elseif swapped.args[1].head === :call && length(fn.args[1].args) >= 3 # Note: >= 3, because the first argument is a function name
        swapped.args[1].args[2] = fn.args[1].args[3]
        swapped.args[1].args[3] = fn.args[1].args[2]
    else
        error("Function method passed for @symmetrical macro must have more than 2 arguments")
    end

    return swapped
end

## Other helpers

# We override this function for some specific types
function is_typeof_equal(left, right)
    _isequal = typeof(left) === typeof(right)
    if !_isequal
        @warn "typeof($left) !== typeof($right)"
    end
    return _isequal
end

function custom_isapprox(left, right; kwargs...)
    _isapprox = isapprox(left, right; kwargs...)
    if !_isapprox
        @warn "$left !≈ $right"
    end
    return _isapprox
end

custom_isapprox(left::NamedTuple, right::NamedTuple; kwargs...) = false

function custom_isapprox(left::NamedTuple{K}, right::NamedTuple{K}; kwargs...) where { K }
    _isapprox = true
    for key in keys(left)
        _isapprox = _isapprox && custom_isapprox(left[key], right[key]; kwargs...)
    end
    if !_isapprox
        @warn "$left !≈ $right"
    end
    return _isapprox
end

##

deep_eltype(::Type{ T })  where { T <: Number } = T
deep_eltype(::Type{ T })  where T               = deep_eltype(eltype(T))
deep_eltype(::T)          where T               = deep_eltype(T)


text:
using KitBase, FluxReconstruction, LinearAlgebra, OrdinaryDiffEq, Plots
using KitBase.OffsetArrays
using KitBase.ProgressMeter: @showprogress
using Base.Threads: @threads

begin
    set = Setup(
        case = "shockvortex",
        space = "2d0f0v",
        flux = "hll",
        collision = "nothing",
        interpOrder = 3,
        limiter = "positivity",
        boundary = "fix",
        cfl = 0.1,
        maxTime = 1.0,
    )
    ps = FRPSpace2D(0.0, 2.0, 100, 0.0, 1.0, 50, set.interpOrder - 1, 1, 1)
    vs = nothing
    gas = Gas(Kn = 1e-6, Ma = 1.12, K = 1.0)
    ib = nothing

    ks = SolverSet(set, ps, vs, gas, ib)
end

function dudt!(du, u, p, t)
    du .= 0.0

    f,
    u_face,
    f_face,
    fx_interaction,
    fy_interaction,
    rhs1,
    rhs2,
    iJ,
    ll,
    lr,
    dhl,
    dhr,
    lpdm,
    γ = p

    nx = size(u, 1) - 2
    ny = size(u, 2) - 2
    nsp = size(u, 3)

    @inbounds @threads for l = 1:nsp
        for k = 1:nsp, j in axes(f, 2), i in axes(f, 1)
            fg, gg = euler_flux(u[i, j, k, l, :], γ)
            for s = 1:4
                f[i, j, k, l, s, :] .= iJ[i, j][k, l] * [fg[s], gg[s]]
            end
        end
    end

    @inbounds @threads for m = 1:4
        for l = 1:nsp, j in axes(u_face, 2), i in axes(u_face, 1)
            u_face[i, j, 1, l, m] = dot(u[i, j, l, :, m], ll)
            u_face[i, j, 2, l, m] = dot(u[i, j, :, l, m], lr)
            u_face[i, j, 3, l, m] = dot(u[i, j, l, :, m], lr)
            u_face[i, j, 4, l, m] = dot(u[i, j, :, l, m], ll)

            for n = 1:2
                f_face[i, j, 1, l, m, n] = dot(f[i, j, l, :, m, n], ll)
                f_face[i, j, 2, l, m, n] = dot(f[i, j, :, l, m, n], lr)
                f_face[i, j, 3, l, m, n] = dot(f[i, j, l, :, m, n], lr)
                f_face[i, j, 4, l, m, n] = dot(f[i, j, :, l, m, n], ll)
            end
        end
    end

    @inbounds @threads for k = 1:nsp
        for j = 1:ny, i = 1:nx+1
            fw = @view fx_interaction[i, j, k, :]
            uL = @view u_face[i-1, j, 2, k, :]
            uR = @view u_face[i, j, 4, k, :]
            flux_hll!(fw, uL, uR, γ, 1.0)
        end
    end
    @inbounds @threads for k = 1:nsp
        for j = 1:ny+1, i = 1:nx
            fw = @view fy_interaction[i, j, k, :]
            uL = local_frame(u_face[i, j-1, 3, k, :], 0.0, 1.0)
            uR = local_frame(u_face[i, j, 1, k, :], 0.0, 1.0)
            flux_hll!(fw, uL, uR, γ, 1.0)
            fw .= global_frame(fw, 0.0, 1.0)
        end
    end

    @inbounds @threads for m = 1:4
        for l = 1:nsp, k = 1:nsp, j = 1:ny, i = 1:nx
            rhs1[i, j, k, l, m] = dot(f[i, j, :, l, m, 1], lpdm[k, :])
        end
    end
    @inbounds @threads for m = 1:4
        for l = 1:nsp, k = 1:nsp, j = 1:ny, i = 1:nx
            rhs2[i, j, k, l, m] = dot(f[i, j, k, :, m, 2], lpdm[l, :])
        end
    end

    @inbounds @threads for m = 1:4
        for l = 1:nsp, k = 1:nsp, j = 1:ny, i = 1:nx
            du[i, j, k, l, m] = -(
                rhs1[i, j, k, l, m] +
                rhs2[i, j, k, l, m] +
                (
                    fx_interaction[i, j, l, m] * iJ[i, j][k, l][1, 1] -
                    f_face[i, j, 4, l, m, 1]
                ) * dhl[k] +
                (
                    fx_interaction[i+1, j, l, m] * iJ[i, j][k, l][1, 1] -
                    f_face[i, j, 2, l, m, 1]
                ) * dhr[k] +
                (
                    fy_interaction[i, j, k, m] * iJ[i, j][k, l][2, 2] -
                    f_face[i, j, 1, k, m, 2]
                ) * dhl[l] +
                (
                    fy_interaction[i, j+1, k, m] * iJ[i, j][k, l][2, 2] -
                    f_face[i, j, 3, k, m, 2]
                ) * dhr[l]
            )
        end
    end

    return nothing
end

begin
    f = OffsetArray{Float64}(
        undef,
        0:ks.ps.nx+1,
        0:ks.ps.ny+1,
        ks.ps.deg + 1,
        ks.ps.deg + 1,
        4,
        2,
    )
    u_face = OffsetArray{Float64}(undef, 0:ks.ps.nx+1, 0:ks.ps.ny+1, 4, ks.ps.deg + 1, 4)
    f_face = OffsetArray{Float64}(undef, 0:ks.ps.nx+1, 0:ks.ps.ny+1, 4, ks.ps.deg + 1, 4, 2)
    fx_interaction = zeros(ks.ps.nx + 1, ks.ps.ny, ks.ps.deg + 1, 4)
    fy_interaction = zeros(ks.ps.nx, ks.ps.ny + 1, ks.ps.deg + 1, 4)
    rhs1 = zeros(ks.ps.nx, ks.ps.ny, ks.ps.deg + 1, ks.ps.deg + 1, 4)
    rhs2 = zeros(ks.ps.nx, ks.ps.ny, ks.ps.deg + 1, ks.ps.deg + 1, 4)
end

p = (
    f,
    u_face,
    f_face,
    fx_interaction,
    fy_interaction,
    rhs1,
    rhs2,
    ps.iJ,
    ps.ll,
    ps.lr,
    ps.dhl,
    ps.dhr,
    ps.dl,
    ks.gas.γ,
)

tspan = (0.0, 0.5)
dt = 0.001
nt = tspan[2] ÷ dt |> Int

# initial condition
u0 = OffsetArray{Float64}(undef, 0:ps.nx+1, 0:ps.ny+1, ps.deg + 1, ps.deg + 1, 4)
for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
    t1, t2 = ib_rh(ks.gas.Ma, ks.gas.γ)

    if ps.x[i, j] <= ps.x1 * 0.25
        prim = [t2[1], t1[2] - t2[2], 0.0, t2[3]]
    else
        prim = [t1[1], 0.0, 0.0, t1[3]]
        s = prim[1]^(1 - ks.gas.γ) / (2 * prim[end])

        κ = 0.2
        μ = 0.204
        rc = 0.05

        x0 = 0.8
        y0 = 0.5

        r = sqrt((ps.xpg[i, j, k, l, 1] - x0)^2 + (ps.xpg[i, j, k, l, 2] - y0)^2)
        η = r / rc

        δu = κ * η * exp(μ * (1 - η^2)) * (ps.xpg[i, j, k, l, 2] - y0) / r
        δv = -κ * η * exp(μ * (1 - η^2)) * (ps.xpg[i, j, k, l, 1] - x0) / r
        δT = -(ks.gas.γ - 1) * κ^2 / (8 * μ * ks.gas.γ) * exp(2 * μ * (1 - η^2))

        T0 = 1 / prim[end]
        ρ = prim[1]^(ks.gas.γ - 1) * (T0 + δT) / T0^(1 / (ks.gas.γ - 1))
        prim1 = [ρ, prim[2] + δu, prim[3] + δv, 1 / (1 / prim[4] + δT)]

        if r <= rc * 8
            prim .= prim1
        end
    end

    u0[i, j, k, l, :] .= prim_conserve(prim, ks.gas.γ)
end

# alternative formulation from dsmc shock-vortex interaction paper
#=for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
    t1, t2 = ib_rh(ks.gas.Ma, ks.gas.γ)

    if ps.x[i, j] <= ps.x1 * 0.25
        prim = [t2[1], t1[2] - t2[2], 0.0, t2[3]]
    else
        prim = [t1[1], 0.0, 0.0, t1[3]]
        s = prim[1]^(1 - ks.gas.γ) / (2 * prim[end])

        x0 = 1.0
        y0 = 0.5
        r = sqrt((ps.xpg[i, j, k, l, 1] - x0)^2 + (ps.xpg[i, j, k, l, 2] - y0)^2)

        vm = 0.1
        r1 = 0.07
        r2 = 2 * r1

        rset = r1 / (r1^2 - r2^2)
        gammaset = (ks.gas.γ - 1) / ks.gas.γ
        cst = 0.1#5
        rho1 =
            (
                (ks.gas.γ - 1) / cst / ks.gas.γ *
                vm^2 *
                rset^2 *
                (
                    (0.5 * r1^2 - 0.5 * r2^4 / r1^2 - 2 * r2^2 * log(r1)) -
                    (0.5 * r2^2 - 0.5 * r2^4 / r2^2 - 2 * r2^2 * log(r2))
                ) + 1^(ks.gas.γ - 1)
            )^(1 / (ks.gas.γ - 1))
        T1 =
            1.0 +
            2 *
            gammaset *
            rset^2 *
            vm^2 *
            (
                (0.5 * r1^2 - 0.5 * r2^4 / r1^2 - 2 * r2^2 * log(r1)) -
                (0.5 * r2^2 - 0.5 * r2^4 / r2^2 - 2 * r2^2 * log(r2))
            )

        if r <= r1
            vθ = vm * r / r1
        elseif r1 < r <= r2
            vθ = vm * (r1) / (r1^2 - r2^2) * (r - r2^2 / r)
        else
            vθ = 0
        end

        @assert vθ >= 0

        cosθ = (ps.xpg[i, j, k, l, 1] - x0) / r
        sinθ = (ps.xpg[i, j, k, l, 2] - y0) / r

        δu = vθ * sinθ
        δv = -vθ * cosθ
        prim[2:3] .= [δu, δv]

        if r <= r1
            prim[1] =
                (
                    ((ks.gas.γ - 1) / cst / ks.gas.γ * vm^2 / 2 / r1^2 * (r^2 - r1^2)) +
                    rho1^(ks.gas.γ - 1)
                )^(1 / (ks.gas.γ - 1))

            prim[4] = T1 + 2 * gammaset * vm^2 / r1^2 / 2 * (r^2 - r1^2)
            prim[4] = 1 / prim[4]
        elseif r1 < r <= r2
            prim[1] =
                (
                    (ks.gas.γ - 1) / cst / ks.gas.γ *
                    vm^2 *
                    rset^2 *
                    (
                        (0.5 * r^2 - 0.5 * r2^4 / r^2 - 2 * r2^2 * log(r)) -
                        (0.5 * r2^2 - 0.5 * r2^4 / r2^2 - 2 * r2^2 * log(r2))
                    ) + 1^(ks.gas.γ - 1)
                )^(1 / (ks.gas.γ - 1))

            prim[4] =
                1.0 +
                2 *
                gammaset *
                rset^2 *
                vm^2 *
                (
                    (0.5 * r^2 - 0.5 * r2^4 / r^2 - 2 * r2^2 * log(r)) -
                    (0.5 * r2^2 - 0.5 * r2^4 / r2^2 - 2 * r2^2 * log(r2))
                )
            prim[4] = 1 / prim[4]
        end
    end

    u0[i, j, k, l, :] .= prim_conserve(prim, ks.gas.γ)
end=#

prob = ODEProblem(dudt!, u0, tspan, p)
itg = init(prob, Midpoint(), save_everystep = false, adaptive = false, dt = dt)

@showprogress for iter = 1:100#nt
    # limiter
    @inbounds @threads for j = 1:ps.ny
        for i = 1:ps.nx
            ũ = @view itg.u[i, j, :, :, :]
            positive_limiter(ũ, ks.gas.γ, ps.wp ./ 4, ps.ll, ps.lr)
        end
    end

    step!(itg)

    # filter
    for i in axes(itg.u, 1), j in axes(itg.u, 2)
        û = ps.iV * itg.u[i, j, 1:3, 1:3, 1][:]
        su = û[end]^2 / sum(û .^ 2)
        isShock = shock_detector(log10(su), ps.deg, -3.0 * log10(ps.deg), 9)

        if isShock
            for s = 1:4
                û = ps.iV * itg.u[i, j, 1:3, 1:3, s][:]
                FR.modal_filter!(û, 5e-4; filter = :l2)
                uNode = reshape(ps.V * û, 3, 3)
                itg.u[i, j, :, :, s] .= uNode
            end
        end
    end

    # boundary
    itg.u[:, 0, :, :, :] .= itg.u[:, 1, :, :, :]
    itg.u[:, ps.ny+1, :, :, :] .= itg.u[:, ps.ny, :, :, :]
    itg.u[ps.nx+1, :, :, :, :] .= itg.u[ps.nx, :, :, :, :]
    #itg.u[:, 0, :, 1:end, :] .= itg.u[:, 1, :, end:-1:1, :]
    #itg.u[:, ps.ny+1, :, 1:end, :] .= itg.u[:, ps.ny, :, end:-1:1, :]
    #itg.u[:, 0, :, 1:end, 3] .*= -1
    #itg.u[:, ps.ny+1, :, end:-1:1, 3] .*= -1
end

begin
    x = zeros(ps.nx * (ps.deg + 1), ps.ny * (ps.deg + 1))
    y = zeros(ps.nx * (ps.deg + 1), ps.ny * (ps.deg + 1))
    sol = zeros(ps.nx * (ps.deg + 1), ps.ny * (ps.deg + 1), 4)

    for i = 1:ps.nx, j = 1:ps.ny
        idx0 = (i - 1) * (ps.deg + 1)
        idy0 = (j - 1) * (ps.deg + 1)

        for k = 1:ps.deg+1, l = 1:ps.deg+1
            idx = idx0 + k
            idy = idy0 + l
            x[idx, idy] = ps.xpg[i, j, k, l, 1]
            y[idx, idy] = ps.xpg[i, j, k, l, 2]

            sol[idx, idy, :] .= conserve_prim(itg.u[i, j, k, l, :], ks.gas.γ)
            sol[idx, idy, 4] = 1 / sol[idx, idy, 4]
        end
    end
end

contourf(x[:, 1], y[1, :], sol[:, :, 1]', aspect_ratio = 1, legend = true)
plot(x[:, 1], sol[:, end÷2+1, 1])


text:
# This file was generated, do not modify it. # hide
roc_curve = roc(ŷ, y[validation])
plt = scatter(roc_curve, legend=false)
plot!(plt, xlab="false positive rate", ylab="true positive rate")
plot!([0, 1], [0, 1], linewidth=2, linestyle=:dash, color=:black)

savefig(joinpath(@OUTPUT, "EX-telco-roc.svg")) # hide

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using .MaxFlow

const data = [
    (in = """
        4 5 0 3
        0 1 10
        1 2 1
        1 3 1
        0 2 1
        2 3 10
        """,
    out = """
        4 3 5
        0 1 2
        0 2 1
        1 2 1
        1 3 1
        2 3 2
        """),
    (in = """
        2 1 0 1
        0 1 100000
        """,
    out = """
        2 100000 1
        0 1 100000
        """),
    (in = """
        2 1 1 0
        0 1 100000
        """,
    out = """
        2 0 0
        """)
]

@testset "MaxFlow" begin
@testset "Sample $i" for i in keys(data)
    ioin = IOBuffer(data[i].in)
    ioout = IOBuffer()
    MaxFlow.main(ioin, ioout)
    @test String(take!(ioout)) == data[i].out
end
end


text:
using Color, FixedPointNumbers, Images

@everywhere function slicefill(I, f, roots, colors, maxiter, eps, bounds, steps)
	remin, remax, immin, immax = bounds
	redelta, imdelta = steps

	yrange = I[1]
	xrange = I[2]

	slice = (size(yrange,1), size(xrange,1))
	tmp = Array(Vector{Float64}, slice)

	xmin = xrange[1]
	ymin = yrange[1]

	for x=xrange, y=yrange
		tmp[y-ymin+1, x-xmin+1] = f((remin + (x - 1/2)*redelta) + (immin + (y - 1/2)*imdelta)*1im, roots, colors, maxiter, eps)
	end

	return tmp
end

function paranewton(roots, colors, f, maxiter, eps, bounds, width, height)
	remin, remax, immin, immax = bounds
	steps = ((remax-remin)/width, (immax-immin)/height)

	mat = DArray(I -> slicefill(I, f, roots, colors, maxiter, eps, bounds, steps), (height, width))

	return convert(Array, mat)
end

function getimage(matrix)
	rgbmatrix = map(rgb -> RGB{Float64}(rgb[1],rgb[2],rgb[3]), flipud(matrix))
	return convert(Image, rgbmatrix)
end

function saveimage(matrix, file)
	imwrite(getimage(matrix), file)
end

println("Number of available workers: ", length(workers()))

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


text:
for T in (Float32, Float64)
    a1 = rand(Complex{T}, 10)
    b1 = fft(a1)
    af1 = AFArray(a1)
    bf1 = fft(af1)
    @test bf1 == fft1(af1)
    cf1 = ifft(bf1)
    @test cf1 == ifft1(bf1)
    @test isa(bf1, AFArray)
    @test isa(cf1, AFArray)
    @test eltype(bf1) == eltype(cf1) == Complex{T}
    @test Array(bf1) ≈ b1
    @test Array(cf1) ≈ a1

    fft!(af1)
    @test eltype(af1) == Complex{T}
    @test Array(af1) ≈ b1

    ifft!(af1)
    @test eltype(af1) == Complex{T}
    @test Array(af1) ≈ a1

    a2 = rand(Complex{T}, 10, 11)
    b2 = fft(a2)
    af2 = AFArray(a2)
    bf2 = fft(af2)
    @test bf2 == fft2(af2)
    cf2 = ifft(bf2)
    @test cf2 == ifft2(bf2)
    @test isa(bf2, AFArray)
    @test isa(cf2, AFArray)
    @test eltype(bf2) == eltype(cf2) == Complex{T}
    @test Array(bf2) ≈ b2
    @test Array(cf2) ≈ a2

    fft2!(af2)
    @test eltype(af2) == Complex{T}
    @test Array(af2) ≈ b2

    ifft2!(af2)
    @test eltype(af2) == Complex{T}
    @test Array(af2) ≈ a2

    a3 = rand(Complex{T}, 10, 11, 12)
    b3 = fft(a3)
    af3 = AFArray(a3)
    bf3 = fft(af3)
    @test bf3 == fft3(af3)
    cf3 = ifft(bf3)
    @test cf3 == ifft3(bf3)
    @test isa(bf3, AFArray)
    @test isa(cf3, AFArray)
    @test eltype(bf3) == eltype(cf3) == Complex{T}
    @test Array(bf3) ≈ b3
    @test Array(cf3) ≈ a3

    fft3!(af3)
    @test eltype(af3) == Complex{T}
    @test Array(af3) ≈ b3

    ifft3!(af3, 1 / length(af3))
    @test eltype(af3) == Complex{T}
    @test Array(af3) ≈ a3

    ar1 = rand(T, 10)
    br1 = rfft(ar1)
    arf1 = AFArray(ar1)
    brf1 = rfft(arf1)
    crf1 = irfft(brf1, 1 / length(arf1), isodd(size(arf1, 1)))
    @test isa(brf1, AFArray)
    @test isa(crf1, AFArray)
    @test eltype(brf1) == Complex{T}
    @test eltype(crf1) == T
    @test Array(brf1) ≈ br1
    @test Array(crf1) ≈ ar1

    ar2 = rand(T, 10, 11)
    arf2 = AFArray(ar2)
    brf2 = rfft2(arf2)
    crf2 = irfft2(brf2, 1 / length(arf2), isodd(size(arf2, 1)))
    @test isa(brf2, AFArray)
    @test isa(crf2, AFArray)
    @test isa(bf1, AFArray)
    @test isa(cf1, AFArray)
    @test eltype(brf2) == Complex{T}
    @test eltype(crf2) == T
    @test Array(brf2) ≈ rfft(ar2)
    @test Array(crf2) ≈ ar2

    ar3 = rand(T, 11, 12, 10)
    arf3 = AFArray(ar3)
    brf3 = rfft3(arf3)
    crf3 = irfft3(brf3, 1 / length(arf3), isodd(size(arf3, 1)))
    @test isa(brf3, AFArray)
    @test isa(crf3, AFArray)
    @test eltype(brf3) == Complex{T}
    @test eltype(crf3) == T
    @test Array(brf3) ≈ rfft(ar3)
    @test Array(crf3) ≈ ar3

    for sz in ((10,), (10, 10), (10, 10, 10))
        a = rand(Complex{T}, sz...)
        ad = AFArray(a)
        b = fft(a)
        bd = fft(ad)
        @test b ≈ Array(bd)
        @test fft!(a) ≈ Array(fft!(ad))
        @test ifft!(a) ≈ Array(ifft!(ad))

        a = rand(T, sz...)
        ad = AFArray(a)
        b = fft(a)
        bd = fft(ad)
        @test b ≈ Array(bd)
        @test rfft(a) ≈ Array(rfft(ad))
    end
end