text:
"""
```
eqcond(m::COTHANK)
```
 Expresses the equilibrium conditions in canonical form using Γ0, Γ1, C, Ψ, and Π matrices.
Using the mappings of states/equations to integers defined in jpt.jl, coefficients are
specified in their proper positions.

Γ0 (n_states x n_states) holds coefficients of current time states.
Γ1 (n_states x n_states) holds coefficients of lagged states.
C  (n_states x 1) is a vector of constants
Ψ  (n_states x n_shocks_exogenous) holds coefficients of iid shocks.
Π  (n_states x n_states_expectational) holds coefficients of expectational states.
"""
function eqcond(m::COTHANK)
    endo = m.endogenous_states
    exo  = m.exogenous_shocks
    ex   = m.expected_shocks
    eq   = m.equilibrium_conditions

    Γ0 = zeros(n_states(m), n_states(m))
    Γ1 = zeros(n_states(m), n_states(m))
    C  = zeros(n_states(m))
    Ψ  = zeros(n_states(m), n_shocks_exogenous(m))
    Π  = zeros(n_states(m), n_shocks_expectational(m))

    ### ENDOGENOUS STATES ###

    ### 1. Intermediate goods producers

    # Sticky prices
    Γ0[eq[:eq_KL_ratio_good1], endo[:klR1_t]] = -1.
    Γ0[eq[:eq_KL_ratio_good1], endo[:w1_t]]   = 1.
    Γ0[eq[:eq_KL_ratio_good1], endo[:ρ_t]]    = -1.

    Γ0[eq[:eq_mc_good1], endo[:mc1_t]]    = -1.
    Γ0[eq[:eq_mc_good1], endo[:ρ_t]]      = m[:α]
    Γ0[eq[:eq_mc_good1], endo[:w1_t]]     = (1. - m[:α])
    Γ0[eq[:eq_mc_good1], endo[:price1_t]] = -1.

    Γ0[eq[:eq_KL_ratio_good2], endo[:klR2_t]] = -1.
    Γ0[eq[:eq_KL_ratio_good2], endo[:w2_t]]   = 1.
    Γ0[eq[:eq_KL_ratio_good2], endo[:ρ_t]]    = -1.

    Γ0[eq[:eq_mc_good2], endo[:mc2_t]]   = -1.
    Γ0[eq[:eq_mc_good2], endo[:ρ_t]]     = m[:α]
    Γ0[eq[:eq_mc_good2], endo[:w2_t]]    = (1. - m[:α])
    Γ0[eq[:eq_mc_good2], endo[:a_t]]     = -(1. - m[:α])
    Γ0[eq[:eq_mc_good2], endo[:price2_t]]     = -1.

#=
    # Flexible prices
    Γ0[eq[:eq_klR1_f], endo[:klR1_t_f]] = -1.
    Γ0[eq[:eq_klR1_f], endo[:w1_t_f]]   = 1.
    Γ0[eq[:eq_klR1_f], endo[:ρ_t_f]]    = -1.

    Γ0[eq[:eq_mc1_f], endo[:mc1_f_t]]     = -1
    Γ0[eq[:eq_mc1_f], endo[:ρ_t_f]]      = m[:α]
    Γ0[eq[:eq_mc1_f], endo[:w1_t_f]]     = (1 - m[:α])

    Γ0[eq[:eq_klR2_f], endo[:klR2_f_t]] = -1
    Γ0[eq[:eq_klR2_f], endo[:w2_f_t]]   = 1
    Γ0[eq[:eq_klR2_f], endo[:ρ_f_t]]    = -1

    Γ0[eq[:eq_mc2_f], endo[:mc2_f_t]]   = -1
    Γ0[eq[:eq_mc2_f], endo[:ρ_f_t]]   = m[:α]
    Γ0[eq[:eq_mc2_f], endo[:w2_f_t]]   = (1 - m[:α])
    Γ0[eq[:eq_mc2_f], endo[:a_f_t]]   = (1 - m[:α])
=#
    ### 2. Production of goods 1 and 2

    # Sticky prices
    Γ0[eq[:eq_production_good1], endo[:y1_t]]   = -1.
    Γ0[eq[:eq_production_good1], endo[:klR1_t]] = (1. + m[:λ_p_ss]) * m[:α]
    Γ0[eq[:eq_production_good1], endo[:L1_t]]   = (1. + m[:λ_p_ss])

    Γ0[eq[:eq_production_good2], endo[:y2_t]]   = -1
    Γ0[eq[:eq_production_good2], endo[:klR2_t]] = (1. + m[:λ_p_ss]) * m[:α]
    Γ0[eq[:eq_production_good2], endo[:L2_t]]   = (1. + m[:λ_p_ss])
    Γ0[eq[:eq_production_good2], endo[:a_t]]    = -m[:λ_p_ss]

#=
    # Flexible prices
    Γ0[eq[:eq_y1_f], endo[:y1_f_t]] = -1
    Γ0[eq[:eq_y1_f], endo[:klR1_f_t]] = (1 + m[:λ_p]) * m[:α]
    Γ0[eq[:eq_y1_f], endo[:L1_f_t]] = (1 + m[:λ_p])
    Γ0[eq[:eq_y1_f], endo[:a_f_t]] = -m[:λ_p]
=#
    ### 3. Production of final gods

    # Sticky prices
    Γ0[eq[:eq_production_finalgood], endo[:y_t]] = -1.
    Γ0[eq[:eq_production_finalgood], endo[:y1_t]] = m[:v]
    Γ0[eq[:eq_production_finalgood], endo[:y2_t]] = (1 - m[:v])
    Γ0[eq[:eq_production_finalgood], endo[:d_t]] = -(1 - m[:v]) / (1 - m[:ζ])

    Γ0[eq[:eq_price_index], endo[:π_t]]  = -1.
    Γ0[eq[:eq_price_index], endo[:d_t]]  = (1. - m[:v]) / (1. - m[:ζ])
    Γ0[eq[:eq_price_index], endo[:π1_t]] = m[:v]
    Γ0[eq[:eq_price_index], endo[:π2_t]] = (1. - m[:v])
    Γ1[eq[:eq_price_index], endo[:d_t]]  = (1. - m[:v]) / (1. - m[:ζ])

    Γ0[eq[:eq_demand_good1], endo[:y1_t]] = 1
    Γ0[eq[:eq_demand_good1], endo[:y2_t]] = -1
    Γ0[eq[:eq_demand_good1], endo[:π1_t]] = m[:ζ]
    Γ0[eq[:eq_demand_good1], endo[:π2_t]] = -m[:ζ]
    Γ0[eq[:eq_demand_good1], endo[:d_t]]  = 1
    Γ1[eq[:eq_demand_good1], endo[:y1_t]] = 1
    Γ1[eq[:eq_demand_good1], endo[:y2_t]] = -1
    Γ1[eq[:eq_demand_good1], endo[:d_t]]  = 1

    Γ0[eq[:eq_π1], endo[:y1_t]]     = 1.
    Γ0[eq[:eq_π1], endo[:price1_t]] = m[:vars]
    Γ0[eq[:eq_π1], endo[:y_t]] = -1.

    Γ0[eq[:eq_π2], endo[:y2_t]]     = 1.
    Γ0[eq[:eq_π2], endo[:d_t]]      = -1.
    Γ0[eq[:eq_π2], endo[:price2_t]] = m[:vars]
    Γ0[eq[:eq_π2], endo[:y_t]]      = -1.

    ### 4. Price Phillips Curve

    # Sticky prices
    κ = (1. - m[:ξ_p] * m[:β]) * (1. - m[:ξ_p]) / ( m[:ξ_p] * (1. + m[:ι_p] * m[:β]))
    Γ0[eq[:eq_price_phillips_good1], endo[:π1_t]] = -1.
    Γ0[eq[:eq_price_phillips_good1], endo[:Eπ1_t]] = m[:β] / (1. + m[:ι_p] * m[:β])
    Γ0[eq[:eq_price_phillips_good1], endo[:mc1_t]] = κ
    Γ0[eq[:eq_price_phillips_good1], endo[:λ_p_t]] = 1. # IN THE NOTES, IT SAYS WE REDEFINE λ̂_p_t = κ λ̂_p_t, DOES THIS APPLY THE REST OF THE DOCUMENT OR JUST FOR THE PRICE_PHILLIPS CURVE EQUATIONS?
    Γ1[eq[:eq_price_phillips_good1], endo[:π1_t]]  = -m[:ι_p] / (1. + m[:ι_p] * m[:β])

    Γ0[eq[:eq_price_phillips_good2], endo[:π2_t]]  = -1.
    Γ0[eq[:eq_price_phillips_good2], endo[:Eπ2_t]] = m[:β] / (1. + m[:ι_p] * m[:β])
    Γ0[eq[:eq_price_phillips_good2], endo[:mc2_t]] = κ
    Γ0[eq[:eq_price_phillips_good2], endo[:λ_p_t]] = 1.
    Γ1[eq[:eq_price_phillips_good2], endo[:π2_t]]  = -m[:ι_p] / (1. + m[:ι_p] * m[:β])

    # Flexible prices
#    Γ0[eq[:eq_R_f], endo[:mc_f_t]] = 1 # THIS CAME FROM THANK, CHECK IF IT'S RIGHT

    ### 5. Consumption
    expγ = exp(m[:γ])

    # Sticky prices
    Γ0[eq[:eq_marg_utility_S], endo[:λ_S_t]] = -1.
    Γ0[eq[:eq_marg_utility_S], endo[:b_t]]   = 1.
    Γ0[eq[:eq_marg_utility_S], endo[:z_t]]   = -m[:h] / (expγ - m[:h])
    Γ0[eq[:eq_marg_utility_S], endo[:c_S_t]] = -expγ / (expγ - m[:h])
    Γ1[eq[:eq_marg_utility_S], endo[:c_S_t]] = -m[:h] / (expγ - m[:h])
#=
    # Flexible
    Γ0[eq[:eq_marg_utility_S_f], endo[:λ_S_f_t]] = -1.
    Γ0[eq[:eq_marg_utility_S_f], endo[:b_f_t]]   = 1.
    Γ0[eq[:eq_marg_utility_S_f], endo[:z_f_t]]   = - m[:h] / (expγ - m[:h])
    Γ0[eq[:eq_marg_utility_S_f], endo[:c_S_f_t]] = -expγ / (expγ - m[:h])
    Γ1[eq[:eq_marg_utility_S_f], endo[:C_S_f_t]] = m[:h] / (expγ - m[:h])
=#
    # Sticky
    Γ0[eq[:eq_marg_utility_H1], endo[:λ_H1_t]] = -1.
    Γ0[eq[:eq_marg_utility_H1], endo[:b_t]]    = 1.
    Γ0[eq[:eq_marg_utility_H1], endo[:z_t]]    = -m[:h] / (expγ - m[:h])
    Γ0[eq[:eq_marg_utility_H1], endo[:c_H1_t]] = -expγ / (expγ - m[:h])
    Γ1[eq[:eq_marg_utility_H1], endo[:c_H1_t]] = -m[:h] / (expγ - m[:h])
#=
    # Flexible
    Γ0[eq[:eq_marg_utility_H1_f], endo[:λ_H1_f_t]] = -1.
    Γ0[eq[:eq_marg_utility_H1_f], endo[:b_f_t]]   = 1.
    Γ0[eq[:eq_marg_utility_H1_f], endo[:z_f_t]]   = - m[:h] / (expγ - m[:h])
    Γ0[eq[:eq_marg_utility_H1_f], endo[:c_H1_f_t]] = -expγ / (expγ - m[:h])
    Γ1[eq[:eq_marg_utility_H1_f], endo[:c_H1_f_t]] = m[:h] / (expγ - m[:h])
=#
    # Sticky
    Γ0[eq[:eq_marg_utility_H2], endo[:λ_H2_t]] = -1.
    Γ0[eq[:eq_marg_utility_H2], endo[:b_t]]    = 1.
    Γ0[eq[:eq_marg_utility_H2], endo[:z_t]]    = -m[:h] / (expγ - m[:h])
    Γ0[eq[:eq_marg_utility_H2], endo[:c_H2_t]] = -expγ / (expγ - m[:h])
    Γ1[eq[:eq_marg_utility_H2], endo[:c_H2_t]] = -m[:h] / (expγ - m[:h])
#=
    # Flexible
    Γ0[eq[:eq_marg_utility_H2_f], endo[:λ_H2_f_t]] = -1
    Γ0[eq[:eq_marg_utility_H2_f], endo[:b_f_t]]   = 1
    Γ0[eq[:eq_marg_utility_H2_f], endo[:z_f_t]]   = - m[:h] / (expγ - m[:h])
    Γ0[eq[:eq_marg_utility_H2_f], endo[:c_H2_f_t]] = -expγ / (expγ - m[:h])
    Γ1[eq[:eq_marg_utility_H2_f], endo[:c_H2_f_t]] = m[:h] / (expγ - m[:h])
=#
    # Sticky

    Γ0[eq[:eq_average_c], endo[:c_t]]   = -1
    Γ0[eq[:eq_average_c], endo[:c_S_t]] = (1 - m[:θ]) * m[:c_S_ss] / m[:c_ss]
    Γ0[eq[:eq_average_c], endo[:c_H1_t]] = m[:θ] * m[:f_H1] * m[:c_H1_ss] / m[:c_ss]
    Γ0[eq[:eq_average_c], endo[:c_H2_t]] = m[:θ] * (1 - m[:f_H1]) * m[:c_H2_ss] / m[:c_ss]
#=
    # Flexible
    Γ0[eq[:eq_average_c_f], endo[:c_f_t]]   = -1
    Γ0[eq[:eq_average_c_f], endo[:c_S_t]] = (1 - m[:θ]) * m[:c_S_ss] / m[:c_ss]
    Γ0[eq[:eq_average_c_f], endo[:c_H1_t]] = m[:θ] * m[:f_H1] * m[:c_H1_ss] / m[:c_ss]
    Γ0[eq[:eq_average_c_f], endo[:c_H2_t]] = m[:θ] * (1 - m[:f_H1]) * m[:c_H2_ss] / m[:c_ss]
=#
    # Sticky
    Γ0[eq[:eq_euler], endo[:λ_S_t]]   = -1. # Euler equation also pins down the bonds demand/self-insurance motive
    Γ0[eq[:eq_euler], endo[:R_t]]     = 1.
    Γ0[eq[:eq_euler], endo[:z_t]]    = -m[:ρ_z]
    Γ0[eq[:eq_euler], endo[:Eπ_t]]    = -1.
    Γ0[eq[:eq_euler], endo[:Eλ_S_t]]  = m[:β] * m[:R_ss] / (expγ * m[:π_ss]) * m[:s]
    Γ0[eq[:eq_euler], endo[:Eλ_H1_t]] =  m[:β] * m[:R_ss] / (expγ * m[:π_ss]) *
                                         m[:f_S1] * m[:λ_H1_ss] * (1. - m[:s]) / m[:λ_S_ss]
    Γ0[eq[:eq_euler], endo[:Eλ_H2_t]] =  m[:β] * m[:R_ss] / (expγ * m[:π_ss]) * (1. - m[:f_S1]) *
                                         m[:λ_H2_ss] * (1. - m[:s]) / m[:λ_S_ss]
    Γ0[eq[:eq_euler], endo[:Ex_t]]    = m[:β] * m[:R_ss] / (expγ * m[:π_ss]) * m[:s] *
                                        (1. - m[:f_S1] * m[:λ_H1_ss] / m[:λ_S_ss] -
                                        (1. - m[:f_S1]) * m[:λ_H2_ss] / m[:λ_S_ss]) *
                                        m[:s_x]
#=
    # Flexible
    Γ0[eq[:eq_euler_f], endo[:λ_S_f_t]] = -1.
    Γ0[eq[:eq_euler_f], endo[:R_f_t]] = 1.
    Γ0[eq[:eq_euler_f], endo[:Ez_f_t]] = -1.
    Γ0[eq[:eq_euler_f], endo[:Eπ_f_t]] = -1.
    Γ0[eq[:eq_euler_f], endo[:Eλ_S_f_t]] = m[:β] * m[:R_ss] / (expγ * m[:π_ss]) * m[:s]
    Γ0[eq[:eq_euler_f], endo[:Eλ_H1_f_t]] =  m[:β] * m[:R_ss] / (expγ * m[:π_ss]) *
                                           m[:f_S1] * m[:λ_H1_ss] * (1. - m[:s]) / m[:λ_S_ss]
    Γ0[eq[:eq_euler_f], endo[:Eλ_H2_f_t]] =  m[:β] * m[:R_ss] / (expγ * m[:π_ss]) * (1. - m[:f_S1]) *
                                           m[:λ_H2_ss] * (1. - m[:s]) / m[:λ_S_ss]
    Γ0[eq[:eq_euler_f], endo[:Ex_f_t]] = m[:β] * m[:R_ss] / (expγ * m[:π_ss]) * m[:s] *
                                       (1. - m[:f_S1] * m[:λ_H1_ss] / m[:λ_S_ss] -
                                       (1. - m[:f_S1]) * m[:λ_H2_ss] / m[:λ_S_ss]) *
                                       m[:s_x]
=#
    # Sticky
    Γ0[eq[:eq_c_H1], endo[:c_H1_t]] = -1.
    Γ0[eq[:eq_c_H1], endo[:w1_t]]   = m[:w1_ss] * m[:H1] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1], endo[:L1_t]]   = m[:w1_ss] * m[:H1] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1], endo[:t_H1_t]] = -m[:x_ss] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1], endo[:τ_H1_t]] = m[:x_ss] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1], endo[:z_t]]    = -m[:f_S1] * (1. - m[:s]) / (m[:θ] * m[:f_H1]) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1], endo[:π_t]]    = -m[:f_S1] * (1. - m[:s]) / (m[:θ] * m[:f_H1]) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1], endo[:x_t]]    = -m[:f_S1] * m[:s] / (m[:θ] * m[:f_H1]) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss] *
                                       m[:s_x]
    Γ1[eq[:eq_c_H1], endo[:R_t]]    = -m[:f_S1] * (1. - m[:s]) / (m[:θ] * m[:f_H1]) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss]
    Γ1[eq[:eq_c_H1], endo[:bᴿ_t]]   = -m[:f_S1] * (1. - m[:s]) / (m[:θ] * m[:f_H1]) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss] *
                                       m[:x_ss] / m[:bᴿ_ss]
#=
    # Flexible
    Γ0[eq[:eq_c_H1_f], endo[:c_f_H1_t]] = -1
    Γ0[eq[:eq_c_H1_f], endo[:w1_f_t]] = m[:w1_ss] * m[:H1] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1_f], endo[:H1_f_t]] = m[:w1_ss] * m[:H1] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1_f], endo[:t_H1_f_t]] = -m[:x_ss] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1_f], endo[:τ_H1_f_t]] = m[:x_ss] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1_f], endo[:z_f_t]] = -m[:f_S1] * (1 - m[:s]) / (m[:θ] * m[:f_H1]) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1_f], endo[:π_f_t]] = -m[:f_S1] * (1 - m[:s]) / (m[:θ] * m[:f_H1]) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss]
    Γ0[eq[:eq_c_H1_f], endo[:x_f_t]] = -m[:f_S1] * m[:s] / (m[:θ] * m[:f_H1]) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss] *
                                        m[:s_x]
    Γ1[eq[:eq_c_H1_f], endo[:R_f_t]] = -m[:f_S1] * (1 - m[:s]) / (m[:θ] * m[:f_H1]) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss]
    Γ1[eq[:eq_c_H1_f], endo[:bᴿ_f_t]] = -m[:f_S1] * (1 - m[:s]) / (m[:θ] * m[:f_H1]) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H1_ss] *
                                        m[:x_ss] / m[:bᴿ_ss]
=#
    # Sticky
    Γ0[eq[:eq_c_H2], endo[:c_H2_t]] = -1
    Γ0[eq[:eq_c_H2], endo[:w2_t]]   = m[:w2_ss] * m[:H2] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2], endo[:L2_t]]   = m[:w2_ss] * m[:H2] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2], endo[:t_H2_t]] = -m[:x_ss] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2], endo[:τ_H2_t]] = m[:x_ss] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2], endo[:z_t]]    = -(1. - m[:f_S1]) * (1 - m[:s]) / (m[:θ] * (1. - m[:f_H1])) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2], endo[:π_t]]    = -(1. - m[:f_S1]) * (1 - m[:s]) / (m[:θ] * (1. - m[:f_H1])) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2], endo[:x_t]]    = -(1. - m[:f_S1]) * m[:s] / (m[:θ] * (1. - m[:f_H1])) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss] *
                                       m[:s_x]
    Γ1[eq[:eq_c_H2], endo[:R_t]]    = -(1. - m[:f_S1]) * (1 - m[:s]) / (m[:θ] * (1. - m[:f_H1])) *
                                       m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss]
    Γ1[eq[:eq_c_H2], endo[:bᴿ_t]]    = -(1. - m[:f_S1]) * (1 - m[:s]) / (m[:θ] * (1. - m[:f_H1])) *
                                        m[:R_ss] / (expγ * m[:π_ss] *  m[:c_H2_ss]) *
                                        m[:x_ss]
#=
    # Flexible
    Γ0[eq[:eq_c_H2_f], endo[:c_f_H2_t]] = -1
    Γ0[eq[:eq_c_H2_f], endo[:w2_f_t]] = m[:w2_ss] * m[:H2] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2_f], endo[:H2_f_t]] = m[:w2_ss] * m[:H2] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2_f], endo[:t_H2_f_t]] = -m[:x_ss] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2_f], endo[:τ_H2_f_t]] = m[:x_ss] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2_f], endo[:z_f_t]] = -(1. - m[:f_S1]) * (1 - m[:s]) / (m[:θ] * (1. - m[:f_H1])) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2_f], endo[:π_f_t]] = -(1. - m[:f_S1]) * (1 - m[:s]) / (m[:θ] * (1. - m[:f_H1])) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss]
    Γ0[eq[:eq_c_H2_f], endo[:x_f_t]] = -(1. - m[:f_S1]) * m[:s] / (m[:θ] * (1. - m[:f_H1])) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss] *
                                        m[:s_x]
    Γ1[eq[:eq_c_H2_f], endo[:R_f_t]] = -(1. - m[:f_S1]) * (1 - m[:s]) / (m[:θ] * (1. - m[:f_H1])) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss]
    Γ1[eq[:eq_c_H2_f], endo[:bᴿ_f_t]] = -(1. - m[:f_S1]) * (1 - m[:s]) / (m[:θ] * (1. - m[:f_H1])) *
                                        m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / m[:c_H2_ss] *
                                        m[:x_ss] / m[:bᴿ_ss]
=#

    ### 6. Investment and capital accumulation

    # Sticky prices
    Γ0[eq[:eq_cap_util], endo[:ρ_t]] = 1.
    Γ0[eq[:eq_cap_util], endo[:u_t]] = -m[:χ]

    Γ0[eq[:eq_ϕ], endo[:ϕ_t]]    = -1.
    Γ0[eq[:eq_ϕ], endo[:Eϕ_t]]   = (1. - m[:δ]) * m[:β] * exp(-m[:γ])
    Γ0[eq[:eq_ϕ], endo[:z_t]]   = -m[:ρ_z]
    Γ0[eq[:eq_ϕ], endo[:Eλ_S_t]] = (1. - (1. - m[:δ]) * m[:β] * exp(-m[:γ]))
    Γ0[eq[:eq_ϕ], endo[:Eρ_t]]   = (1. - (1. - m[:δ]) * m[:β] * exp(-m[:γ]))

    Γ0[eq[:eq_i], endo[:λ_S_t]]  = -1.
    Γ0[eq[:eq_i], endo[:ϕ_t]]    = 1.
    Γ0[eq[:eq_i], endo[:μ_t]]    = 1.
    Γ0[eq[:eq_i], endo[:i_S_t]]  = -exp(2. * m[:γ]) * m[:S′′] * (1. + m[:β])
    Γ0[eq[:eq_i], endo[:z_t]]    = -exp(2. * m[:γ]) * m[:S′′] * (1 - m[:β] * m[:ρ_z])
    Γ0[eq[:eq_i], endo[:Ei_S_t]] = m[:β] * exp(2. * m[:γ]) * m[:S′′]
    Γ1[eq[:eq_i], endo[:i_S_t]]  = -exp(2. * m[:γ]) * m[:S′′]

    Γ0[eq[:eq_cap_input], endo[:k_t]]   = -1.
    Γ0[eq[:eq_cap_input], endo[:u_t]]   = 1.
    Γ0[eq[:eq_cap_input], endo[:z_t]]   = -1.
    Γ1[eq[:eq_cap_input], endo[:k_S_t]] = -1.

    Γ0[eq[:eq_cap_accum], endo[:k_S_t]] = 1.
    Γ0[eq[:eq_cap_accum], endo[:z_t]]   = (1. - m[:δ]) * exp(-m[:γ])
    Γ0[eq[:eq_cap_accum], endo[:μ_t]]   = -(1. - (1. - m[:δ]) * exp(-m[:γ]))
    Γ0[eq[:eq_cap_accum], endo[:i_S_t]] = -(1. - (1. - m[:δ]) * exp(-m[:γ]))
    Γ1[eq[:eq_cap_accum], endo[:k_S_t]] = (1. - m[:δ]) * exp(-m[:γ])
#=
    # Flexible prices
    Γ0[eq[:eq_k_S_f], endo[:k_S_f_t]] = 1
    Γ0[eq[:eq_k_S_f], endo[:z_t]] = (1 - m[:δ]) * exp(-m[:γ])
    Γ0[eq[:eq_k_S_f], endo[:μ_t]] = (1 - (1 - m[:δ]) * exp(-m[:γ]))
    Γ0[eq[:eq_k_S_f], endo[:i_S_f_t]] = (1 - (1 - m[:δ]) * exp(-m[:γ]))
    Γ1[eq[:eq_k_S_f], endo[:k_S_f_t]] = (1 - m[:δ]) * exp(-m[:γ])
=#
    ### 7. Wage Phillips Curve

    # Sticky prices
    κ_w = (1. - m[:ξ_w] * m[:β]) * (1. - m[:ξ_w]) / (m[:ξ_w] * (1. + m[:β]) * (1. +
              m[:ν] * (1. + 1. / m[:λ_w_ss])))
    Γ0[eq[:eq_wage_phillips_good1], endo[:w1_t]]   = -1.
    Γ0[eq[:eq_wage_phillips_good1], endo[:Ew1_t]]  = m[:β] / (1. + m[:β])
    Γ0[eq[:eq_wage_phillips_good1], endo[:g_w1_t]] = -κ_w
    Γ0[eq[:eq_wage_phillips_good1], endo[:π_t]]    = -(1. + m[:β] * m[:ι_w]) / (1. + m[:β])
    Γ0[eq[:eq_wage_phillips_good1], endo[:Eπ_t]]   = m[:β] / (1. + m[:β])
    Γ0[eq[:eq_wage_phillips_good1], endo[:z_t]]    = -(1. + m[:β] * m[:ι_w] - m[:ρ_z] * m[:β]) / (1. + m[:β])
    Γ0[eq[:eq_wage_phillips_good1], endo[:λ_w_t]]  = 1.
    Γ1[eq[:eq_wage_phillips_good1], endo[:w1_t]]   = -1. / (1. + m[:β])
    Γ1[eq[:eq_wage_phillips_good1], endo[:π_t]]    = -m[:ι_w] / (1. + m[:β])
    Γ1[eq[:eq_wage_phillips_good1], endo[:z_t]]    = -m[:ι_w] / (1. + m[:β])

    Γ0[eq[:eq_marg_sub_good1], endo[:g_w1_t]]  = 1.
    Γ0[eq[:eq_marg_sub_good1], endo[:w1_t]]    = -1.
    Γ0[eq[:eq_marg_sub_good1], endo[:L1_t]]    = m[:ν]
    Γ0[eq[:eq_marg_sub_good1], endo[:b_t]]     = 1.
    Γ0[eq[:eq_marg_sub_good1], endo[:λ_SH1_t]] = -1.

    Γ0[eq[:eq_avg_marg_util_SH1], endo[:λ_SH1_t]] = -1.
    Γ0[eq[:eq_avg_marg_util_SH1], endo[:λ_S_t]]   = (1. - m[:θ]) * m[:f_S1] / ((1. - m[:θ]) * m[:f_S1] + m[:θ] * m[:f_H1]) *
                                                    m[:λ_S_ss] / m[:λ_SH1_ss]
    Γ0[eq[:eq_avg_marg_util_SH1], endo[:λ_H1_t]]  = m[:θ] * m[:f_H1] / ((1. - m[:θ]) * m[:f_S1] + m[:θ] * m[:f_H1]) *
                                                    m[:λ_H1_ss] / m[:λ_SH1_ss]

    Γ0[eq[:eq_wage_phillips_good2], endo[:w2_t]]   = -1.
    Γ0[eq[:eq_wage_phillips_good2], endo[:Ew2_t]]  = m[:β] / (1. + m[:β])
    Γ0[eq[:eq_wage_phillips_good2], endo[:g_w2_t]] = -κ_w
    Γ0[eq[:eq_wage_phillips_good2], endo[:π_t]]    = -(1. + m[:β] * m[:ι_w]) / (1. + m[:β])
    Γ0[eq[:eq_wage_phillips_good2], endo[:Eπ_t]]   = m[:β] / (1. + m[:β])
    Γ0[eq[:eq_wage_phillips_good2], endo[:z_t]]    = -(1. + m[:β] * m[:ι_w] - m[:ρ_z] * m[:β]) / (1. + m[:β])
    Γ0[eq[:eq_wage_phillips_good2], endo[:λ_w_t]]  = 1.
    Γ1[eq[:eq_wage_phillips_good2], endo[:w2_t]]   = -1. / (1. + m[:β])
    Γ1[eq[:eq_wage_phillips_good2], endo[:π_t]]    = -m[:ι_w] / (1. + m[:β])
    Γ1[eq[:eq_wage_phillips_good2], endo[:z_t]]    = -m[:ι_w] / (1. + m[:β])

    Γ0[eq[:eq_marg_sub_good2], endo[:g_w2_t]]  = 1.
    Γ0[eq[:eq_marg_sub_good2], endo[:w2_t]]    = -1.
    Γ0[eq[:eq_marg_sub_good2], endo[:L2_t]]    = m[:ν]
    Γ0[eq[:eq_marg_sub_good2], endo[:b_t]]     = 1.
    Γ0[eq[:eq_marg_sub_good2], endo[:λ_SH2_t]] = -1.

    Γ0[eq[:eq_avg_marg_util_SH2], endo[:λ_SH2_t]] = -1.
    Γ0[eq[:eq_avg_marg_util_SH2], endo[:λ_S_t]]   = (1. - m[:θ]) * (1. - m[:f_S1]) / ((1. - m[:θ]) *
                                                    (1. - m[:f_S1]) + m[:θ] * (1 - m[:f_H1])) *
                                                    m[:λ_S_ss] / m[:λ_SH2_ss]
    Γ0[eq[:eq_avg_marg_util_SH2], endo[:λ_H2_t]]  = m[:θ] * (1. - m[:f_H1]) / ((1. - m[:θ]) *
                                                    (1. - m[:f_S1]) + m[:θ] * (1. - m[:f_H1])) *
                                                    m[:λ_H2_ss] / m[:λ_SH2_ss]

    ### 8. Monetary and fiscal policy
    Γ0[eq[:eq_mp], endo[:R_t]]    = -1.
    Γ1[eq[:eq_mp], endo[:R_t]]    = -m[:ρ_R]
    Γ0[eq[:eq_mp], endo[:π_t]]    = (1. - m[:ρ_R]) * m[:ϕ_1]
    Γ0[eq[:eq_mp], endo[:x_t]]    = (1. - m[:ρ_R]) * m[:ϕ_3]
    Γ0[eq[:eq_mp], endo[:η_mp_t]] = 1.
    Γ1[eq[:eq_mp], endo[:x_t]]    = (1. - m[:ρ_R]) * m[:ϕ_3]

    Γ0[eq[:eq_revenue], endo[:t_t]]    = -1.
    Γ0[eq[:eq_revenue], endo[:t_S_t]]  = (1. - m[:θ])
    Γ0[eq[:eq_revenue], endo[:t_H1_t]] = m[:θ] * m[:f_H1]
    Γ0[eq[:eq_revenue], endo[:t_H2_t]] = m[:θ] * (1. - m[:f_H1])

    Γ0[eq[:eq_transfer], endo[:τ_t]]    = -1.
    Γ0[eq[:eq_transfer], endo[:τ_S_t]]  = (1. - m[:θ])
    Γ0[eq[:eq_transfer], endo[:τ_H1_t]] = m[:θ] * m[:f_H1]
    Γ0[eq[:eq_transfer], endo[:τ_H2_t]] = m[:θ] * (1. - m[:f_H1])

#=
    Γ0[eq[:eq_realdebt], endo[:bᴿ_t]] = m[:x_ss] / m[:bᴿ_ss]
    Γ0[eq[:eq_realdebt], endo[:π_t]]  = m[:R_ss] / (expγ * m[:π_ss])
    Γ0[eq[:eq_realdebt], endo[:z_t]]  = m[:R_ss] / (expγ * m[:π_ss])
    Γ0[eq[:eq_realdebt], endo[:g_t]]  = -(1. - m[:R_ss] / (expγ * m[:π_ss])) * m[:x_ss] / (m[:g_ss] - m[:t_ss] + m[:τ_ss])
    Γ0[eq[:eq_realdebt], endo[:t_t]]  = (1. - m[:R_ss] / (expγ * m[:π_ss])) * m[:x_ss] / (m[:g_ss] - m[:t_ss] + m[:τ_ss])
    Γ0[eq[:eq_realdebt], endo[:τ_t]]  = -(1. - m[:R_ss] / (expγ * m[:π_ss])) * m[:x_ss] / (m[:g_ss] - m[:t_ss] + m[:τ_ss])
    Γ1[eq[:eq_realdebt], endo[:R_t]]  = m[:R_ss] / (expγ * m[:π_ss])
    Γ1[eq[:eq_realdebt], endo[:bᴿ_t]] = m[:R_ss] / (expγ * m[:π_ss]) * m[:x_ss] / m[:bᴿ_ss]
=#
    Γ0[eq[:eq_realdebt], endo[:bᴿ_t]] = 1.
    Γ0[eq[:eq_realdebt], endo[:π_t]]  = m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / (m[:x_ss])
    Γ0[eq[:eq_realdebt], endo[:z_t]]  = m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / (m[:x_ss])
    Γ0[eq[:eq_realdebt], endo[:g_t]]  = -1.
    Γ0[eq[:eq_realdebt], endo[:t_t]]  = 1.
    Γ0[eq[:eq_realdebt], endo[:τ_t]]  = -1.
    Γ1[eq[:eq_realdebt], endo[:R_t]]  = m[:R_ss] / (expγ * m[:π_ss]) * m[:bᴿ_ss] / (m[:x_ss])
    Γ1[eq[:eq_realdebt], endo[:bᴿ_t]] = m[:R_ss] / (expγ * m[:π_ss])

    Γ0[eq[:eq_fiscal_rule], endo[:t_t]] = -m[:x_ss] / m[:t_ss]
    Γ0[eq[:eq_fiscal_rule], endo[:x_t]] = (1. + (1. - m[:ρ_T]) * m[:ζ_X] - (1. - m[:ρ_T]) * m[:ζ_G])
    Γ0[eq[:eq_fiscal_rule], endo[:g_t]] = (1. - m[:ρ_T]) * m[:ζ_G] * m[:x_ss] / m[:g_ss]
    Γ1[eq[:eq_fiscal_rule], endo[:t_t]] = -m[:ρ_T] * m[:x_ss] / m[:t_ss]
    Γ1[eq[:eq_fiscal_rule], endo[:bᴿ_t]] = -(1. - m[:ρ_T]) * m[:ζ_B]
    Γ1[eq[:eq_fiscal_rule], endo[:x_t]] = m[:ρ_T] + (1. - m[:ρ_T]) * m[:ζ_B]  * m[:bᴿ_ss] / m[:x_ss]

    Γ0[eq[:eq_redistribute_H1], endo[:t_H1_t]] = m[:θ] * m[:f_H1] # IN THE TEXT, THIS IS f_H, but I think it should be f_H1
    Γ0[eq[:eq_redistribute_H1], endo[:t_t]]    = -m[:ψ_H1]

    Γ0[eq[:eq_redistribute_H2], endo[:t_H2_t]] = m[:θ] * (1. - m[:f_H1]) # IN THE TEXT, THIS IS f_H, but I think it should be f_H1
    Γ0[eq[:eq_redistribute_H2], endo[:t_t]]    = -m[:ψ_H2]

    ### 9. Market clearing and aggregate resource constraint
    # commenting out following eq condition b/c of presence of π_S_ss, unsure of parameter's meaning
    Γ0[eq[:eq_real_profits], endo[:π_S_t]]     = -1.
    Γ0[eq[:eq_real_profits], endo[:y_t]]       = m[:y_ss] / ((1. - m[:θ]) * m[:x_ss])
    Γ0[eq[:eq_real_profits], endo[:price1_t]]  = -m[:price1_ss] * m[:mc1_ss] * m[:klR1_ss] * m[:L1] / ((1. - m[:θ]) * m[:x_ss])
    Γ0[eq[:eq_real_profits], endo[:mc1_t]]     = -m[:price1_ss] * m[:mc1_ss] * m[:klR1_ss] * m[:L1] / ((1. - m[:θ]) * m[:x_ss])
    Γ0[eq[:eq_real_profits], endo[:klR1_t]]    = -m[:price1_ss] * m[:mc1_ss] * m[:klR1_ss] * m[:L1] / ((1. - m[:θ]) * m[:π_S_ss])
    Γ0[eq[:eq_real_profits], endo[:L1_t]]      = -m[:price1_ss] * m[:mc1_ss] * m[:klR1_ss] * m[:L1] / ((1. - m[:θ]) * m[:π_S_ss])
    Γ0[eq[:eq_real_profits], endo[:price2_t]]  = -m[:price2_ss] * m[:A2] ^ (1. - m[:α]) * m[:mc2_ss] * m[:klR2_ss] * m[:L2] /
                                               ((1. - m[:θ]) * m[:x_ss])
    Γ0[eq[:eq_real_profits], endo[:mc2_t]]     = -m[:price2_ss] * m[:A2] ^ (1. - m[:α]) * m[:mc2_ss] * m[:klR2_ss] * m[:L2] /
                                               ((1. - m[:θ]) * m[:x_ss])
    Γ0[eq[:eq_real_profits], endo[:klR2_t]]    = -m[:price2_ss] * m[:A2] ^ (1. - m[:α]) * m[:mc2_ss] * m[:klR2_ss] * m[:L2] /
                                               ((1. - m[:θ]) * m[:x_ss])
    Γ0[eq[:eq_real_profits], endo[:L2_t]]      = -m[:price2_ss] * m[:A2] ^ (1. - m[:α]) * m[:mc2_ss] * m[:klR2_ss] * m[:L2] /
                                               ((1. - m[:θ]) * m[:x_ss])

    Γ0[eq[:eq_cap_market_clear], endo[:klR1_t]] = m[:klR1_ss] * m[:L1] / m[:k_ss]
    Γ0[eq[:eq_cap_market_clear], endo[:L1_t]]   = m[:klR1_ss] * m[:L1] / m[:k_ss]
    Γ0[eq[:eq_cap_market_clear], endo[:klR2_t]] = m[:klR2_ss] * m[:L2] / m[:k_ss]
    Γ0[eq[:eq_cap_market_clear], endo[:L2_t]]   = m[:klR2_ss] * m[:L2] / m[:k_ss]
    Γ0[eq[:eq_cap_market_clear], endo[:k_t]]    = -1.

    Γ0[eq[:eq_goods_market_clear], endo[:c_t]]   = m[:c_ss] / m[:y_ss]
    Γ0[eq[:eq_goods_market_clear], endo[:i_S_t]] = (1 - m[:θ]) * m[:i_S_ss] / m[:y_ss]
    Γ0[eq[:eq_goods_market_clear], endo[:g_t]]   = m[:x_ss] / m[:y_ss]
    Γ0[eq[:eq_goods_market_clear], endo[:u_t]]   = m[:ρ_ss] * m[:k_ss] / m[:y_ss]
    Γ0[eq[:eq_goods_market_clear], endo[:y_t]]   = -1

    Γ0[eq[:eq_gdp], endo[:x_t]] = -1.
    Γ0[eq[:eq_gdp], endo[:y_t]] = 1.
    Γ0[eq[:eq_gdp], endo[:u_t]] = -m[:ρ_ss] * m[:k_ss] / m[:y_ss]

    ### 10. Exogenous processes
# May need to re-sign these values b/c want positive MP shock (for example) to be raise in rates.

    # z
    Γ0[eq[:eq_z], endo[:z_t]] = 1
    Γ1[eq[:eq_z], endo[:z_t]] = m[:ρ_z]
    Ψ[eq[:eq_z], exo[:z_sh]]  = 1

    # a
    Γ0[eq[:eq_a], endo[:a_t]] = 1.
    Γ1[eq[:eq_a], endo[:a_t]] = m[:ρ_a]
    Ψ[eq[:eq_a],  exo[:a_sh]] = 1.

    # d
    Γ0[eq[:eq_d], endo[:d_t]] = 1.
    Γ1[eq[:eq_d], endo[:d_t]] = m[:ρ_d]
    Ψ[eq[:eq_d],  exo[:d_sh]] = 1.
#=
    # φ
    Γ0[eq[:eq_φ], endo[:φ_t]] = 1.
    Γ1[eq[:eq_φ], endo[:φ_t]] = m[:ρ_φ]
    Ψ[eq[:eq_φ],  exo[:φ_sh]] = 1.
=#
    # μ
    Γ0[eq[:eq_μ], endo[:μ_t]] = 1.
    Γ1[eq[:eq_μ], endo[:μ_t]] = m[:ρ_μ]
    Ψ[eq[:eq_μ],  exo[:μ_sh]] = 1.

    # b
    Γ0[eq[:eq_b], endo[:b_t]] = 1.
    Γ1[eq[:eq_b], endo[:b_t]] = m[:ρ_b]
    Ψ[eq[:eq_b],  exo[:b_sh]] = 1.

    # λ_p
    Γ0[eq[:eq_λ_p], endo[:λ_p_t]]   = 1.
    Γ1[eq[:eq_λ_p], endo[:λ_p_t]]   = m[:ρ_λ_p]
#=  Γ0[eq[:eq_λ_p], endo[:λ_p_t1]]  = -1
    Γ1[eq[:eq_λ_p], endo[:λ_p_t1]]  = -m[:η_λ_p]
    Γ0[eq[:eq_λ_p1], endo[:λ_p_t1]] = 1 =#
    Ψ[eq[:eq_λ_p], exo[:λ_p_sh]]   = 1.

    # λ_w
    Γ0[eq[:eq_λ_w], endo[:λ_w_t]] = 1.
    Γ1[eq[:eq_λ_w], endo[:λ_w_t]] = m[:ρ_λ_w]
#=  Γ0[eq[:eq_λ_w], endo[:λ_w_t1]]  = -1
    Γ1[eq[:eq_λ_w], endo[:λ_w_t1]]  = -m[:η_λ_w]
    Γ0[eq[:eq_λ_w1], endo[:λ_w_t1]] = 1 =#
    Ψ[eq[:eq_λ_w],  exo[:λ_w_sh]] = 1.

    # η_mp
    Γ0[eq[:eq_η_mp], endo[:η_mp_t]] = 1.
    Γ1[eq[:eq_η_mp], endo[:η_mp_t]] = m[:ρ_η_mp]
    Ψ[eq[:eq_η_mp],  exo[:η_mp_sh]] = 1.

    # g
    Γ0[eq[:eq_g], endo[:g_t]] = 1.
    Γ1[eq[:eq_g], endo[:g_t]] = m[:ρ_g]
    Ψ[eq[:eq_g],  exo[:g_sh]] = 1.

    # τ_S
    Γ0[eq[:eq_τ_S], endo[:τ_S_t]] = 1.
    Γ1[eq[:eq_τ_S], endo[:τ_S_t]] = m[:ρ_τ_S]
    Ψ[eq[:eq_τ_S],  exo[:τ_S_sh]] = 1.

    # τ_H1
    Γ0[eq[:eq_τ_H1], endo[:τ_H1_t]] = 1.
    Γ1[eq[:eq_τ_H1], endo[:τ_H1_t]] = m[:ρ_τ_H1]
    Ψ[eq[:eq_τ_H1],  exo[:τ_H1_sh]] = 1.

    # τ_H2
    Γ0[eq[:eq_τ_H2], endo[:τ_H2_t]] = 1.
    Γ1[eq[:eq_τ_H2], endo[:τ_H2_t]] = m[:ρ_τ_H2]
    Ψ[eq[:eq_τ_H2],  exo[:τ_H2_sh]] = 1.

    ### Expectational terms

    # π
    Γ0[eq[:eq_Eπ], endo[:π_t]]  = 1.
    Γ1[eq[:eq_Eπ], endo[:Eπ_t]] = 1.
    Π[eq[:eq_Eπ],  ex[:Eπ_sh]]  = 1.

    # π1
    Γ0[eq[:eq_Eπ1], endo[:π1_t]]  = 1.
    Γ1[eq[:eq_Eπ1], endo[:Eπ1_t]] = 1.
    Π[eq[:eq_Eπ1],  ex[:Eπ1_sh]]  = 1.

    # π2
    Γ0[eq[:eq_Eπ2], endo[:π2_t]]  = 1.
    Γ1[eq[:eq_Eπ2], endo[:Eπ2_t]] = 1.
    Π[eq[:eq_Eπ2],  ex[:Eπ2_sh]]  = 1.

    # λ_S
    Γ0[eq[:eq_Eλ_S], endo[:λ_S_t]]  = 1.
    Γ1[eq[:eq_Eλ_S], endo[:Eλ_S_t]] = 1.
    Π[eq[:eq_Eλ_S],  ex[:Eλ_S_sh]]  = 1.

    # λ_H1
    Γ0[eq[:eq_Eλ_H1], endo[:λ_H1_t]]  = 1.
    Γ1[eq[:eq_Eλ_H1], endo[:Eλ_H1_t]] = 1.
    Π[eq[:eq_Eλ_H1],  ex[:Eλ_H1_sh]]  = 1.

    # λ_H2
    Γ0[eq[:eq_Eλ_H2], endo[:λ_H2_t]]  = 1.
    Γ1[eq[:eq_Eλ_H2], endo[:Eλ_H2_t]] = 1.
    Π[eq[:eq_Eλ_H2],  ex[:Eλ_H2_sh]]  = 1.

    # x
    Γ0[eq[:eq_Ex], endo[:x_t]]  = 1
    Γ1[eq[:eq_Ex], endo[:Ex_t]] = 1
    Π[eq[:eq_Ex],  ex[:Ex_sh]]  = 1

    # ϕ
    Γ0[eq[:eq_Eϕ], endo[:ϕ_t]]  = 1.
    Γ1[eq[:eq_Eϕ], endo[:Eϕ_t]] = 1.
    Π[eq[:eq_Eϕ],  ex[:Eϕ_sh]]  = 1.

    # ρ
    Γ0[eq[:eq_Eρ], endo[:ρ_t]]  = 1.
    Γ1[eq[:eq_Eρ], endo[:Eρ_t]] = 1.
    Π[eq[:eq_Eρ],  ex[:Eρ_sh]]  = 1.

    # i_S
    Γ0[eq[:eq_Ei_S], endo[:i_S_t]]  = 1.
    Γ1[eq[:eq_Ei_S], endo[:Ei_S_t]] = 1.
    Π[eq[:eq_Ei_S],  ex[:Ei_S_sh]]  = 1.

    # w1
    Γ0[eq[:eq_Ew1], endo[:w1_t]]  = 1.
    Γ1[eq[:eq_Ew1], endo[:Ew1_t]] = 1.
    Π[eq[:eq_Ew1],  ex[:Ew1_sh]]  = 1.

    # w2
    Γ0[eq[:eq_Ew2], endo[:w2_t]]  = 1.
    Γ1[eq[:eq_Ew2], endo[:Ew2_t]] = 1.
    Π[eq[:eq_Ew2],  ex[:Ew2_sh]]  = 1.

#=    # c
    Γ0[eq[:eq_Ec], endo[:c_t]] = 1
    Γ1[eq[:eq_Ec], endo[:Ec_t]] = 1
    Π[eq[:eq_Ec], endo[:Ec_t]] = 1

    # c_H1
    Γ0[eq[:eq_Ec_H1], endo[:c_H1_t]] = 1
    Γ1[eq[:eq_Ec_H1], endo[:Ec_H1_t]] = 1
    Π[eq[:eq_Ec_H1], endo[:Ec_H1_t]] = 1

    # c_H2
    Γ0[eq[:eq_Ec_H2], endo[:c_H2_t]] = 1
    Γ1[eq[:eq_Ec_H2], endo[:Ec_H2_t]] = 1
    Π[eq[:eq_Ec_H2], endo[:Ec_H2_t]] = 1

    # c_H2
    Γ0[eq[:eq_Ec_H2], endo[:c_H2_t]] = 1
    Γ1[eq[:eq_Ec_H2], endo[:Ec_H2_t]] = 1
    Π[eq[:eq_Ec_H2], endo[:Ec_H2_t]] = 1
=#

    return Γ0, Γ1, C, Ψ, Π
end


text:
using ProgressMeter;
using UnicodePlots;

using AdalmPluto;

function bench_refill(radio::PlutoSDR, rate::Int)
    updateSamplingRate!(radio, rate; doLog=false);

    to_refill = 2 * rate;
    refilled = 0;
    start = time();
    while refilled < to_refill
        C_iio_buffer_refill(radio.rx.buf.C_ptr);
        refilled += radio.rx.buf.C_size;
    end
    stop = time();

    refill_rate = to_refill / (stop - start);

    return refill_rate;
end

function bench_refill_read(radio::PlutoSDR, rate::Int)
    updateSamplingRate!(radio, rate; doLog=false);

    to_read = 2 * rate;
    read = 0;
    start = time();
    while read < to_read
        C_iio_buffer_refill(radio.rx.buf.C_ptr);
        C_iio_channel_read!(radio.rx.iio.rx0_i, radio.rx.buf.C_ptr, radio.rx.buf.i_raw_samples);
        C_iio_channel_read!(radio.rx.iio.rx0_q, radio.rx.buf.C_ptr, radio.rx.buf.q_raw_samples);
        read += radio.rx.buf.C_size;
    end
    stop = time();

    reading_rate = to_read / (stop - start);

    return reading_rate;
end

function bench_refillJuliaBufferRX(radio::PlutoSDR, rate::Int)
    updateSamplingRate!(radio, rate; doLog=false);

    to_refill = 2 * rate;
    refilled = 0;
    start = time();
    while refilled < to_refill
        AdalmPluto.refillJuliaBufferRX(radio);
        refilled += radio.rx.buf.C_size;
    end
    stop = time();

    refill_rate = to_refill / (stop - start);

    return refill_rate;
end

function bench_recv(radio::PlutoSDR, samplingRate::Int)
    updateSamplingRate!(radio, samplingRate; doLog=false);

    # 2 seconds (if effective rate == radio rate)
    to_process = 2 * samplingRate;
    processed = 0;
    sig = zeros(ComplexF32, radio.rx.buf.C_size); # C_size is in samples
    start = time();
    while processed < to_process
        received = recv!(sig, radio);
        processed += received;
    end
    stop = time();

    radio_rate = radio.rx.effectiveSamplingRate;  # effective configuration (≠ effective rate)
    effective_rate = to_process / (stop - start);

    return radio_rate, effective_rate;
end

function bench_all(mode=:rx, doPlot=true)
    carrierFreq::Int = 100e6;
    gain = 64;
    bandwidth::Int = 20e6;

    rates           = Int.(6e5:1e6:5e7);
    radio_rates     = zeros(Int, length(rates));
    effective_rates = zeros(Float64, length(rates));

    radio = openPluto(carrierFreq, Int(1e6), gain; bandwidth=bandwidth);
    progress = Progress(length(rates), barglyphs=BarGlyphs("[=> ]"));
    if mode == :rx
        refill_rates    = zeros(Float64, length(rates));
        read_rates      = zeros(Float64, length(rates));
        refillJL_rates  = zeros(Float64, length(rates));
        for (idx, rate) in enumerate(rates)
            # benches
            refill_rate                = bench_refill(radio, rate);
            read_rate                  = bench_refill_read(radio, rate);
            refillJL_rate              = bench_refillJuliaBufferRX(radio, rate);
            radio_rate, effective_rate = bench_recv(radio, rate);
            # store values
            refill_rates[idx]    = refill_rate;
            read_rates[idx]      = read_rate;
            refillJL_rates[idx]  = refillJL_rate;
            radio_rates[idx]     = radio_rate;
            effective_rates[idx] = effective_rate;

            # progress bar
            update!(progress, idx);
        end
        res = Dict{Symbol, Array{Float64}}();
        res[:refill]    = refill_rates;
        res[:read]      = read_rates;
        res[:refillJL]  = refillJL_rates;
        res[:effective] = effective_rates;
        res[:cfg]       = radio_rates;
    elseif mode == :tx
        # placeholder
    end
    println("");

    close(radio);

    if doPlot
        plt = lineplot(rates, rates; width=2*length(rates), height=length(rates), name="desired");
        lineplot!(plt, rates, refill_rates; name="C_iio_buffer_refill");
        lineplot!(plt, rates, read_rates; name="refill + C_iio_channel_read");
        #  lineplot!(plt, rates, refillJL_rates; name="refill + read + reinterpret");
        lineplot!(plt, rates, radio_rates; name="effective cfg");
        lineplot!(plt, rates, effective_rates; name="effective rate (recv!)");
        return plt, res;
    end

    return res;
end

ext: jl


text:
module TestExamples2DCurved

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Curved Mesh" begin
  @testset "elixir_advection_basic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_curved.jl"),
      l2   = [9.14468177884088e-6],
      linf = [6.437440532947036e-5])
  end

  @testset "elixir_advection_extended_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended_curved.jl"),
      l2   = [4.842990962468553e-6],
      linf = [3.47372094784415e-5])
  end

  @testset "elixir_advection_extended_curved.jl with polydeg=4" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended_curved.jl"),
      l2   = [2.0610527374594057e-6],
      linf = [7.21793425673134e-6],
      atol = 1e-12, # required to make CI tests pass on macOS
      cells_per_dimension = (16, 23),
      polydeg = 4,
      cfl = 1.4)
  end

  @testset "elixir_advection_rotated.jl" begin
    @testset "elixir_advection_rotated.jl with α = 0.0" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_rotated.jl"),
        l2   = [7.013143474176369e-6],
        linf = [4.906526503622999e-5],
        alpha = 0.0)
    end

    @testset "elixir_advection_rotated.jl with α = 0.1" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_rotated.jl"),
        l2   = [7.013143474176369e-6],
        linf = [4.906526503622999e-5],
        alpha = 0.1)
    end

    @testset "elixir_advection_rotated.jl with α = 0.5 * pi" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_rotated.jl"),
        l2   = [7.013143474176369e-6],
        linf = [4.906526503622999e-5],
        alpha = 0.5 * pi)
    end
  end

  @testset "elixir_advection_parallelogram.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_parallelogram.jl"),
      l2   = [0.0005165995033861579],
      linf = [0.002506176163321605])
  end

  @testset "elixir_advection_waving_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_waving_flag.jl"),
      l2   = [0.00017007823031108628],
      linf = [0.0015264963372674245])
  end

  @testset "elixir_advection_free_stream_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_free_stream_curved.jl"),
      l2   = [5.984863383701255e-15],
      linf = [1.8207657603852567e-13])
  end

  @testset "elixir_advection_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic_curved.jl"),
      l2   = [0.00023766972629056245],
      linf = [0.004142508319267935])
  end

  @testset "elixir_advection_restart_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_curved.jl"),
      l2   = [6.398955192910044e-6],
      linf = [3.474337336717426e-5])
  end

  @testset "elixir_advection_restart_curved.jl with waving flag mesh" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_curved.jl"),
      l2   = [0.00017274040834067234],
      linf = [0.0015435741643734513],
      elixir_file="elixir_advection_waving_flag.jl",
      restart_file="restart_000041.h5")
  end

  @testset "elixir_advection_restart_curved.jl with free stream mesh" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_curved.jl"),
      l2   = [6.639325026542281e-15],
      linf = [1.829647544582258e-13],
      elixir_file="elixir_advection_free_stream_curved.jl",
      restart_file="restart_000068.h5")
  end

  @testset "elixir_euler_source_terms_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_curved.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5])
  end

  @testset "elixir_euler_source_terms_rotated.jl" begin
    @testset "elixir_euler_source_terms_rotated.jl with α = 0.0" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5],
      alpha = 0.0)
    end

    @testset "elixir_euler_source_terms_rotated.jl with α = 0.1" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"),
      l2   = [8.517816067144339e-7, 1.1545335192659009e-6, 1.3105743360575469e-6, 4.2778880887284874e-6],
      linf = [8.357837601113971e-6, 9.243793987812055e-6, 1.1305611891110345e-5, 4.496194865932779e-5],
      alpha = 0.1)
    end

    @testset "elixir_euler_source_terms_rotated.jl with α = 0.2 * pi" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"),
      l2   = [8.517916398386897e-7, 7.775434972297097e-7, 1.5639649498758719e-6, 4.277897693664215e-6],
      linf = [8.357487707666422e-6, 4.287882448716918e-6, 1.4423290043641401e-5, 4.496211536153538e-5],
      alpha = 0.2 * pi)
    end

    @testset "elixir_euler_source_terms_rotated.jl with α = 0.5 * pi" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5],
      alpha = 0.5 * pi)
    end
  end

  @testset "elixir_euler_source_terms_parallelogram.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_parallelogram.jl"),
      l2   = [1.0893247337495366e-5, 1.0425688625462763e-5, 1.755222105014883e-5, 5.136512290929069e-5],
      linf = [7.449791413693951e-5, 7.621073881791673e-5, 0.00011093303834863733, 0.00039625209916493986])
  end

  @testset "elixir_euler_source_terms_waving_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_waving_flag.jl"),
      l2   = [3.2280277059090205e-5, 3.481614479752735e-5, 2.8784017747658748e-5, 0.00011549476000734391],
      linf = [0.00025339087459608223, 0.0003425481056145152, 0.0002454647901921625, 0.0012806891514367535])
  end

  @testset "elixir_euler_free_stream_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_curved.jl"),
      l2   = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14],
      linf = [1.9539925233402755e-14, 2.9791447087035294e-13, 6.502853810985698e-13, 2.7000623958883807e-13])
  end

  @testset "elixir_euler_free_stream_curved.jl with FluxRotated(flux_lax_friedrichs)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_curved.jl"),
      surface_flux=FluxRotated(flux_lax_friedrichs),
      l2   = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14],
      linf = [1.9539925233402755e-14, 2.9791447087035294e-13, 6.502853810985698e-13, 2.7000623958883807e-13])
  end

  @testset "elixir_euler_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_nonperiodic_curved.jl"),
      l2   = [2.3653424742684444e-6, 2.1388875095440695e-6, 2.1388875095548492e-6, 6.010896863397195e-6],
      linf = [1.4080465931654018e-5, 1.7579850587257084e-5, 1.7579850592586155e-5, 5.956893531156027e-5])
  end
end

end # module


text:
module XAM

export
    SAM,
    BAM

include("sam/sam.jl")
include("bam/bam.jl")

using .SAM
using .BAM

end # module


text:
module TestNondeterministicThreading

include("preamble.jl")

function slow_identity(x)
    sleep(x)
    return x
end

throw_if_nothing(x) = x === nothing ? error("GOT NOTHING!") : x

@testset "happy path" begin
    @testset for seed in 1:100
        rng = MersenneTwister(seed)
        ntasks = rand(rng, 2:32)
        len = rand(rng, (1:10) .* ntasks)
        xs = 1e-3 .+ rand(rng, Float64, len) .* 9e-3
        ys =
            xs |>
            NondeterministicThreading(basesize = 1, ntasks = ntasks) |>
            Map(slow_identity)
        @test collect(ys) == xs
    end
end

@testset "Reduced" begin
    y =
        1:100 |>
        NondeterministicThreading(basesize = 1, ntasks = 10) |>
        Map() do i
            sleep(i * 1e-6)
            return i
        end |>
        ReduceIf() do i
            i == 1
        end |>
        foldxl(right)
    @test y == 1
end

@testset "error" begin
    @testset for seed in 1:100
        rng = MersenneTwister(seed)
        ntasks = rand(rng, 2:32)
        len = rand(rng, (1:10) .* ntasks)
        xs = 1e-3 .+ rand(rng, Float64, len) .* 9e-3
        xs = collect(Union{Float64,Nothing}, xs)
        xs[rand(rng, 1:len)] = nothing
        ys =
            xs |>
            NondeterministicThreading(basesize = 1, ntasks = ntasks) |>
            Map(throw_if_nothing) |>
            Map(slow_identity)
        err = @test_error collect(ys)
        @test occursin("GOT NOTHING!", sprint(showerror, err))
    end
end

end  # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -51.708688423784224 <= q <= 53.49505517188359)

                     add_NL_constraint(m, :(gelu(0.12651985253467002 + 0.8800990177661969*gelu(0.9367289958343452 + 0.010085610540752743*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.9059509477084742*gelu(-0.9076856947519341 + 0.6427398296015978*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.7662733112940292*gelu(0.014219790982957559 + -0.1830460372609095*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.3066544631255512*gelu(0.37873800689777815 + 0.3786798527526116*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.673637701918695*gelu(0.41028851710489267 + 0.3789873528571355*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + gelu(-0.8950204772394628 + -0.19603919659654778*gelu(0.9367289958343452 + 0.010085610540752743*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.902347490498435*gelu(-0.9076856947519341 + 0.6427398296015978*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.30440629677970144*gelu(0.014219790982957559 + -0.1830460372609095*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.6770808003403048*gelu(0.37873800689777815 + 0.3786798527526116*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.03232239144278681*gelu(0.41028851710489267 + 0.3789873528571355*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + gelu(-0.11076390433121786 + 0.7661103402200116*gelu(0.9367289958343452 + 0.010085610540752743*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.9009199905889229*gelu(-0.9076856947519341 + 0.6427398296015978*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.5631507797475948*gelu(0.014219790982957559 + -0.1830460372609095*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.4163802251888926*gelu(0.37873800689777815 + 0.3786798527526116*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.17662226846197893*gelu(0.41028851710489267 + 0.3789873528571355*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + gelu(0.8242960442017702 + -0.7506191216499181*gelu(0.9367289958343452 + 0.010085610540752743*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.31521888864028513*gelu(-0.9076856947519341 + 0.6427398296015978*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.5308637261158768*gelu(0.014219790982957559 + -0.1830460372609095*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.6654841763336146*gelu(0.37873800689777815 + 0.3786798527526116*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.2645412007156045*gelu(0.41028851710489267 + 0.3789873528571355*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + gelu(-0.6657453597374303 + 0.46323784626870257*gelu(0.9367289958343452 + 0.010085610540752743*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.7128966868012561*gelu(-0.9076856947519341 + 0.6427398296015978*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.43240821263487783*gelu(0.014219790982957559 + -0.1830460372609095*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.6998749584965918*gelu(0.37873800689777815 + 0.3786798527526116*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.6177298271304656*gelu(0.41028851710489267 + 0.3789873528571355*gelu(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*gelu(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*gelu(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*gelu(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*gelu(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Plain Float Type Unions
# -----------------------

"""
`PMTY{𝗽} = Union{Measurement{𝗽}} where 𝗽<:Union{Float16,Float32,Float64,BigFloat}`\n
𝗣lain 𝗠easurement 𝗧𝗬pe: Plain (unitless) `Measurement`s.
"""
PMTY{𝗽} = Union{Measurement{𝗽}} where 𝗽<:PREC

# United Float Type Unions
# ------------------------

"""
`UETY{𝗽,𝗱,𝘂} = Union{Quantity{𝗽,𝗱,𝘂}} where {𝗽<:PREC,𝗱,𝘂}`\n
𝗨nited 𝗘xact 𝗧𝗬pe: `PREC`ision-parametric united `Quantity`(ie)s.
"""
UETY{𝗽,𝗱,𝘂} = Union{Quantity{𝗽,𝗱,𝘂}} where {𝗽<:PREC,𝗱,𝘂}

"""
`UMTY{𝗽,𝗱,𝘂} = Union{Quantity{Measurement{𝗽},𝗱,𝘂}} where {𝗽<:PREC,𝗱,𝘂}`\n
𝗨nited 𝗠easurement 𝗧𝗬pe: `PREC`ision-parametric, `Measurement` united `Quantity`(ie)s.
"""
UMTY{𝗽,𝗱,𝘂} = Union{Quantity{Measurement{𝗽},𝗱,𝘂}} where {𝗽<:PREC,𝗱,𝘂}

"""
`UATY{𝗽,𝗱,𝘂} = Union{UETY{𝗽,𝗱,𝘂},UMTY{𝗽,𝗱,𝘂}} where {𝗽<:PREC,𝗱,𝘂}`
𝗨nited 𝗔mount 𝗧𝗬pe: `PREC`ision and `EXAC`tness-parametric, united `Quantity`(ie)s — the default
underlying data type for `EngTherm` `AMOUNTS`.
"""
UATY{𝗽,𝗱,𝘂} = Union{UETY{𝗽,𝗱,𝘂},UMTY{𝗽,𝗱,𝘂}} where {𝗽<:PREC,𝗱,𝘂}

#----------------------------------------------------------------------------------------------#
#                                        Auxiliar Types                                        #
#----------------------------------------------------------------------------------------------#

# The 4 type quadrants are:
#
#                  |    plain       united      |
#  ----------------+----------------------------+
#   float-based    |    plnF{𝗽}     UATY{𝗽,𝗱,𝘂} |
#   non-float reals|    REAL        uniR{𝘁,𝗱,𝘂} |

# REAL: plain Julia Reals other than `PREC` (since Unitful.Quantity <: Number)
REAL = Union{AbstractIrrational,Integer,Rational}

# Plain, unitless floats
plnF{𝗽} = Union{𝗽, Measurement{𝗽}} where 𝗽<:PREC

# Unit-ed reals
uniR{𝘁,𝗱,𝘂} = Union{Quantity{𝘁,𝗱,𝘂}, Quantity{Measurement{𝘁},𝗱,𝘂}} where {𝘁<:REAL,𝗱,𝘂}

#----------------------------------------------------------------------------------------------#
#                                   Raw Quantity Constructor                                   #
#----------------------------------------------------------------------------------------------#

# Dimensionless and Unitless constants
const DLESS = NoDims
const ULESS = Unitful.FreeUnits{(),NoDims,nothing}

# Adapted from https://github.com/PainterQubits/Unitful.jl/issues/283#issuecomment-552285299
_qty(x::𝗽) where 𝗽<:PREC = Quantity{𝗽, DLESS, ULESS}(x)
_qty(x::PMTY{𝗽}) where 𝗽<:PREC = Quantity{Measurement{𝗽}, DLESS, ULESS}(x)
_qty(x::UETY{𝗽,𝗱,𝘂}) where {𝗽<:PREC,𝗱,𝘂} = Quantity{𝗽,𝗱,𝘂}(x)
_qty(x::UMTY{𝗽,𝗱,𝘂}) where {𝗽<:PREC,𝗱,𝘂} = Quantity{Measurement{𝗽},𝗱,𝘂}(x)

#----------------------------------------------------------------------------------------------#
#                                        Promote Rules                                         #
#----------------------------------------------------------------------------------------------#

import Base: promote_rule

promote_rule(::Type{EX}, ::Type{EX}) = EX
promote_rule(::Type{EX}, ::Type{MM}) = MM
promote_rule(::Type{MM}, ::Type{EX}) = MM
promote_rule(::Type{MM}, ::Type{MM}) = MM


text:
module CommunalHelperConnectedMoveBlock
using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/ConnectedMoveBlock" ConnectedMoveBlock(
                            x::Integer, y::Integer,
                            width::Integer = Maple.defaultBlockWidth,
                            height::Integer = Maple.defaultBlockWidth,
                            direction::String="Right",
                            fast::Bool=false)

const placements = Ahorn.PlacementDict(
    "Connected Move Block ($direction) (Communal Helper)" => Ahorn.EntityPlacement(
        ConnectedMoveBlock,
        "rectangle",
        Dict{String, Any}(
            "direction" => direction
        )
    ) for direction in Maple.move_block_directions
)

Ahorn.editingOptions(entity::ConnectedMoveBlock) = Dict{String, Any}(
    "direction" => Maple.move_block_directions
)
Ahorn.minimumSize(entity::ConnectedMoveBlock) = 16, 16
Ahorn.resizable(entity::ConnectedMoveBlock) = true, true

Ahorn.selection(entity::ConnectedMoveBlock) = Ahorn.getEntityRectangle(entity)

const arrows = Dict{String, String}(
    "up" => "objects/moveBlock/arrow02",
    "left" => "objects/moveBlock/arrow04",
    "right" => "objects/moveBlock/arrow00",
    "down" => "objects/moveBlock/arrow06",
)
const buttonColor = (71, 64, 112, 255) ./ 255
const button = "objects/moveBlock/button"

const midColor = (4, 3, 23, 255) ./ 255
const highlightColor = (59, 50, 101, 255) ./ 255

# Gets rectangles from Solid Extensions
function getExtensionRectangles(room::Maple.Room)
    entities = filter(e -> e.name == "CommunalHelper/SolidExtension", room.entities)
    rects = []

    for e in entities
        push!(rects, Ahorn.Rectangle(
            Int(get(e.data, "x", 0)),
            Int(get(e.data, "y", 0)),
            Int(get(e.data, "width", 8)),
            Int(get(e.data, "height", 8))
        ))
    end

    return rects
end

# Checks for collision with an array of rectangles at specified tile position
function notAdjacent(entity::ConnectedMoveBlock, ox, oy, rects)
    x, y = Ahorn.position(entity)
    rect = Ahorn.Rectangle(x + ox + 4, y + oy + 4, 1, 1)

    for r in rects
        if Ahorn.checkCollision(r, rect)
            return false
        end
    end

    return true
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ConnectedMoveBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tileWidth = div(width, 8)
    tileHeight = div(height, 8)

    Ahorn.drawRectangle(ctx, 1, 1, width - 2, height - 2, highlightColor, highlightColor)
    Ahorn.drawRectangle(ctx, 8, 8, width - 16, height - 16, midColor)

    direction = lowercase(get(entity.data, "direction", "up"))
    arrowSprite = Ahorn.getSprite(arrows[lowercase(direction)], "Gameplay")

    block, innerCorners = "objects/moveBlock/base",  "objects/CommunalHelper/connectedMoveBlock/innerCorners"

    rects = getExtensionRectangles(room)
    rect = Ahorn.Rectangle(x, y, width, height)
	if !(rect in rects)
        push!(rects, rect)
    end

    for i in 1:tileWidth, j in 1:tileHeight
		drawX, drawY = (i - 1) * 8, (j - 1) * 8

        closedLeft = !notAdjacent(entity, drawX - 8, drawY, rects)
        closedRight = !notAdjacent(entity, drawX + 8, drawY, rects)
        closedUp = !notAdjacent(entity, drawX, drawY - 8, rects)
        closedDown = !notAdjacent(entity, drawX, drawY + 8, rects)
        completelyClosed = closedLeft && closedRight && closedUp && closedDown

        if completelyClosed

            if notAdjacent(entity, drawX + 8, drawY - 8, rects)
                # up right
                Ahorn.drawImage(ctx, innerCorners, drawX, drawY, 8, 0, 8, 8)

            elseif notAdjacent(entity, drawX - 8, drawY - 8, rects)
                # up left
                Ahorn.drawImage(ctx, innerCorners, drawX, drawY, 0, 0, 8, 8)

            elseif notAdjacent(entity, drawX + 8, drawY + 8, rects)
                # down right
                Ahorn.drawImage(ctx, innerCorners, drawX, drawY, 8, 8, 8, 8)

            elseif notAdjacent(entity, drawX - 8, drawY + 8, rects)
                # down left
                Ahorn.drawImage(ctx, innerCorners, drawX, drawY, 0, 8, 8, 8)

            else
                # entirely surrounded, fill tile
                Ahorn.drawImage(ctx, block, drawX, drawY, 8, 8, 8, 8)

            end
        else
            if closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, block, drawX, drawY, 8, 0, 8, 8)

            elseif closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, block, drawX, drawY, 8, 16, 8, 8)

            elseif closedLeft && !closedRight && closedUp && closedDown
                Ahorn.drawImage(ctx, block, drawX, drawY, 16, 8, 8, 8)

            elseif !closedLeft && closedRight && closedUp && closedDown
                Ahorn.drawImage(ctx, block, drawX, drawY, 0, 8, 8, 8)

            elseif closedLeft && !closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, block, drawX, drawY, 16, 0, 8, 8)

            elseif !closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, block, drawX, drawY, 0, 0, 8, 8)

            elseif !closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, block, drawX, drawY, 0, 16, 8, 8)

            elseif closedLeft && !closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, block, drawX, drawY, 16, 16, 8, 8)
            end
        end
    end

    Ahorn.drawRectangle(ctx, div(width - arrowSprite.width, 2) + 1, div(height - arrowSprite.height, 2) + 1, 8, 8, highlightColor, highlightColor)
    Ahorn.drawImage(ctx, arrowSprite, div(width - arrowSprite.width, 2), div(height - arrowSprite.height, 2))
end

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text: latexraw(i::Missings.Missing) = "\\textrm{NA}"


text:
using Documenter
using IntervalArithmetic, IntervalRootFinding

makedocs(
    modules = [IntervalRootFinding],
    doctest = true,
    format = :html,
    authors = "David P. Sanders",
    sitename = "IntervalRootFinding.jl",

    pages = Any[
        "Home" => "index.md",
        "API" => "api.md"
    ]
    )

# deploydocs(
#     deps = Deps.pip("pygments", "mkdocs", "mkdocs-cinder", "python-markdown-math"),
#     repo   = "github.com/dpsanders/IntervalConstraintProgramming.jl.git",
#     julia = "0.5",
#     osname = "linux"
# )

deploydocs(
    repo = "github.com/JuliaIntervals/IntervalRootFinding.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
    julia = "release",
    osname = "linux"
)


text:
import Base: @deprecate

@deprecate DataFrame!(args...; kwargs...) DataFrame(args...; copycols=false, kwargs...)

# TODO: remove these definitions in year 2021
by(args...; kwargs...) = throw(ArgumentError("by function was removed from DataFrames.jl. " *
                                             "Use the `combine(groupby(...), ...)` or `combine(f, groupby(...))` instead."))

aggregate(args...; kwargs...) = throw(ArgumentError("aggregate function was removed from DataFrames.jl. " *
                                                    "Use the `combine` function instead."))


text:
using DynamicAxisWarping
using MFCC
using Distances
using DataFrames
using Statistics

"""
    acdist(s1, s2; [method=:dtw, dist=SqEuclidean(), radius=10])

Calculate the acoustic distance between `s1` and `s2` with `method` version of dynamic time warping and `dist` as the interior distance function. Using `method=:dtw` uses vanilla dynamic time warping, while `method=:fastdtw` uses the fast dtw approximation. Note that this is not a true mathematical distance metric because dynamic time warping does not necessarily satisfy the triangle inequality, nor does it guarantee the identity of indiscernibles.

Args
=====

* `s1` Features-by-time array of first sound to compare
* `s2` Features-by-time array of second sound to compare
* `method` (keyword) Which method of dynamic time warping to use
* `dist` (keyword) Any distance function implementing the `SemiMetric` interface from the `Distances` package
* `dtwradius` (keyword) maximum warping radius for vanilla dynamic timew warping; if no value passed, no warping constraint is used argument unused when method=:fastdtw
* `fastradius` (keyword) The radius to use for the fast dtw method; argument unused when method=:dtw
"""
function acdist(s1, s2; method=:dtw, dist=SqEuclidean(), dtwradius=nothing, fastradius=10)

  if method == :dtw
    if isnothing(dtwradius)
      dtwradius = max(size(s1, 2), size(s2, 2))
    elseif size(s2, 2) < size(s1, 2)
      s1, s2 = s2, s1
    end
    imin, imax = radiuslimits(dtwradius, size(s1, 2), size(s2, 2))
    if last(imax) < size(s2, 2) imax[end] = size(s2, 2) end
    return dtw(s1, s2, dist, imin, imax)[1]
  elseif method == :fastdtw
    return fastdtw(s1, s2, dist, fastradius)[1]
  else
    error("Unsupported method argument.")
  end
end

"""
    acdist(s1::Sound, s2::Sound, rep=:mfcc; [method=:dtw, dist=SqEuclidean(), radius=10])

Convert `s1` and `s2` to a frequency representation specified by `rep`, then calculate acoustic distance between `s1` and `s2`. Currently only `:mfcc` is supported for `rep`, using defaults from the `MFCC` package except that the first coefficient for each frame is removed and replaced with the sum of the log energy of the filterbank in that frame, as is standard in ASR.
"""
function acdist(s1::Sound, s2::Sound, rep=:mfcc; method=:dtw, dist=SqEuclidean(), dtwradius=nothing, fastradius=10)

  if rep == :mfcc
    r1 = sound2mfcc(s1)
    r2 = sound2mfcc(s2)
  else
    error("Unsupported rep argument. Plesae consult documentation with ?acdist")
  end

  return acdist(r1, r2, method=method, dist=dist, dtwradius=dtwradius, fastradius=fastradius)
end

"""
    avgseq(S; [method=:dtw, dist=SqEuclidean(), radius=10, center=:medoid, dtwradius=nothing, progress=false])

Return a sequence representing the average of the sequences in `S` using the dba method for sequence averaging. Supports `method=:dtw` for vanilla dtw and `method=:fastdtw` for fast dtw approximation when performing the sequence comparisons. With `center=:medoid`, finds the medoid as the sequence to use as the initial center, and with `center=:rand` selects a random element in `S` as the initial center.

Args
======

* `S` An array of sequences to average
* `method` (keyword) The method of dynamic time warping to use
* `dist` (keyword) Any distance function implementing the `SemiMetric` interface from the `Distances` package
* `radius` (keyword) The radius to use for the fast dtw method; argument unused when method=:dtw
* `center` (keyword) The method used to select the initial center of the sequences in `S`
* `dtwradius` (keyword) How far a time step can be mapped when comparing sequences; passed directly to `DTW` function from `DynamicAxisWarping`; if set to `nothing`, the length of the longest sequence will be used, effectively removing the radius restriction
* `progress` Whether to show the progress coming from `dba`
"""
function avgseq(S; method=:dtw, dist=SqEuclidean(), fastradius=10, center=:medoid, dtwradius=nothing, progress=false)

  if isnothing(dtwradius)
    dtwradius = maximum(size(s, 2) for s in S)
  end

  if method == :dtw
    d = DTW(dtwradius, dist)
  elseif method == :fastdtw
    d = FastDTW(radius, dist)
  else
    error("Unsupported method argument.")
  end

  if center == :medoid
    D = zeros(length(S), length(S))
    for idx in CartesianIndices(D)
      D[idx] = acdist(S[idx[1]], S[idx[2]], method=method, dist=dist, dtwradius=dtwradius, fastradius=fastradius)[1]
    end
    sums = vec(sum(D, dims=1))
    init = S[argmin(sums)]
  elseif center == :rand
    init = rand(S)
  end

  return dba(S, d, init_center=init, show_progress=progress)[1]
end

"""
    avgseq(S::Array{Sound}, rep=:mfcc; [method=:dtw, dist=SqEuclidean(), radius=10, center=:medoid, dtwradius=nothing, progress=false])

Convert the `Sound` objects in `S` to a representation designated by `rep`, then find the average sequence of them. Currently only `:mfcc` is supported for `rep`, using defaults from the `MFCC` package except that the first coefficient for each frame is removed and replaced with the sum of the log energy of the filterbank in that frame, as is standard in ASR.
"""
function avgseq(S::Array{Sound}, rep=:mfcc; method=:dtw, dist=SqEuclidean(), fastradius=10, center=:medoid, dtwradius=nothing, progress=false)
  if rep == :mfcc
    S = sound2mfcc.(S)
  else
    error("Unsupported rep argument. Please consult documentation with ?acdist")
  end

  return avgseq(S, method=method, dist=dist, fastradius=fastradius, center=center, dtwradius=dtwradius, progress=false)
end

"""
    distinctiveness(s, corpus; [method=:dtw, dist=SqEuclidean(), radius=10, reduction=mean])

Calculates the acoustic distinctiveness of `s` given the corpus `corpus`. The `method`, `dist`, and `radius` arguments are passed into `acdist`. The `reduction` argument can be any function that reduces an iterable to one number, such as `mean`, `sum`, or `median`.

For more information, see Kelley (2018, September, How acoustic distinctiveness affects spoken word recognition: A pilot study, DOI: 10.7939/R39G5GV9Q) and Kelley & Tucker (2018, Using acoustic distance to quantify lexical competition, DOI: 10.7939/r3-wbhs-kr84).
  """
function distinctiveness(s, corpus; method=:dtw, dist=SqEuclidean(), dtwradius=nothing, fastradius=10, reduction=mean)
  return reduction(map(x -> acdist(s, x, method=method, dist=dist, dtwradius=dtwradius, fastradius=fastradius), corpus))
end

"""
    distinctiveness(s::Sound, corpus::Array{Sound}, rep=:mfcc; [method=:dtw, dist=SqEuclidean(), radius=10, reduction=mean])

Converts `s` and `corpus` to a representation specified by `rep`, then calculates the acoustic distinctiveness of `s` given `corpus`. Currently only `:mfcc` is supported for `rep`, using defaults from the `MFCC` package except that the first coefficient for each frame is removed and replaced with the sum of the log energy of the filterbank in that frame, as is standard in ASR.
"""
function distinctiveness(s::Sound, corpus::Array{Sound}, rep=:mfcc; method=:dtw, dist=SqEuclidean(), dtwradius=nothing, fastradius=10, reduction=mean)

  if rep == :mfcc
    s = sound2mfcc(s)
    corpus = sound2mfcc.(corpus)
  else
    error("Unsupported rep argument")
  end

  return distinctiveness(s, corpus, method=method, dist=dist, dtwradius=dtwradius, fastradius=fastradius, reduction=reduction)
end

function sound2mfcc(s::Sound; useFrameEngery=true, kw...)
  m, e, _ = mfcc(s.samples, s.sr; kw...)
  if useFrameEngery
    e = log.(sum(e, dims=2))
    m[:,1] = e
  end
  return Array(m')
end


text:
using HMCExamples

HMCExamples.main_rbc_1_kalman(ARGS)

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
using Gadfly

ydata=[19, 19, 19, 50, 50, 50, 49, 68, 68, 28, 28, 28, 49, 49, 49, 16, 16, 48, 48, 48, 28, 28, 28, 21, 16, 22, 15, 41, 41, 31, 31, 21, 21, 21, 11, 11, 16, 17, 14, 53, 15, 15, 15, 46, 21, 21, 21, 41, 22, 22, 22, 38, 38, 38, 38, 16, 16, 16, 16, 47, 47, 21, 21, 23, 23, 47, 47, 47, 15, 15, 15, 26, 14, 46, 46, 46, 16, 16, 16, 48, 36, 47, 11, 11, 45, 45, 45, 45, 5]

plot(x=ones(length(ydata)), y=ydata, Geom.boxplot)


text:
using GridInterpolations
using Base.Test

# use import otherwise Interpolations.interpolate conflicts with GridInterpolations.interpolate
import Interpolations: BSpline, Linear, OnGrid, Flat
import Interpolations

function compareToGrid(testType::Symbol=:random, numDims::Int=3, pointsPerDim::Int=5, testPointsDim::Int=5; numRandomTests::Int=1000, eps::Float64=1e-10)
    # Compare interpolation results from GridInterpolations with the Grid Module:
    # create a multi-dimensional test Grid with data, cut points, etc.
    # (dataM, gridM, gridI) = compareToGrid(2, 3, 5);
    cutsDim = pointsPerDim*ones(Int,numDims)
    gridData = rand(cutsDim...)

    # Set up the data structures:
    gridI = Interpolations.interpolate(gridData, BSpline(Linear()), OnGrid())  		# This is for the Interpolations package
    gridI = Interpolations.extrapolate(gridI, Flat())
    cutPointMat = 0:(1/(pointsPerDim-1)):1
    cutPoints = [cutPointMat for i=1:numDims]
    gridM = RectangleGrid(tuple(cutPoints...)...) 				# Mykel's interpolation package

    # Package the data in a way that can be read by the interpolation package:
    dataM = Array{Float64}(length(gridM))
    for i=1:length(gridM)
        # x = ind2x(gridM, i)
        x = ind2x(gridM, i)
        dataM[i] = gridData[getFractionalIndexes(gridM,x)...]
    end

    # Check a finer grid of coordinates, for comprehensiveness?

    if testType==:random
        # Interpolate in GridInterpolations with the same points, report any discrepancy
        # Select a number of random points and compare the interpolation results.
        testPoint = zeros(numDims)
        for i=1:numRandomTests
            rand!(testPoint)
            testI = gridI[((pointsPerDim-1)*testPoint+1)...]
            testM = GridInterpolations.interpolate(gridM,dataM,testPoint)

            if (abs(testI-testM)>eps)
                display("Failed Random Point Interpolation Test")
                display("Non-matching interpolation points:")
                display("Test point = $testPoint")
                display("Grid Package interpolation result = $testI")
                display("Interpolation Module result = $testM")
                return false
            end

        end
        return true
    end

    if testType==:extrapNeg
        # Test behavior outside the grid, i.e. extrapolation:
        testPoint = zeros(numDims)
        for i=1:numRandomTests
            testPoint = -rand!(testPoint)
            testI = gridI[((pointsPerDim-1)*testPoint+1)...]
            testM = interpolate(gridM,dataM,testPoint)

            if (abs(testI-testM)>eps)
                display("Failed Random Negative Point Extrapolation Test")
                display("Non-matching interpolation points:")
                display("Test point = $testPoint")
                display("Grid Package interpolation result = $testI")
                display("Interpolation Module result = $testM")
                return false
            end

        end
        return true
    end

    if testType==:extrapPos
        for i=1:numRandomTests
            testPoint = rand(numDims)+1
            testI = gridI[((pointsPerDim-1)*testPoint+1)...]
            testM = interpolate(gridM,dataM,testPoint)

            if (abs(testI-testM)>eps)
                display("Failed Random Positive Point Extrapolation Test")
                display("Non-matching interpolation points:")
                display("Test point = $testPoint")
                display("Grid Package interpolation result = $testI")
                display("Interpolation Module result = $testM")
                return false
            end

        end
        return true
    end

    return true

end

function getFractionalIndexes(g::AbstractGrid, s::Array)
    # Returns the fractional index of sprime within the grid-defined discretization.

    fracInd = Array{Int64}(length(g.cutPoints))

    for i=1:length(g.cutPoints)
        gridDisc = g.cutPoints[i]
        fracInd[i] = getFracIndex(gridDisc, s[i])
    end

    return fracInd

end

function getFracIndex(vararray::Array, value::Float64)
    # Searches through the array vararray and returns the fractional index position
    # that value lies in.  It CAPassumes vararray increases monotonically, and returns
    # the first or last element if the value is outside the bounds of the array.

    if value <= vararray[1]
        return 1
    end

    if value >= vararray[end]
        return indmax(vararray)
    end

    i=1
    while value > vararray[i+1]
        i+=1
    end

    return i+(value-vararray[i])÷(vararray[i+1] - vararray[i])

end

function simplexMagic(NDISC::Int=20, NPOINTS::Int=3, checkFileName::AbstractString="", eps::Float64=1e-10)
    if isempty(checkFileName)
        checkFileName = joinpath(@__DIR__, "simplexMagicTest20.txt")
    end

    val = transpose([ [8.,1,6] [3,5,7] [4,9,2] ]) # transposed magic(3) from matlab

    sGrid = SimplexGrid(1.:NPOINTS, 1.:NPOINTS)
    sInterpVal = zeros(NDISC, NDISC)

    for i = 1:NDISC, j = 1:NDISC
        valX = mapPt(i, NPOINTS, NDISC)
        valY = mapPt(j, NPOINTS, NDISC)
        sInterpVal[i, j] = interpolate(sGrid, vcat(val...), [valX, valY])
    end

    sInterpValTest = readdlm(checkFileName)

    testErr = sum(abs(sInterpVal-sInterpValTest))

    if (testErr > eps)
        display("Failed Simplex Comparison Test")
        display("Test error = $testErr")
    end

    return testErr<eps

end

# constructs a new grid from repr output and tests to see if it's the same
function reprConstruct()
    g1 = SimplexGrid([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
    g2 = eval(parse(repr(g1)))
    return g1.cutPoints == g2.cutPoints
end

function mapPt(pt::Int, NPOINTS::Int, NDISC::Int)
    return pt * (NPOINTS-1) / NDISC + 1
end

function checkCounters()
    dims = 2
    nPts = 100
    grid = RectangleGrid(1.:10, 1.:10)
    if length(grid) != nPts || dimensions(grid) != dims
        return false
    end
    data = rand(10,10)
    temp = interpolate(grid, data, [-1.5,1.5])
    temp = interpolate(grid, data, [1.,1.])
    grid = SimplexGrid(1.:10, 1.:10)
    if length(grid) != nPts || dimensions(grid) != dims
        return false
    end
    temp = interpolate(grid, data, [-1.5,1.5])
    temp = interpolate(grid, data, [1.,1.])
    temp = interpolate(grid, data, [1.5,1.5])
#    showcompact(grid)
#    show(grid)
    return true
end

function testMask()
    x = [1.5, 1.5]
    vals = collect(1.:9)
    mask = falses(9)
    grid = RectangleGrid(1.:3, 1.:3)
    interped = interpolate(grid, vals, x)
    masked   = maskedInterpolate(grid, vals, x, mask)
    @test interped ≈ masked atol = 1e-3
    mask[1] = true # mask the first entry
    truth  = 3.6666
    masked = maskedInterpolate(grid, vals, x, mask)
    @test truth ≈ masked atol = 1e-3
    return true
end

@test compareToGrid(:random) == true
@test compareToGrid(:extrapPos) == true
#@test compareToGrid(:extrapNeg)==true

@test simplexMagic() == true

@test reprConstruct() == true

@test checkCounters() == true

@test testMask() == true

# Check that duplicates are disallowed
@test_throws ErrorException RectangleGrid([1,1])
@test_throws ErrorException SimplexGrid([1,1])

# Sanity check overall implementations
function test_rect_implemented()
    rgrid = RectangleGrid([2,5],[2,5])

	@test length(rgrid) == 4
	@test GridInterpolations.dimensions(rgrid) == 2

	@test ind2x(rgrid,1) == [2,2]
	@test ind2x(rgrid,2) == [5,2]
	@test ind2x(rgrid,3) == [2,5]
	@test ind2x(rgrid,4) == [5,5]
	x = [0,0]
	ind2x!(rgrid, 4, x)
	@test x == [5,5]

	@test interpolants(rgrid, [1,1]) == ([1],[1.0])
	@test interpolants(rgrid, [2,5]) == ([3],[1.0])

	indices, weights = interpolants(rgrid, [1.5,3])
	@test indices == [1,3]
	@test isapprox(weights, [0.666667, 0.333333], rtol=1e-5)

	@test interpolate(rgrid, [1,2,3,4], [1,1]) == 1.0
	@test maskedInterpolate(rgrid, [1,2,3,4], [1,1], BitArray([false, false, false, false])) == 1.0
	@test isnan(maskedInterpolate(rgrid, [1,2,3,4], [1,1], BitArray([true, false, false, false])))

	@test isapprox(interpolate(rgrid, [1,2,3,4], [1.5,3]), 1.66666666, rtol=1e-5)
	@test maskedInterpolate(rgrid, [1,2,3,4], [1.5,3], BitArray([true, false, false, false])) == 3
	@test maskedInterpolate(rgrid, [1,2,3,4], [1.5,3], BitArray([false, false, true, false])) == 1
end
test_rect_implemented()

# Sanity check overall implementations
function test_simplex_implemented()
    grid = SimplexGrid([2,5],[2,5])

	@test length(grid) == 4
	@test GridInterpolations.dimensions(grid) == 2

	@test ind2x(grid,1) == [2,2]
	@test ind2x(grid,2) == [5,2]
	@test ind2x(grid,3) == [2,5]
	@test ind2x(grid,4) == [5,5]
	x = [0,0]
	ind2x!(grid, 4, x)
	@test x == [5,5]

	indices, weights = interpolants(grid, [1,1])
	full_weight_indices = find(x -> x == 1, weights)
	@test  length(full_weight_indices) == 1
	@test  weights[full_weight_indices[1]] == 1.0

	indices, weights = interpolants(grid, [1.5,3])
	full_weight_indices = find(x -> isapprox(x, 0.666667, rtol=1e-5), weights)
	@test  length(full_weight_indices) == 1
	@test  isapprox(weights[full_weight_indices[1]], 0.666667, rtol=1e-5) == true
	full_weight_indices = find(x -> isapprox(x, 0.333333, rtol=1e-5), weights)
	@test  length(full_weight_indices) == 1
	@test  isapprox(weights[full_weight_indices[1]], 0.333333, rtol=1e-5)  == true

	@test interpolate(grid, [1,2,3,4], [1,1]) == 1.0
	@test maskedInterpolate(grid, [1,2,3,4], [1,1], BitArray([false, false, false, false])) == 1.0
	@test isnan(maskedInterpolate(grid, [1,2,3,4], [1,1], BitArray([true, false, false, false])))

	@test isapprox(interpolate(grid, [1,2,3,4], [1.5,3]), 1.66666666, rtol=1e-5)
	@test maskedInterpolate(grid, [1,2,3,4], [1.5,3], BitArray([true, false, false, false])) == 3
	@test maskedInterpolate(grid, [1,2,3,4], [1.5,3], BitArray([false, false, true, false])) == 1
end
test_simplex_implemented()

# check whether rectangle & simplex expect sorted order of cutpoints on dimensions
function test_ordering(grid)
	@test ind2x(grid,1) == [2,18] # 1
	@test ind2x(grid,2) == [5,18] # 2
	@test ind2x(grid,3) == [2,15] # 3
	@test ind2x(grid,4) == [5,15] # 4
	@test ind2x(grid,5) == [2,12] # 5
	@test ind2x(grid,6) == [5,12] # 6

	return interpolate(grid, [1,2,3,4,5,6], [6, 20]) == 2.0
end
@test_throws ErrorException test_ordering( RectangleGrid([2,5], [18,15,12]) )
@test_throws ErrorException test_ordering( SimplexGrid([2,5], [18,15,12]) )

include(joinpath(@__DIR__, "..", "src", "interpBenchmarks.jl"))
compareBenchmarks(4, 10, 100, quiet=false)

println("All tests complete")


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ 47554e79-497d-411d-ab2f-182cc915af48
begin
	filter!(path -> path != "@v#.#", LOAD_PATH)
	using Pkg
	if occursin("/srv/julia/pkg/pluto_notebooks", @__DIR__)
		Pkg.develop(url = "https://github.com//likanzhan/StatisticsWithJulia.jl")
		Pkg.activate(joinpath(Pkg.devdir(), "StatisticsWithJulia.jl"))
	else
		Pkg.activate(Base.current_project())
	end
	Pkg.instantiate()
	using StatisticsWithJulia
end

# ╔═╡ a888c516-dba1-4cda-9ada-2e516fc106e9
StatisticsWithJulia.header(title = "Julia 语言与数据科学初步", author = "战立侃")

# ╔═╡ 29bad371-6dd9-4ec1-8f9e-95b034a07a21
md"""

- [期末试卷](https://raw.githubusercontent.com/likanzhan/StatisticsWithJulia.jl/main/notebooks/期末考试试题.pdf)

Julia 语言是一门面向未来的技术型编程语言。 本课程将以 [Julia语言](https://www.julialang.org) 为基本工具， 介绍统计和数据科学领域的一些基本知识。 本课程的主要参考书为：

- Nazarathy, Y., & Klok, H. (2021). *Statistics with Julia: Fundamentals for data science, machine learning and artificial intelligence*. Springer.

想要了解更多 Julia 语言本身的一般知识， 请参看 Julia 语言 的[官方文档](https://docs.julialang.org)。
"""

# ╔═╡ 82b7a596-1b6c-11ec-1731-5b8c4d67cc2c
StatisticsWithJulia.list_notebooks(@__FILE__, "notebooks")

# ╔═╡ 31b32b9b-326e-4bd4-9ca2-fbd249394a87
md"""
要在本地运行该课程，请在 Julia REPL 下执行以下代码：

```julia
julia> using Pkg
julia> Pkg.develop(url = "https://github.com//likanzhan/StatisticsWithJulia.jl")
julia> Pkg.activate(joinpath(Pkg.devdir(), "StatisticsWithJulia.jl"))
julia> Pkg.instantiate()
julia> using StatisticsWithJulia
julia> StatisticsWithJulia.start()
```

"""

# ╔═╡ Cell order:
# ╟─a888c516-dba1-4cda-9ada-2e516fc106e9
# ╟─47554e79-497d-411d-ab2f-182cc915af48
# ╟─29bad371-6dd9-4ec1-8f9e-95b034a07a21
# ╟─82b7a596-1b6c-11ec-1731-5b8c4d67cc2c
# ╟─31b32b9b-326e-4bd4-9ca2-fbd249394a87