text:
"""
```
COTHANK{T} <: AbstractRepModel{T}
```

The `THANK` type defines the structure of THANK, an implementation
of the DSGE model in Bilbiie, Primiceri, Tambalotti (2019)

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed
  as a function of elements of `parameters`.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for all model
  parameters and steady-states to their indices in `parameters` and
  `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

* `endogenous_states::OrderedDict{Symbol,Int}`: Maps each state to a column in
  the measurement and equilibrium condition matrices.

* `exogenous_shocks::OrderedDict{Symbol,Int}`: Maps each shock to a column in
  the measurement and equilibrium condition matrices.

* `expected_shocks::OrderedDict{Symbol,Int}`: Maps each expected shock to a
  column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::OrderedDict{Symbol,Int}`: Maps each equlibrium
  condition to a row in the model's equilibrium condition matrices.

* `endogenous_states_augmented::OrderedDict{Symbol,Int}`: Maps lagged states to
  their columns in the measurement and equilibrium condition equations. These
  are added after `gensys` solves the model.

* `observables::OrderedDict{Symbol,Int}`: Maps each observable to a row in the
  model's measurement equation matrices.

* `pseudo_observables::OrderedDict{Symbol,Int}`: Maps each pseudo-observable to
  a row in the model's pseudo-measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"m1002\", cached here for
  filepath computation.

* `subspec::String`: The model subspecification number, indicating that
  some parameters from the original model spec (\"ss10\") are initialized
  differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation
  without changing the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model
  units. DSGE.jl will fetch data from the Federal Reserve Bank of St. Louis's
  FRED database; all other data must be downloaded by the user. See `load_data`
  and `Observable` for further details.

* `pseudo_observable_mappings::OrderedDict{Symbol,PseudoObservable}`: A
  dictionary that stores names and transformations to/from model units. See
  `PseudoObservable` for further details.
"""
mutable struct COTHANK{T} <: AbstractRepModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                       # model steady-state values
    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
                                                           # parameters and steady-states

    endogenous_states::OrderedDict{Symbol,Int}             # these fields used to create matrices in the
    exogenous_shocks::OrderedDict{Symbol,Int}              # measurement and equilibrium condition equations.
    expected_shocks::OrderedDict{Symbol,Int}               #
    equilibrium_conditions::OrderedDict{Symbol,Int}        #
    endogenous_states_augmented::OrderedDict{Symbol,Int}   #
    observables::OrderedDict{Symbol,Int}                   #
    pseudo_observables::OrderedDict{Symbol,Int}            #

    spec::String                                           # Model specification number (eg "m990")
    subspec::String                                        # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                         # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                    # Settings/flags for testing mode
    rng::MersenneTwister                                   # Random number generator
    testing::Bool                                          # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
    pseudo_observable_mappings::OrderedDict{Symbol, PseudoObservable}
end

description(m::COTHANK) = "Bilbiie, Primiceri, Tambalotti (2019), $(m.subspec)."

"""
`init_model_indices!(m::COTHANK)`

Arguments:
`m:: COTHANK`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::COTHANK)
    # Endogenous states
    endogenous_states = [:y_t, :y1_t, :y2_t, :k_t, :k_S_t, :L1_t,  :L2_t, :klR1_t, :klR2_t,   # check if missing ι terms
                         :w1_t, :w2_t, :π_t, :π1_t, :π2_t, :mc1_t, :mc2_t, :λ_p_t, :λ_w_t,      # sticky prices
                         :λ_S_t, :λ_H1_t, :λ_H2_t, :λ_SH1_t, :λ_SH2_t, :c_t,  :c_H1_t, :c_H2_t, :c_S_t, :i_S_t,
                         :R_t, :u_t, :ϕ_t, :z_t, :ρ_t, :a_t, :d_t, :x_t, :b_t, :t_t, #:φ_t,
                         :g_t, :g_w1_t, :g_w2_t, :μ_t, :bᴿ_t, :η_mp_t, :t_S_t, :t_H1_t, :t_H2_t,
                         :τ_t, :τ_S_t, :τ_H1_t, :τ_H2_t, :π_S_t, :price1_t, :price2_t, :Eπ1_t, :Eπ2_t, :Eλ_S_t,
                         :Eλ_H1_t, :Eλ_H2_t, :Ex_t, :Eϕ_t, :Eρ_t, :Ei_S_t, :Ew1_t, :Eπ_t, :Ew2_t]
                                                                                # CHECK BACK IF ANY CAN BE AUGMENTED VARIABLES
#=                         :y_f_t, :k_f_t, :L_f_t, :Rk_f_t, :w_f_t, :mc_f_t,       # flexible prices
                         :λ_f_t, :c_f_t, :R_f_t, :u_f_t, :ϕ_f_t,
                         :i_f_t, :kbar_f_t, :wgap_f_t, :gdp_f_t, :Ec_f_t, :Eλ_f_t,
                         :Eϕ_f_t, :ERk_f_t, :Ei_f_t, :λ_h_f_t, :λ_s_f_t, :c_h_f_t, :c_s_f_t,
                         :t_h_f_t, :Ey_f_t, :Eλ_h_f_t, :Eλ_s_f_t]
=#
    # Exogenous shocks
    exogenous_shocks = [:a_sh, :d_sh, :λ_p_sh, :b_sh, :z_sh, :τ_H1_sh, :τ_H2_sh, :μ_sh, :λ_w_sh,
                        :τ_S_sh, :g_sh, :η_mp_sh]

    # Expectations shocks
    expected_shocks = [:Eπ1_sh, :Eπ2_sh, :Eλ_S_sh, :Eλ_H1_sh, :Eλ_H2_sh,
                       :Ex_sh, :Eϕ_sh, :Eρ_sh, :Ei_S_sh, :Ew1_sh, :Eπ_sh, :Ew2_sh] # sticky prices

    # Equilibrium conditions
    equilibrium_conditions = [:eq_KL_ratio_good1, :eq_mc_good1, :eq_KL_ratio_good2, :eq_mc_good2, :eq_production_good1, :eq_production_good2,
                              :eq_production_finalgood, :eq_price_index, :eq_demand_good1, :eq_π1, :eq_π2, :eq_price_phillips_good1,
                              :eq_price_phillips_good2, :eq_marg_utility_S, :eq_marg_utility_H1, :eq_marg_utility_H2, :eq_average_c,
                              :eq_euler, :eq_c_H1, :eq_c_H2, :eq_cap_util, :eq_ϕ, :eq_i, :eq_cap_input, :eq_cap_accum, :eq_wage_phillips_good1,
                              :eq_marg_sub_good1, :eq_avg_marg_util_SH1, :eq_wage_phillips_good2, :eq_marg_sub_good2, :eq_avg_marg_util_SH2,
                              :eq_mp, :eq_revenue, :eq_transfer, :eq_realdebt, :eq_fiscal_rule, :eq_redistribute_H1, :eq_redistribute_H2,
                              :eq_real_profits,
                              :eq_cap_market_clear, :eq_goods_market_clear, :eq_gdp,
                              :eq_a, :eq_d, :eq_λ_p, :eq_b, :eq_z, :eq_τ_H1, :eq_τ_H2, :eq_μ, :eq_λ_w, :eq_τ_S, :eq_g, :eq_η_mp, #eq_φ
                               :eq_Eπ1, :eq_Eπ2, :eq_Eλ_S, :eq_Eλ_H1, :eq_Eλ_H2, :eq_Ex, :eq_Eϕ,
                              :eq_Eρ, :eq_Ei_S, :eq_Ew1, :eq_Eπ, :eq_Ew2]
#=
                              :eq_y_f, :eq_k_f, :eq_L_f, :eq_Rk_f, :eq_w_f, :eq_mc_f,       # flexible prices
                              :eq_λ_f, :eq_c_f, :eq_R_f, :eq_u_f, :eq_ϕ_f,
                              :eq_i_f, :eq_kbar_f, :eq_wgap_f, :eq_gdp_f, :eq_Ec_f, :eq_Eλ_f,
                              :eq_Eϕ_f, :eq_ERk_f, :eq_Ei_f, :eq_λ_h_f, :eq_λ_s_f, :eq_c_h_f,
                              :eq_c_s_f, :eq_t_h_f, :eq_Ey_f, :eq_Eλ_h_f, :eq_Eλ_s_f]
=#
    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = []

    # Observables
    observables = keys(m.observable_mappings)

    # Pseudo-observables
    pseudo_observables = keys(m.pseudo_observable_mappings)

    for (i,k) in enumerate(endogenous_states);           m.endogenous_states[k]           = i end
    for (i,k) in enumerate(exogenous_shocks);            m.exogenous_shocks[k]            = i end
    for (i,k) in enumerate(expected_shocks);             m.expected_shocks[k]             = i end
    for (i,k) in enumerate(equilibrium_conditions);      m.equilibrium_conditions[k]      = i end
    for (i,k) in enumerate(endogenous_states);           m.endogenous_states[k]           = i end
    for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i+length(endogenous_states) end
    for (i,k) in enumerate(observables);                 m.observables[k]                 = i end
    for (i,k) in enumerate(pseudo_observables);          m.pseudo_observables[k]          = i end
end

function COTHANK(subspec::String="ss1";
                   custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                   testing = false)

    # Model-specific specifications
    spec               = split(basename(@__FILE__),'.')[1]
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # Initialize empty model
    m = COTHANK{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(), OrderedDict{Symbol,Int}(),

            # model indices
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}(),
            OrderedDict{Symbol,PseudoObservable}())

    # Set settings
    model_settings!(m)
    DSGE.default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Set observable and pseudo-observable transformations
    init_observable_mappings!(m)
    #init_pseudo_observable_mappings!(m)

    # Initialize parameters
    init_parameters!(m)

    init_model_indices!(m)
    init_subspec!(m)
    steadystate!(m)

    return m
end

"""
```
init_parameters!(m::COTHANK)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::COTHANK)
    # Calibrated parameters#
    m <= parameter(:δ, 0.025, fixed=true,
                   description="δ: The capital depreciation rate.",
                   tex_label = "\\delta" ) # from Giorgio's matlab script

    # Non standard devation parameters
    m <= parameter(:α, 1-0.6/(1.15), (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.05), fixed = false,
                   description="α: Capital elasticity in the intermediate goods sector's production function (also known as the capital share).",
                   tex_label = "\\alpha") # from Giorgio's matlab script

    m <= parameter(:ι_p, 0.21, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.5, 0.15), fixed = false,
                   description = "ι_p: The weight attributed to last period's inflation in price indexation. " *
                   "(1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label = "\\iota_p") # from THANK

    m <= parameter(:ι_w, 0.15, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.5, 0.15), fixed = false,
                   description="ι_w: The weight attributed to last period's wage in wage indexation. "
                    * "(1-ι_w) is the weight attributed to steady-state wages.",
                   tex_label = "\\iota_w") # from THANK

    m <= parameter(:h, 0.5, (-5.0, 5.0), (-5., 5.), ModelConstructors.Untransformed(), Beta(0.5, 0.1), fixed = false,
                   description="h: habit formation parameter.",
                   tex_label = "h") # from Giorgio's matlab script

    m <= parameter(:λ_p_ss, 0.15,(-1000.,1000.),(-1000., 1000.), ModelConstructors.Untransformed(), Normal(0.15, 0.05),  fixed = false,
                   description="λ_p_ss: The steady state net price markup.",
                   tex_label = "\\lambda_p") # from Giorgio's matlab script

    m <= parameter(:λ_w_ss, 0.12,(-1000.,1000.),(-1000., 1000.), ModelConstructors.Untransformed(), Normal(0.15, 0.05), fixed = false,
                   description = "λ_w_ss: The steady state net wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label = "\\lambda_w") # from Giorgio's matlab script

    m <= parameter(:L_ss, 0.5, (-1000., 1000.), (-1e3, 1e3), ModelConstructors.Untransformed(), Normal(0.0,0.5), fixed=false,
                   description="L_ss: The steady state for log hours.",
                   tex_label = "\\log(L)^{ss}") # from Giogio's matlab script

    m <= parameter(:Fβ, 0.1, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Gamma(0.25, 0.1), fixed=false,
                   description = "Fβ: Discount rate transformed.",

                   tex_label = "100(\\beta^{-1} - 1)") # from THANK

    m <= parameter(:ν, 4., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Gamma(2, 0.75), fixed = false,
                   description="ν_l: The inverse Frisch elasticity.",
                   tex_label = "\\nu") # from THANK

    m <= parameter(:ξ_p, 0.85, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.66, 0.1), fixed=false,
                   description="ξ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ξ_p). "
                   * "With probability ξ_p, prices are adjusted according to a weighted average of the previous period's inflation "
                   * "(π_t1) and steady-state inflation (π_ss).",
                   tex_label = "\\xi_p") # from THANK

    m <= parameter(:ξ_w, 0.75, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.66, 0.1), fixed=false,
                   description="ξ_w: (1-ξ_w) is the probability with which households can freely choose wages in each period. "
                   * "With probability ξ_w, wages increase at a geometrically weighted average of the steady state rate of "
                   * "wage increases and last period's productivity times last period's inflation.",
                   tex_label = "\\xi_w") # from THANK

    m <= parameter(:χ, 5., (0., 10.), (1e-5, 0.), ModelConstructors.Exponential(), Gamma(5., 1.), fixed = false,
                   description="χ: The elasticity of the capital utilization cost function.",
                   tex_label = "\\chi") # from THANK

    m <= parameter(:S′′, 2.5, (-15., 15.), (-15., 15.), ModelConstructors.Untransformed(), Gamma(4., 1), fixed = false,
                   description="S′′: The investment adjust cost.",
                   tex_label = "S^{\\prime\\prime}") # from THANK

    m <= parameter(:ϕ_1, 2., (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Normal(1.7, 0.3), fixed = false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_1") # from THANK

    m <= parameter(:ϕ_2, 0.1, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.125, 0.05), fixed = false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label = "\\psi_2") # from THANK

    m <= parameter(:ϕ_3, 0.25, (-0.5, 0.5), (-0.5, 0.5), ModelConstructors.Untransformed(), Normal(0.125, 0.05), fixed = false,

                   description="ψ₃: Weight on ouptut growth in the monetary policy rule.",
                   tex_label = "\\psi_3") # from THANK

    # the following new parameters (from θ to v) are inputs into the steady state calculations given by Giorgio's matlab script
    m <= parameter(:θ, 0.25, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script
    m <= parameter(:τ_x, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script
    m <= parameter(:t_x, 0.2, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script

    m <= parameter(:f_H1, 1/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) #from Giorgio's MATLAB script
    m <= parameter(:f_S1, 3/4, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script
    m <= parameter(:H1, 1.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script
    m <= parameter(:H2, 1.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script
    m <= parameter(:A2, 0.7, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script
    m <= parameter(:s, 0.995, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script
    m <= parameter(:g_x, 0.18, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # from Giorgio's MATLAB script
    m <= parameter(:v, 0.7, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true)
    m <= parameter(:vars, 1.5, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true)
    m <= parameter(:ψ_H1, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true)
    m <= parameter(:ψ_H2, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true)

    # the following parameters (from ζ to sx) are new but had no given value, so current value is random placeholder
    m <= parameter(:ζ, 1.5, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # ???
         # the above parameter appears as \varzeta in the final good producers equilibrium conditions
         # but no \varzeta unicode characters exists in Julia so we use \zeta instead
    m <= parameter(:ζ_X, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # ???
    m <= parameter(:ζ_G, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # ???
    m <= parameter(:ζ_B, 1.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # ???
    m <= parameter(:s_x, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Untransformed(), RootInverseGamma(0.5, 1), fixed = true) # ???

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_R, 0.8, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_R: persistence in the monetary policy rule.",
                   tex_label = "\\rho_{R}") # from THANK

    m <= parameter(:ρ_η_mp, 0.8, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_R: persistence in the monetary policy rule.",
                   tex_label = "\\rho_{R}") # from THANK

    m <= parameter(:ρ_z, 0.25, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_{z}") # from THANK

    m <= parameter(:ρ_a, 0.15, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.4, 0.2), fixed = false,
                   description="ρ_a: AR(1) coefficient", # fix description
                   tex_label = "\\rho_{a}") # ???

    m <= parameter(:ρ_d, 0.15, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.4, 0.2), fixed = false,
                   description="ρ_d: AR(1) coefficient", # fix description
                   tex_label = "\\rho_{d}") # ???

    m <= parameter(:ρ_φ, 0.15, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.4, 0.2), fixed = false,
                   description="ρ_φ: AR(1) coefficient.", # fix description
                   tex_label = "\\rho_{\\varphi}") # ???

    m <= parameter(:ρ_μ, 0.7, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}") # from THANK

    m <= parameter(:ρ_b, 0.7, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_b") # from THANK

    m <= parameter(:ρ_λ_p, 0.95, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_λ_p: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_p}") # from THANK

    m <= parameter(:ρ_λ_w, 0.98, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}") # from THANK

    m <= parameter(:ρ_g, 0.98, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g") # from THANK

    m <= parameter(:ρ_τ_S, 0.98, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_τ_S: AR(1) coefficient.", # fix description
                   tex_label = "\\rho_\\tau_S") # ???

    m <= parameter(:ρ_τ_H1, 0.98, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_τ_H1: AR(1) coefficient.", # fix description
                   tex_label = "\\rho_\\tau_H1") # ???

    m <= parameter(:ρ_τ_H2, 0.98, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_τ_H2: AR(1) coefficient.", # fix description
                   tex_label = "\\rho_g") # ???

    m <= parameter(:ρ_T, 0.8, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_T: parameter governing fiscal rule.", # Update this description later
                   tex_label = "\\rho_{R}") # ???

    m <= parameter(:η_λ_p, 0.0, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.50, 0.20), fixed = false,
                   description="η_λ_p: Moving average component in the price markup shock.",
                   tex_label = "\\eta_{\\lambda_p}") # from THANK

    m <= parameter(:η_λ_w, 0.0, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.50, 0.20), fixed = false,
                   description="η_λ_w: Moving average component in the wage markup shock.",
                   tex_label = "\\eta_{\\lambda_w}") # from THANK

    # exogenous processes - standard deviation
    m <= parameter(:σ_z, 0.9, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_z: The standard deviation of the technology process.",
                   tex_label = "\\sigma_{z}") # from THANK

    m <= parameter(:σ_a, 0.9, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_a: The standard deviation of the COVID-19 shock",
                   tex_label = "\\sigma_{a}") # ???

    m <= parameter(:σ_d, 0.9, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_d: ", # fix description
                   tex_label = "\\sigma_{d}") # ???

    m <= parameter(:σ_φ, 4., (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_φ: The standard deviation of the shock to disutility from labor supply", # fix description
                   tex_label = "\\sigma_{\\mu}") # ???, set to 4, following some values used in literature (e.g. Del Negro, Smets, and Schorfheide (2007) (DSGE-VAR paper))

    m <= parameter(:σ_μ, 5., (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label = "\\sigma_{\\mu}") # from THANK

    m <= parameter(:σ_b, 0.05, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.1, 1), fixed = false,
                   description="σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}") # from THANK

    m <= parameter(:σ_λ_p, 0.15, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.1, 1), fixed = false,
                   description="σ_λ_p: The standard deviation of the price markup shock",
                   tex_label = "\\sigma_{\\lambda_p}") # from THANK

    m <= parameter(:σ_λ_w, 0.2, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.1, 1), fixed = false,
                   description="σ_λ_w: The standard deviation of the wage markup shock",
                   tex_label = "\\sigma_{\\lambda_w}") # from THANK

    m <= parameter(:σ_η_mp, 0.238, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed = false,
                   description="σ_η_mp: The standard deviation of the monetary policy shock.", # fix description
                   tex_label = "\\sigma_{\\eta_mp}") # ???, set to a reasonable calibration from the FRBNY DSGE model

    m <= parameter(:σ_g, 0.35, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}") # from THANK

    m <= parameter(:σ_τ_S, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_τ_S.", # fix description
                   tex_label = "\\sigma_{\\tau_S}") # ???, zeroed out for now

    m <= parameter(:σ_τ_H1, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_τ_H1:.", # fix description
                   tex_label = "\\sigma_{\\tau_{H1}}") # ???, zeroed out for now

    m <= parameter(:σ_τ_H2, 0.0, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_τ_H2: .", # fix description
                   tex_label = "\\sigma_{\\tau_{H2}}") # ???, zeroed out for now

    # steady states
    m <= SteadyStateParameter(:γ, NaN, tex_label = "\\gamma")
    m <= SteadyStateParameter(:β, NaN, tex_label = "\\beta")
    m <= SteadyStateParameter(:r_ss, NaN, tex_label = "r_ss")
    m <= SteadyStateParameter(:π_ss, NaN, tex_label = "\\pi_ss")
    m <= SteadyStateParameter(:π_S_ss, NaN, tex_label = "\\pi_ss")
    m <= SteadyStateParameter(:expγ, NaN, tex_label = "exp^\\gamma")
#=
    m <= SteadyStateParameter(:ψ_H1, NaN, tex_label = "\\psi_H1")
    m <= SteadyStateParameter(:ψ_H2, NaN, tex_label = "\\psi_H2")
=#
    m <= SteadyStateParameter(:L1, NaN, tex_label = "L1")
    m <= SteadyStateParameter(:L2, NaN, tex_label = "L2")
    m <= SteadyStateParameter(:ρ_ss, NaN, tex_label = "\\rho")
    m <= SteadyStateParameter(:w1_ss, NaN, tex_label = "w1")
    m <= SteadyStateParameter(:w2_ss, NaN, tex_label = "w2")
    m <= SteadyStateParameter(:klR1_ss, NaN, tex_label = "klr1_{ss}")
    m <= SteadyStateParameter(:klR2_ss, NaN, tex_label = "klr2_{ss}")
    m <= SteadyStateParameter(:k_ss, NaN, tex_label = "k_{ss}")
    m <= SteadyStateParameter(:F1_ss, NaN, tex_label = "F1_{ss}")
    m <= SteadyStateParameter(:y1_ss, NaN, tex_label = "y1_ss")
    m <= SteadyStateParameter(:y_ss, NaN, tex_label = "y_{ss}")
    m <= SteadyStateParameter(:y2_ss, NaN, tex_label = "y2_{ss}")
#    m <= SteadyStateParameter(:F2_ss, NaN, tex_label = "F2_{ss}")
    m <= SteadyStateParameter(:x_ss, NaN, tex_label = "x_{ss}")
    m <= SteadyStateParameter(:i_S_ss, NaN, tex_label = "i_s_{ss}")
    m <= SteadyStateParameter(:g_ss, NaN, tex_label = "g_{ss}")
    m <= SteadyStateParameter(:c_ss, NaN, tex_label = "c_{ss}")
    m <= SteadyStateParameter(:τ_ss, NaN, tex_label = "\\tau_{ss}")
    m <= SteadyStateParameter(:τ_H1_ss, NaN, tex_label = "\\tau_h1_{ss}")
    m <= SteadyStateParameter(:τ_H2_ss, NaN, tex_label = "\\tau_H2_{ss}")
    m <= SteadyStateParameter(:τ_s_ss, NaN, tex_label = "\\tau_s_{ss}")

    m <= SteadyStateParameter(:t_ss, NaN, tex_label = "t_{ss}")
    m <= SteadyStateParameter(:t_H1_ss, NaN, tex_label = "t_H1_{ss}")
    m <= SteadyStateParameter(:t_H2_ss, NaN, tex_label = "t_H2_{ss}")
    m <= SteadyStateParameter(:t_S_ss, NaN, tex_label = "t_s_{ss}")

    m <= SteadyStateParameter(:c_H1_ss, NaN, tex_label = "c_h10_{ss}")
    m <= SteadyStateParameter(:c_H2_ss, NaN, tex_label = "c_h20_{ss}")
    m <= SteadyStateParameter(:c_S_ss, NaN, tex_label = "c_s0_{ss}")
    m <= SteadyStateParameter(:λ_H1_ss, NaN, tex_label = "\\lambda_h10_{ss}")
    m <= SteadyStateParameter(:λ_H2_ss, NaN, tex_label = "\\lambda_h20_{ss}")
    m <= SteadyStateParameter(:λ_S_ss, NaN, tex_label = "\\lambda_s0_{ss}")
    m <= SteadyStateParameter(:λ_SH1_ss, NaN, tex_label = "\\lambda_sh1_{ss}")
    m <= SteadyStateParameter(:λ_SH2_ss, NaN, tex_label = "\\lambda_sh1_{ss}")
    m <= SteadyStateParameter(:bᴿ_ss, NaN, tex_label = "bR0_{ss}")
    m <= SteadyStateParameter(:R_ss, NaN, tex_label = "R0_{ss}")
    m <= SteadyStateParameter(:mc1_ss, NaN, tex_label = "mc1_{ss}")
    m <= SteadyStateParameter(:mc2_ss, NaN, tex_label = "mc2_{ss}")
    m <= SteadyStateParameter(:price1_ss, NaN, tex_label = "price1_{ss}")
    m <= SteadyStateParameter(:price2_ss, NaN, tex_label = "price2_{ss}")

end

"""
```
steadystate!(m::COTHANK)
```

Calculates the model's steady-state values. `steadystate!(m)` must be called whenever
the parameters of `m` are updated.
"""
function steadystate!(m::COTHANK)

    # These steady state values are inputs given by Giorgio's matlab script
    m[:γ]        = log(1.02)/4
    m[:expγ]     = exp(m[:γ])
    m[:β]        = 0.99
    m[:π_ss]     = 1.02^(1/4)
    m[:π_S_ss]   = 1.02^(1/4)
    m[:r_ss]     = exp(m[:γ]) / m[:β] - 1.

    # Lines 553-595 are taken from Giorgio's MATLAB script for solving steady state, steady states that we
    # are unsure about are in lines 598-604
#    m[:ψ_H1]     = m[:θ] * m[:f_H1] / 1.1 # from Giorgio's MATLAB script for solving steady state
#    m[:ψ_H2]     = m[:θ] * (1 - m[:f_H1]) / 1.1 # from Giorgio's MATLAB script for solving steady state

    m[:L1]   = ((1-m[:θ]) * m[:f_S1] + m[:θ] * m[:f_H1]) * m[:H1] # from Giorgio's MATLAB script for solving steady state
    m[:L2]   = ((1-m[:θ]) * (1 - m[:f_S1]) + m[:θ] * (1 - m[:f_H1])) *m[:H2] # from Giorgio's MATLAB script for solving steady state

    price1ss_price2ss = (m[:v] * m[:L2] * m[:A2] / m[:L1] / (1 - m[:v])) ^ ((1 - m[:α]) / (m[:α] + m[:vars] - m[:α] * m[:vars]))
    m[:price1_ss]     = ( m[:v] + (1 - m[:v]) * price1ss_price2ss ^ (m[:vars] - 1)) ^ (-1 / (1 - m[:vars]))
    m[:price2_ss]     = m[:price1_ss] / price1ss_price2ss

    m[:ρ_ss]    = exp(m[:γ]) / m[:β] - (1 - m[:δ]) # from Giorgio's MATLAB script for solving steady state

    m[:w1_ss]   = (m[:price1_ss] * m[:α] ^ m[:α] * (1 - m[:α]) ^ (1 - m[:α]) / (1 + m[:λ_p_ss]) / (m[:ρ_ss] ^ m[:α])) ^ (1 / (1 - m[:α])) # from Giorgio's MATLAB script for solving steady state
    m[:w2_ss]   = m[:A2] * m[:w1_ss] * price1ss_price2ss ^ (1 / (m[:α] - 1)) # from Giorgio's MATLAB script for solving steady state

    m[:klR1_ss] = m[:w1_ss] / m[:ρ_ss] * m[:α] / (1 - m[:α]) # from Giorgio's MATLAB script for solving steady state
    m[:klR2_ss] = m[:A2] * m[:klR1_ss] * price1ss_price2ss ^ (1 / (m[:α] - 1)) # from Giorgio's MATLAB script for solving steady state
    m[:k_ss]    = m[:klR1_ss] * m[:L1] + m[:klR2_ss] * m[:L2] # from Giorgio's MATLAB script for solving steady state

    m[:F1_ss]   = m[:λ_p_ss] / (1 + m[:λ_p_ss]) * m[:L1] * m[:klR1_ss] ^ m[:α] # from Giorgio's MATLAB script for solving steady state
    m[:y1_ss]   = m[:L1] * m[:klR1_ss] ^ m[:α] - m[:F1_ss] # from Giorgio's MATLAB script for solving steady state
    m[:y_ss]    = m[:y1_ss] / m[:v] * m[:price1_ss] ^ m[:vars] # from Giorgio's MATLAB script for solving steady state
    m[:y2_ss]   = (1 - m[:v]) * m[:y_ss] * m[:price2_ss] ^ (-m[:vars]) # from Giorgio's MATLAB script for solving steady state
#    m[:F2_ss]   = m[:λ_p_ss] * m[:y2_ss] / m[:A2] # unneeded for now
    m[:x_ss]    = m[:y_ss] # from Giorgio's MATLAB script for solving steady state

    m[:i_S_ss]   = (1 - (1 - m[:δ]) * exp(-m[:γ])) * exp(m[:γ]) * m[:k_ss] / (1 - m[:θ]) # from Giorgio's MATLAB script for solving steady state
    m[:g_ss]    = m[:g_x] * m[:x_ss] # from Giorgio's MATLAB script for solving steady state
    m[:c_ss]    = m[:y_ss] - m[:g_ss] - (1 - m[:θ]) * m[:i_S_ss] # from Giorgio's MATLAB script for solving steady state

    m[:τ_ss]    = m[:τ_x] * m[:x_ss] # from Giorgio's MATLAB script for solving steady state
    m[:τ_H1_ss] = m[:τ_ss] * m[:ψ_H1] / m[:θ] / m[:f_H1] # from Giorgio's MATLAB script for solving steady state
    m[:τ_H2_ss] = m[:τ_ss] *m[:ψ_H2] / m[:θ] / (1 - m[:f_H1]) # from Giorgio's MATLAB script for solving steady state
    m[:τ_s_ss]  = (m[:τ_ss] - m[:θ] * (m[:f_H1] * m[:τ_H1_ss] + (1 - m[:f_H1]) * m[:τ_H2_ss])) / (1 - m[:θ])

    m[:t_ss]    = m[:t_x] * m[:x_ss] # from Giorgio's MATLAB script for solving steady state
    m[:t_H1_ss] = m[:t_ss] * m[:ψ_H1] / m[:θ] / m[:f_H1] # from Giorgio's MATLAB script for solving steady state
    m[:t_H2_ss] = m[:t_ss] * m[:ψ_H2] / m[:θ] / (1 - m[:f_H1]) # from Giorgio's MATLAB script for solving steady state
    m[:t_S_ss]  = (m[:t_ss] - m[:θ] * ( m[:f_H1] * m[:t_H1_ss] + (1 - m[:f_H1]) * m[:t_H2_ss])) / (1 - m[:θ]) # from Giorgio's MATLAB script for solving steady state

    # the values below will be used as initial values for solving the fixed point problem
    m[:c_H1_ss] = m[:w1_ss] * m[:H1] - m[:t_H1_ss] + m[:τ_H1_ss] # from Giorgio's MATLAB script for solving steady state
    m[:c_H2_ss] =  m[:w2_ss] * m[:H2] - m[:t_H2_ss] + m[:τ_H2_ss] # from Giorgio's MATLAB script for solving steady state
    m[:c_S_ss]  = (m[:c_ss] - m[:c_H1_ss] * m[:θ] * m[:f_H1] - m[:c_H2_ss] * m[:θ] * (1-m[:f_H1])) / (1 - m[:θ]) # from Giorgio's MATLAB script for solving steady state
    m[:λ_H1_ss] = exp(m[:γ]) / (exp(m[:γ]) - m[:h]) / m[:c_H1_ss] # from Giorgio's MATLAB script for solving steady state
    m[:λ_H2_ss] = exp(m[:γ]) / (exp(m[:γ]) - m[:h]) / m[:c_H2_ss] # from Giorgio's MATLAB script for solving steady state
    m[:λ_S_ss]  = exp(m[:γ]) / (exp(m[:γ]) - m[:h]) / m[:c_S_ss] # from Giorgio's MATLAB script for solving steady state
    m[:R_ss]    = 1.025 / 4.
    m[:bᴿ_ss]   = (m[:g_ss] - m[:t_ss] + m[:τ_ss]) / (1 - m[:R_ss] / exp(m[:γ]) / m[:π_ss]) # from Giorgio's MATLAB script for solving steady state

    # PROBLEM STEADY STATES
    m[:mc1_ss]  = m[:ρ_ss] ^ m[:α] * m[:w1_ss] ^ (1 - m[:α]) / (m[:α] ^ m[:α] * (1 - m[:α]) ^ (1 - m[:α])) / m[:price1_ss]
    m[:mc2_ss]  = m[:ρ_ss] ^ m[:α] * (m[:w2_ss] / m[:A2]) ^ (1 - m[:α]) / (m[:α] ^ m[:α] * (1 - m[:α]) ^ (1 - m[:α])) / m[:price2_ss]
    m[:π_S_ss]  = (m[:y_ss] - m[:price1_ss] * m[:mc1_ss] * m[:L1] * m[:klR1_ss] ^ m[:α] - m[:price2_ss] * m[:mc2_ss] * m[:L2] * m[:klR2_ss] ^ m[:α] * m[:A2] ^
                      (1 - m[:α])) / (1 - m[:θ])

    # function to compute residuals for solving system of nonlinear equations
    function distSS!(r::AbstractVector{S}, x::AbstractVector{S}) where {S <: Real}

        c_H1 = x[1]
        c_H2 = x[2]
        c_S  = x[3]
        λ_H1 = x[4]
        λ_H2 = x[5]
        λ_S  = x[6]
        bR   = x[7]
        R    = x[8]

        # Define residual in r
        r[1] = m[:w1_ss] * m[:H1] - m[:t_H1_ss] + m[:τ_H1_ss] + R * bR * m[:f_S1] * (1 - m[:s]) / m[:θ] / m[:f_H1] / exp(m[:γ]) / m[:π_ss] - c_H1
        r[2] = m[:w2_ss] * m[:H2] - m[:t_H2_ss] + m[:τ_H2_ss] + R * bR * (1 - m[:f_S1]) * (1 - m[:s]) / m[:θ] /(1-m[:f_H1]) / exp(m[:γ]) / m[:π_ss] - c_H2
        r[3] = m[:c_ss] / (1 - m[:θ]) - m[:θ] * (m[:f_H1] * c_H1 + (1 - m[:f_H1]) * c_H2) / (1 - m[:θ]) - c_S

        r[4] = exp(m[:γ]) / (exp(m[:γ]) - m[:h]) / c_H1 - λ_H1
        r[5] = exp(m[:γ]) / (exp(m[:γ]) - m[:h]) / c_H2 - λ_H2
        r[6] = exp(m[:γ]) / (exp(m[:γ]) - m[:h]) / c_S - λ_S

        r[7] = (m[:g_ss] - m[:t_ss] + m[:τ_ss]) / (1 - R / exp(m[:γ]) / m[:π_ss]) - bR

        r[8] = exp(m[:γ]) * m[:π_ss] / m[:β] * λ_S / (m[:s] * λ_S + (1 - m[:s]) * (m[:f_S1] * λ_H1 + (1 - m[:f_S1]) * λ_H2)) - R

        return r
    end

    init = [m[:c_H1_ss].value, m[:c_H2_ss].value, m[:c_S_ss].value, m[:λ_H1_ss].value, m[:λ_H2_ss].value, m[:λ_S_ss].value,
            m[:bᴿ_ss].value, m[:R_ss].value]
    ss = nlsolve(distSS!, init)

    @assert ss.f_converged "Steady state did not converge"

    solution = ss.zero

    m[:c_H1_ss], m[:c_H2_ss], m[:c_S_ss], m[:λ_H1_ss], m[:λ_H2_ss], m[:λ_S_ss], m[:bᴿ_ss], m[:R_ss] = solution

    m[:λ_SH1_ss] = (1 - m[:θ]) * m[:f_S1] / ((1 - m[:θ]) * m[:f_S1] + m[:θ] * m[:f_H1]) * m[:λ_S_ss] + m[:θ] *
                       m[:f_H1] / ((1 - m[:θ]) * m[:f_S1] + m[:θ] * m[:f_H1]) * m[:λ_H1_ss]
    m[:λ_SH2_ss] = (1 - m[:θ]) * (1 - m[:f_S1]) / ((1 - m[:θ]) * (1 - m[:f_S1]) + m[:θ] * (1 - m[:f_H1])) * m[:λ_S_ss] +
                       m[:θ] * (1 - m[:f_H1]) / ((1 - m[:θ]) * (1 - m[:f_S1]) + m[:θ] * (1 - m[:f_H1])) * m[:λ_H2_ss]
    return m
end

function model_settings!(m::COTHANK)

    DSGE.default_settings!(m)

    # Anticipated shocks
    m <= Setting(:n_anticipated_shocks, 0,
                 "Number of anticipated policy shocks")
    m <= Setting(:n_anticipated_shocks_padding, 0,
                 "Padding for anticipated policy shocks")

    # Data
    m <= Setting(:data_id, 3, "Dataset identifier")
    if get_setting(m, :cond_id) in collect(1:5)
        m <= Setting(:cond_full_names, [:obs_gdp, :obs_corepce, :obs_spread, :obs_nominalrate, :obs_longrate],
                     "Observables used in conditional forecasts")
    elseif get_setting(m, :cond_id) == 6
        m <= Setting(:cond_full_names, [:obs_gdp, :obs_corepce, :obs_spread, :obs_nominalrate, :obs_longrate,
                                        :obs_gdpdeflator], "Observables used in conditional forecasts")
    end
    m <= Setting(:cond_semi_names, [:obs_spread, :obs_nominalrate, :obs_longrate],
                 "Observables used in semiconditional forecasts")

    # Forecast
    m <= Setting(:use_population_forecast, true,
                 "Whether to use population forecasts as data")
    m <= Setting(:shockdec_startdate, Nullable(quartertodate("2007-Q1")),
                 "Date of start of shock decomposition output period. If null, then shockdec starts at date_mainsample_start")
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule DuckDB_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("DuckDB")
JLLWrappers.@generate_main_file("DuckDB", UUID("2cbbab25-fc8b-58cf-88d4-687a02676033"))
end  # module DuckDB_jll

ext: jl


text:
include("header.jl")

@testset "sparse" begin

    # Issue #114 (a): plus for Sparse.
    using AutoGrad: Sparse, full, addto!
    a = Sparse(zeros(3,4), Any[ [1.,1.], [1.,1.], 1., 1. ], Any[ ([1,2],), (3:4,), (2,2), (1,) ])
    b = a + a
    @test b isa Sparse
    @test full(b) == full(a) + full(a)
    addto!(a, a)
    @test a isa Sparse
    @test full(a) == full(b)
    b = a + a
    a .+= a
    @test full(a) == full(b)

    # Issue #114 (b): lmul! ambiguous for Sparse, breaks gclip.
    using LinearAlgebra
    foo(w) = (s = 0.0; for i=1:length(w); s+=w[i]; end; s)
    w = Param(randn(2,2))
    J = @diff foo(w)
    @test lmul!(0.5, grad(J,w)) isa Sparse

    # Issue #116: Sparse values and indices should be Vector{Any}
    foo(w) = sum(w[:,[1,2]]) + sum(w[:,[1 2]])
    w = Param(randn(2,2))
    @test @gcheck foo(w)

    # Issue #121: Sparse values and indices should be Vector{Any}
    w,x,y = Param(randn(7)), randn(10), randn(10)
    pred(w,x)=w[5:6]'*tanh.(w[1:2]*x' .+ w[3:4]) .+ w[7]
    @test @gcheck mean(abs2,y-pred(w,x)[1,:])
end


text:
# 1. 生成网络函数
function generateNet(net0,Ntotal,k)
    net = zeros(Int64,Ntotal,k);
    for hang=1:Ntotal
        for lie=1:k
            hl = convert(Int64, net0[hang,lie]);
            net[hang,lie] = hl;
        end
    end
    return net
end

# 2. 更新策略函数
function updateStr(S,strategy,death_num,Ntotal,b,c,k,w)

    # 死亡节点的邻居序列
    death_neigh = S[death_num,:];
    death_neigh_strategy = zeros(Int64,1,4);
    # 死亡节点邻居的策略序列
    death_neigh_strategy = strategy[death_neigh];
    sum_dnei_str = sum(death_neigh_strategy);

    # 更新策略

    # 如果死亡邻居全为D或全为C，则学习该策略
    if sum_dnei_str == 0
        new_strategy = 0;
    elseif sum_dnei_str == k
        new_strategy = 1;
    else
        # 说明死亡节点邻居中既有C也有D，要分别计算邻居的收益
        # 计算死亡节点邻居的收益，就需要找出死亡节点邻居的邻居
        death_neigh_Payoff=zeros(1,k);
        for nn=1:k
            neighNow = death_neigh[nn];
            death_nei_nei = S[neighNow,:];
            # 计算出邻居的邻居的策略序列
            death_nei_nei_strategyC =  sum(strategy[death_nei_nei]);
            # 如果邻居节点i为C策略
            if strategy[neighNow] == 1
                death_neigh_Payoff[nn] = -k * c + b * death_nei_nei_strategyC;
            else
                death_neigh_Payoff[nn] = b * death_nei_nei_strategyC;
            end
        end
        death_neigh_fitness = 1 - w + w * death_neigh_Payoff;
        fitc = sum(death_neigh_fitness .* death_neigh_strategy);
        fitTotal = sum(death_neigh_fitness);
        Pc=fitc / fitTotal;
        probC = rand();
        if probC<=Pc
            new_strategy = 1;
        else
            new_strategy = 0;
        end

    end

    return new_strategy

end

# 3. 主函数
println(Dates.format(now(), "HH:MM:SS"))

# 1.1 变量定义
w = 0.01;
c = 1;
Ncwin = 0;
b = 4.0;

# 1.2 生成网络
net0 = readdlm("N500K4.txt" );
(Ntotal,k) = size(net0);
net = generateNet(net0,Ntotal,k);

# 1.3 轮回循环
for i=1:10000

    # 初始化网络，及其策略分布：1个C，499个D
    S = net;
    strategy = zeros(Int64,Ntotal,1);
    # 随机选择一个个体使其策略为C
    a = rand(1:Ntotal,1);
    strategy[a] = 1;
    # 初始c策略的人数
    nc = 1;

    while ((nc>0) && (nc<Ntotal))

        # 在每一步随机选择一个个体死亡，邻居并根据适应度竞争空位
        death_num = rand(1:Ntotal,1);
        old_str= strategy[death_num];  # 原策略
        old_strategy = old_str[1];
        # 更新策略
        @time new_strategy = updateStr(S,strategy,death_num,Ntotal,b,c,k,w);
        strategy[death_num] = new_strategy;
        nc = nc + new_strategy - old_strategy;

    end

    if nc == Ntotal
        Ncwin = Ncwin + 1;
    end

end

fix_p = Ncwin/10000;
println(fix_p);
println(Dates.format(now(), "HH:MM:SS"))

max_stars_repo_name: dev-juliacn/sailing

max_issues_repo_name: dev-juliacn/sailing

max_forks_repo_name: dev-juliacn/sailing


text:
# ============================================================
# Geometrical Methods
# ============================================================

export global_frame, local_frame, local_velocity, unit_normal, point_distance
export PSpace1D, PSpace2D, CSpace2D, PSpace3D
export uniform_mesh, ndgrid, meshgrid, find_idx
export UnstructPSpace

include("geo_general.jl")
include("geo_struct.jl")
include("geo_unstruct.jl")


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

using Test

include("MOI/solvers.jl")
include("MOI/wrapper.jl")
include("MOI/minlptests.jl")

# GAMS is really fragile on macos and windows
if Sys.islinux()
  GC.gc()

  include("MOI/solvers_fragile.jl")
  include("MOI/wrapper.jl")
end


text:
#options.jl  - sets up the options dictionary and methods to easily set options
#from within command line or using a unums.jl preprocessing block.
doc"""
  *Unums.options* is the associative dictionary of options for this module.
  options can be acessed or set using the associated symbols.

  * :devmode => true          sets developer mode protections on.
  * :usegnum => false         sets it so that operands return gnums
  * :env_ESS => 4             sets the default environment esizesize
  * :env_FSS => 6             sets the default environment fsizesize
  * :longform => false        outputs unum types as their long form

  NB: usegnum will eventually be set to a default "true"

  *overriding using ~/.unums.jl*

  You can specify a julia script in ~/.unums.jl which will be executed before running any unums script.  While
  this can be used for anything, you can override options using the following code:

    RexAssembly.options[*:option*] = *value*

  *overriding using command line arguments*

  RexAssembly automatically parses the command line arguments to see if the options should be set, this
  takes precedence over internal defaults and contents of ~/.unums.jl.  For example:

    >script.jl devmode=false

    >julia unumscript.jl devmode=false

  both would disable developer mode checks.
"""
const options = Dict{Symbol, Any}(
  :devmode => true,
  :usegnum => false,
  :env_ESS => 4,
  :env_FSS => 6,
  :longform => false)

#unixey operating systems should look for a .unums.jl file in the home directary
#and use that as a way to set preferences.
if (is_linux() || is_apple() || is_bsd())
  isfile("$(homedir)/.unums.jl") && include("$(homedir)/.unums.jl")
  isfile("$(homedir)/.unums.jl") && include("$(homedir)/.unums.jl")
end

function argparse(s)
  (s == "true") && return true
  (s == "false") && return false
  isnumber(s) && return parse(s)
end

for arg in ARGS
  arglist = split(arg,"=")
  if length(arglist) == 1
    options[Symbol(arglist[1])] = true
  else
    options[Symbol(arglist[1])] = argparse(arglist[2])
  end
end

################################################################################
## include code for options.

include("./environment.jl")
include("./devsafety.jl")
include("./usegnum.jl")
include("./longform.jl")

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using YAML
using Test
using Streamfall

TEST_DIR = @__DIR__

@testset "Bare node creation" begin
    # Test direct running of a single time step
    ihacres = create_node(BilinearNode, "IHACRES", 100.0)
    @test Streamfall.run_timestep!(ihacres, 6.0, 3.0, 1) isa Any

    # Expuh form does not yet support time stepping
    expuh = create_node(ExpuhNode, "Expuh", 100.0)
    @test Streamfall.run_node!(expuh, 6.0, 3.0, 50.0, 10.0, 5.0) isa Any

    gr4j = create_node(GR4JNode, "GR4J", 100.0)
    @test Streamfall.run_timestep!(gr4j, 6.0, 3.0, 1) isa Any

    hymod = create_node(SimpleHyModNode, "HyMod", 100.0)
    @test Streamfall.run_timestep!(hymod, 6.0, 3.0, 1) isa Any

    symhyd = create_node(SYMHYDNode, "SYMHYD", 100.0)
    @test Streamfall.run_timestep!(symhyd, 6.0, 3.0, 1) isa Any
end

@testset "Ensure no NaN outputs" begin
    test_node = BilinearNode(
        "Test",  # name/id
        1985.73,  # area
        200.0,  # d
        2.0,  # d2
        1.0,  # e
        1.675,  # f
        54.35254,  # a
        0.187,  # b
        2.9,  # gw storage_coef
        0.727,  # alpha
        100.0,  # storage
        0.0,  # initial quickflow
        0.0,  # initial slowflow
        0.0  # initial gw store
    )

    rain = 7.96848605e+01
    evap = 3.32467909e+00
    inflow = 9.19583373e-11
    extraction = 1.59313987e-12
    gw_exchange = 2.47076926e-11
    current_store = 1.03467364e+02
    quick_store = 8.46269687e+02
    slow_store = 3.67133471e+02

    res = run_step!(test_node, rain, evap, inflow, extraction, gw_exchange;
                    current_store=current_store, quick_store=quick_store, slow_store=slow_store)

    @test !any(isnan, res)
end

@testset "Network creation" begin
    # Ensure specified parameter values are being assigned on node creation
    # Load and generate stream network
    network = YAML.load_file(joinpath(TEST_DIR, "data/campaspe/campaspe_network.yml"))
    sn = create_network("Example Network", network)

    target_node = get_prop(sn, 1, :node)

    @test target_node.area == 1985.73
    @test target_node.level_params[1] == -3.3502
end

@testset "Interim CMD (no NaN poisoning)" begin
    params = (214.6561105573191, 76.6251447, 200.0, 2.0, 0.727)
    current_store, rain, d, d2, alpha = params

    interim_results = [0.0, 0.0, 0.0]
    @ccall IHACRES.calc_ft_interim_cmd(interim_results::Ptr{Cdouble},
                                       current_store::Cdouble,
                                       rain::Cdouble,
                                       d::Cdouble,
                                       d2::Cdouble,
                                       alpha::Cdouble)::Cvoid

    (mf, e_rainfall, recharge) = interim_results

    @test !isnan(mf)
    @test !isnan(e_rainfall)
    @test !isnan(recharge)
end

@testset "Catchment Moisture Deficit" begin
    cmd = 100.0
    et = 6.22
    e_rain = 6.83380027058404E-06
    recharge = 3.84930005080411E-06
    rain = 0.0000188

    n_cmd = @ccall IHACRES.calc_cmd(cmd::Cdouble, rain::Cdouble, et::Cdouble, e_rain::Cdouble, recharge::Cdouble)::Float64

    @test isapprox(n_cmd, 106.22, atol=0.001)
end

@testset "IHACRES calculations" begin
    area = 1985.73
    a = 54.352
    b = 0.187
    e_rain = 3.421537294474909e-6
    recharge = 3.2121031313153022e-6

    prev_quick = 100.0
    prev_slow = 100.0

    flow_results = [0.0, 0.0, 0.0]
    @ccall IHACRES.calc_ft_flows(
        flow_results::Ptr{Cdouble},
        prev_quick::Cdouble,
        prev_slow::Cdouble,
        e_rain::Cdouble,
        recharge::Cdouble,
        area::Cdouble,
        a::Cdouble,
        b::Cdouble
    )::Cvoid

    quickflow = (prev_quick + (e_rain * area))
    α = exp(-a)
    quickflow = α * quickflow

    @test flow_results[1] == quickflow

    slow_store = prev_slow + (recharge * area)
    α = exp(-b)
    slow_store = α * slow_store
    @test flow_results[2] == slow_store

    e_rain = 0.0
    recharge = 0.0

    prev_quick = 3.3317177943791187
    prev_slow = 144.32012122323678

    flow_results = [0.0, 0.0, 0.0]
    @ccall IHACRES.calc_ft_flows(
        flow_results::Ptr{Cdouble},
        prev_quick::Cdouble,
        prev_slow::Cdouble,
        e_rain::Cdouble,
        recharge::Cdouble,
        area::Cdouble,
        a::Cdouble,
        b::Cdouble
    )::Cvoid

    quickflow = (prev_quick + (e_rain * area))
    α = exp(-a)
    @test flow_results[1] == (α * quickflow)

    slow_store = prev_slow + (recharge * area)
    α = exp(-b)
    β = (1.0 - α) * slow_store
    slow_store = α * slow_store
    @test flow_results[2] == slow_store
end

include("test_metrics.jl")
include("test_data_op.jl")
include("test_networks.jl")
include("test_calibration.jl")


text:
module K8sClusterManagers

using DataStructures: DefaultOrderedDict, OrderedDict
using Distributed: Distributed, ClusterManager, WorkerConfig, cluster_cookie
using JSON: JSON
using Mocking: Mocking, @mock
using kubectl_jll

export K8sClusterManager, KubeError, isk8s

function __init__()
    if !kubectl_jll.is_available()
        error("kubectl_jll does not support the current platform. See: ",
              "https://github.com/JuliaBinaryWrappers/kubectl_jll.jl#platforms")
    end
end

include("namespace.jl")
include("pod.jl")
include("native_driver.jl")

end


text:
module Types

abstract type AbstractDelayed end

const SPACE_TYPE = Union{
    Dict{Symbol, T} where T,
    AbstractDelayed,
    AbstractVector,
}

end # module


text:
#module TestResampling

using Distributed
import ComputationalResources: CPU1, CPUProcesses, CPUThreads
using .TestUtilities
using ProgressMeter
import Tables

@everywhere begin
    using .Models
    using StableRNGs
    rng = StableRNG(1513515)
    const verb = 0
end

using Test
using MLJBase
import Distributions
import StatsBase
@static if VERSION >= v"1.3.0-DEV.573"
    using .Threads
end

struct DummyInterval <: Interval end
dummy_interval=DummyInterval()

dummy_measure_det(yhat, y) = 42
MLJBase.target_scitype(::typeof(dummy_measure_det)) = Table(MLJBase.Textual)
MLJBase.prediction_type(::typeof(dummy_measure_det)) = :deterministic

dummy_measure_interval(yhat, y) = [123, 456]
MLJBase.target_scitype(::typeof(dummy_measure_interval)) =
    Table(MLJBase.Textual)
MLJBase.prediction_type(::typeof(dummy_measure_interval)) = :interval

@testset "_actual_operations" begin
    clf = ConstantClassifier()
    rgs = ConstantRegressor()
    clf_det = DeterministicConstantClassifier()
    rgs_det = DeterministicConstantRegressor()
    measures = [LogLoss(), Accuracy(), BrierScore()] # mixed prob/determ
    measures_det = [Accuracy(), FScore()]
    operations = [predict, predict_mode, predict]

    # single measure gets replicated to match length of `measures`:
    @test MLJBase._actual_operations(predict_mean,
                                   [Accuracy(), FScore()],
                                    clf,
                                    1) ==
                                   [predict_mean, predict_mean]

    # handling of a measure with `:unknown` `prediction_type` (eg,
    # custom measure):
    my_mae(yhat, y) = abs.(yhat - y)
    @test(
        MLJBase._actual_operations(nothing, [my_mae, LPLoss()], rgs_det , 1) ==
        [predict, predict])

    @test MLJBase._actual_operations(predict, [LogLoss(),], clf, 1) ==
        [predict,]
    @test MLJBase._actual_operations(operations, measures, clf, 1) == operations
    @test_throws MLJBase.ERR_OPERATION_MEASURE_MISMATCH _ =
        MLJBase._actual_operations([predict, predict_mode], measures, clf, 1)
    @test_throws MLJBase.ERR_OPERATION_MEASURE_MISMATCH junk =
        MLJBase._actual_operations([predict,], measures, clf, 1)
    @test_throws MLJBase.ERR_INVALID_OPERATION _ =
        MLJBase._actual_operations(transform, [LogLoss(),], clf, 1)
    @test MLJBase._actual_operations(nothing, measures, clf, 1) == operations
    @test(
       @test_logs MLJBase._actual_operations(nothing, [Accuracy(),], clf, 1) ==
           [predict_mode])
    @test MLJBase._actual_operations(nothing, [l2,], rgs, 1) ==
        [predict_mean, ]
    @test_throws(MLJBase.err_incompatible_prediction_types(clf_det, LogLoss()),
                 MLJBase._actual_operations(nothing, [LogLoss(),], clf_det, 1))
    @test MLJBase._actual_operations(nothing, measures_det, clf_det, 1) ==
        [predict, predict]

    # measure/model differ in prediction type but weird target_scitype:
    @test_throws(
        MLJBase.err_ambiguous_operation(clf, dummy_measure_det),
        MLJBase._actual_operations(nothing, [dummy_measure_det, ], clf, 1))

    # measure has :interval prediction type but model does not (2 cases):
    @test_throws(
        MLJBase.err_ambiguous_operation(clf, dummy_measure_interval),
        MLJBase._actual_operations(nothing,
                                   [dummy_measure_interval, ], clf, 1))
    @test_throws(
        MLJBase.err_ambiguous_operation(clf_det, dummy_measure_interval),
        MLJBase._actual_operations(nothing,
                                   [dummy_measure_interval, ], clf_det, 1))

    # both measure and model have :interval prediction type:
    @test MLJBase._actual_operations(nothing,
                               [dummy_measure_interval, ],
                               dummy_interval, 1) == [predict, ]

    # model has :interval prediction type but measure does not:
    @test_throws(
        MLJBase.err_ambiguous_operation(dummy_interval, LogLoss()),
        MLJBase._actual_operations(nothing,
                                   [LogLoss(), ], dummy_interval, 1))
end

@testset "_feature_dependencies_exist" begin
    measures = Any[rms, rsq, log_loss, brier_score]
    @test !MLJBase._feature_dependencies_exist(measures)
    my_feature_dependent_loss(ŷ, X, y) =
        sum(abs.(ŷ - y) .* X.penalty)/sum(X.penalty);
    MLJBase.is_feature_dependent(::typeof(my_feature_dependent_loss)) = true
    push!(measures, my_feature_dependent_loss)
    @test MLJBase._feature_dependencies_exist(measures)
end

@testset_accelerated "dispatch of resources and progress meter" accel begin

    @info "Checking progress bars:"

    X = (x = [1, ],)
    y = [2.0, ]

    @everywhere begin
        nfolds = 6
        nmeasures = 2
        func(mach, k) = (
            (sleep(MLJBase.PROG_METER_DT*rand(rng)); fill(1:k, nmeasures)),
            :fitted_params,
            :report,
        )
    end
    mach = machine(ConstantRegressor(), X, y)
    if accel isa CPUThreads
        result = MLJBase._evaluate!(
            func,
            mach,
            CPUThreads(Threads.nthreads()),
            nfolds,
            1
        )
    else
        result = MLJBase._evaluate!(func, mach, accel, nfolds, 1)
    end
    measurements = vcat(result[1]...)
    @test measurements ==
        [1:1, 1:1, 1:2, 1:2, 1:3, 1:3, 1:4, 1:4, 1:5, 1:5, 1:6, 1:6]
    @test collect(result[2]) == fill(:fitted_params, nfolds)
end

@test CV(nfolds=6) == CV(nfolds=6)
@test CV(nfolds=5) != CV(nfolds=6)
@test MLJBase.train_test_pairs(CV(), 1:10) !=
     MLJBase.train_test_pairs(CV(shuffle=true), 1:10)
@test MLJBase.train_test_pairs(Holdout(), 1:10) !=
     MLJBase.train_test_pairs(Holdout(shuffle=true), 1:10)

@testset "train test pairs" begin
    cv = CV(nfolds=5)

    pairs = MLJBase.train_test_pairs(cv, 1:24)
    @test pairs == [
        (6:24, 1:5),
        ([1:5..., 11:24...], 6:10),
        ([1:10..., 16:24...], 11:15),
        ([1:15..., 21:24...], 16:20),
        (1:20, 21:24)
    ]

    # Not enough data for the number of folds.
    @test_throws ArgumentError MLJBase.train_test_pairs(cv, 1:4)
end

@testset "checking measure/model compatibility" begin
    model = ConstantRegressor()
    y = rand(rng,4)

    # model prediction type is Probablistic but measure is Deterministic:
    @test_throws(ArgumentError,
                  MLJBase._check_measure(rms, predict, model, y))

    @test MLJBase._check_measure(rms, predict_mean, model, y)

    @test MLJBase._check_measure(rms, predict_median, model, y)

    # has `y`  `Finite` elscityp but measure `rms` is for `Continuous`:
    y=categorical(collect("abc"))
    @test_throws(ArgumentError,
                 MLJBase._check_measure(rms, predict_median, model, y))
    model = ConstantClassifier()
    # model prediction type is Probablistic but measure is Deterministic:
    @test_throws(ArgumentError,
                 MLJBase._check_measure(mcr, predict, model, y))

    @test MLJBase._check_measure(mcr, predict_mode, model, y)

    # `Determistic` model but `Probablistic` measure:
    model = DeterministicConstantClassifier()
    @test_throws(ArgumentError,
                 MLJBase._check_measure(cross_entropy, predict, model, y))

    # measure with wrong target_scitype:
    @test_throws(ArgumentError,
                 MLJBase._check_measures([brier_score, rms],
                                         [predict_mode, predict_mean],
                                         model, y))

    model = ConstantClassifier()
    @test MLJBase._check_measures([brier_score, cross_entropy, accuracy],
                                  [predict, predict, predict_mode],
                                  model, coerce(y, Multiclass))
end

@testset "check weights" begin
    @test_throws(MLJBase.ERR_WEIGHTS_REAL,
                 MLJBase._check_weights([:junk, :junk], 2))
    @test_throws(MLJBase.ERR_WEIGHTS_LENGTH,
                 MLJBase._check_weights([0.5, 0.5], 3))
    @test MLJBase._check_weights([0.5, 0.5], 2)
end

@testset "check class weights" begin
    w = Dict('a'=> 0.2, 'b'=>0.8)
    @test_throws(MLJBase.ERR_WEIGHTS_DICT,
                 MLJBase._check_class_weights([0.1, 0.4], ['a', 'b']))
    @test_throws(MLJBase.ERR_WEIGHTS_CLASSES,
                 MLJBase._check_class_weights(w, ['a', 'c']))
    @test MLJBase._check_class_weights(w, ['b', 'a'])
end

@testset_accelerated "folds specified" accel begin
    x1 = ones(10)
    x2 = ones(10)
    X  = (x1=x1, x2=x2)
    y  = [1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0]

    my_rms(yhat, y) = sqrt(mean((yhat -y).^2))
    my_mae(yhat, y) = abs.(yhat - y)
    MLJBase.reports_each_observation(::typeof(my_mae)) = true
    MLJBase.prediction_type(::typeof(my_rms)) = :deterministic
    MLJBase.prediction_type(::typeof(my_mae)) = :deterministic

    resampling = [(3:10, 1:2),
                  ([1, 2, 5, 6, 7, 8, 9, 10], 3:4),
                  ([1, 2, 3, 4, 7, 8, 9, 10], 5:6),
                  ([1, 2, 3, 4, 5, 6, 9, 10], 7:8),
                  (1:8, 9:10)]

    for cache in [true, false]

        model = DeterministicConstantRegressor()
        mach  = machine(model, X, y, cache=cache)

        # check detection of incompatible measure (cross_entropy):
        @test_throws ArgumentError evaluate!(mach, resampling=resampling,
                                             measure=[cross_entropy, rmslp1],
                                             verbosity=verb,
                                             acceleration=accel)
        result = evaluate!(mach, resampling=resampling, verbosity=verb,
                           measure=[my_rms, my_mae, rmslp1], acceleration=accel)

        v = [1/2, 3/4, 1/2, 3/4, 1/2]

        @test result.per_fold[1] ≈ v
        @test result.per_fold[2] ≈ v
        @test result.per_fold[3][1] ≈ abs(log(2) - log(2.5))
        @test ismissing(result.per_observation[1])
        @test result.per_observation[2][1] ≈ [1/2, 1/2]
        @test result.per_observation[2][2] ≈ [3/4, 3/4]
        @test result.measurement[1] ≈ mean(v)
        @test result.measurement[2] ≈ mean(v)

        # fitted_params and report per fold:
        @test map(fp->fp.fitresult, result.fitted_params_per_fold) ≈
            [1.5, 1.25, 1.5, 1.25, 1.5]
        @test all(==(NamedTuple()), result.report_per_fold)
    end
end

@testset "repeated resampling" begin
    x1 = ones(20)
    x2 = ones(20)
    X = (x1=x1, x2=x2)
    y = rand(rng, 20)

    holdout = Holdout(fraction_train=0.75, rng=rng)
    model = Models.DeterministicConstantRegressor()

    for cache in [true, false]
        mach = machine(model, X, y, cache=cache)
        result = evaluate!(mach, resampling=holdout, verbosity=verb,
                       measure=[rms, rmslp1], repeats=6)
        per_fold = result.per_fold[1]
        @test unique(per_fold) |> length == 6
        @test abs(mean(per_fold) - std(y)) < 0.06 # very rough check

        cv = CV(nfolds=3, rng=rng)
        result = evaluate!(mach, resampling=cv, verbosity=verb,
                       measure=[rms, rmslp1], repeats=6)
        per_fold = result.per_fold[1]
        @test unique(per_fold) |> length == 18
        @test abs(mean(per_fold) - std(y)) < 0.06 # very rough check
    end
end

@testset_accelerated "holdout" accel begin
    x1 = ones(4)
    x2 = ones(4)
    X = (x1=x1, x2=x2)
    y = [1.0, 1.0, 2.0, 2.0]

    @test MLJBase.show_as_constructed(Holdout)
    holdout = Holdout(fraction_train=0.75)
    model = Models.DeterministicConstantRegressor()
    for cache in [true, false]
        mach = machine(model, X, y, cache=cache)
        result = evaluate!(mach, resampling=holdout, verbosity=verb,
                           measure=[rms, rmslp1], acceleration=accel)
        result = evaluate!(mach, resampling=holdout, verbosity=verb,
                           acceleration=accel)
        @test result.measurement[1] ≈ 2/3

        # test direct evaluation of a model + data:
        result = evaluate(model, X, y, verbosity=0,
                          resampling=holdout, measure=rms, cache=cache)
        @test result.measurement[1] ≈ 2/3
    end

    X = (x=rand(rng,100),)
    y = rand(rng,100)

    for cache in [true, false]
        mach = machine(model, X, y, cache=cache)
        evaluate!(mach, verbosity=verb,
                  resampling=Holdout(shuffle=true, rng=rng), acceleration=accel)
        e1 = evaluate!(mach, verbosity=verb,
                       resampling=Holdout(shuffle=true),
                       acceleration=accel).measurement[1]
        @test e1 != evaluate!(mach, verbosity=verb,
                              resampling=Holdout(),
                              acceleration=accel).measurement[1]
    end
end

@testset_accelerated "Exception handling (see issue 235)" accel begin
    X, y = make_moons(50)
    model = ConstantClassifier()

    bad_loss(yhat, y) = throw(Exception())
    @test_throws Exception evaluate(model, X, y, measure=bad_loss, verbosity=0)
end

@testset_accelerated "cv" accel begin
    x1 = ones(10)
    x2 = ones(10)
    X = (x1=x1, x2=x2)
    y = [1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0]

    @test MLJBase.show_as_constructed(CV)
    cv=CV(nfolds=5)
    for cache in [true, false]
        model = Models.DeterministicConstantRegressor()
        mach = machine(model, X, y, cache=cache)
        result = evaluate!(mach, resampling=cv, measure=[rms, rsq, rmslp1],
                           acceleration=accel, verbosity=verb)

        @test result.per_fold[1] ≈ [1/2, 3/4, 1/2, 3/4, 1/2]

        shuffled = evaluate!(mach, resampling=CV(shuffle=true), verbosity=verb,
                             acceleration=accel) # using rms default
        @test shuffled.measurement[1] != result.measurement[1]
    end
end

@testset "TimeSeriesCV" begin
    tscv = TimeSeriesCV(; nfolds=3)

    pairs = MLJBase.train_test_pairs(tscv, 1:10)
    @test pairs == [
        (1:4, [5, 6]),
        (1:6, [7, 8]),
        (1:8, [9, 10])
    ]

    pairs = MLJBase.train_test_pairs(tscv, 1:2:15)
    @test pairs == [
        ([1, 3], [5, 7])
        ([1, 3, 5, 7], [9, 11])
        ([1, 3, 5, 7, 9, 11], [13, 15])
    ]

    @test_logs(
        (:warn, "TimeSeriesCV is being applied to `rows` not in sequence. "),
        MLJBase.train_test_pairs(tscv, reverse(1:10))
    )

    # Not enough data for the number of folds.
    @test_throws ArgumentError MLJBase.train_test_pairs(TimeSeriesCV(10), 1:8)
end

@testset "stratified_cv" begin

    # check in explicit example:
    y = categorical(
        ['b', 'c', 'a', 'b', 'b', 'b', 'c', 'c', 'c', 'a', 'a', 'a']
    )
    scv = StratifiedCV(nfolds=3)
    rows = 1:12
    pairs = MLJBase.train_test_pairs(scv, rows, y)
    expected_pairs = [
        ([2, 4, 9, 11, 3, 5, 7, 12], [1, 6, 8, 10]),
        ([1, 6, 8, 10, 3, 5, 7, 12], [2, 4, 9, 11]),
        ([1, 6, 8, 10, 2, 4, 9, 11], [3, 5, 7, 12])
    ]
    # Explanation of expected_pairs: The row indices are processed one at
    # a time. The test fold that a row index is placed in is determined
    # by this lookup:
    #
    # b b b b c c c c a a a a
    # 1 2 3 1 2 3 1 2 3 1 2 3
    #
    # For example, the first row such that y[row] == 'c' is placed in the
    # second fold, and the second row such that y[row] == 'c' is placed in
    # the third fold.
    @test pairs == expected_pairs

    # test invariance to label renaming:
    z = replace(y, 'a' => 'b', 'b' => 'c', 'c' => 'a')
    pairs = MLJBase.train_test_pairs(scv, rows, z)
    @test pairs == expected_pairs

    # test the case where rows is a shuffled subset of y:
    y = categorical(['a', 'b', 'c', 'a', 'b', 'c', 'a', 'b', 'c'])
    rows = 8:-1:2
    pairs = MLJBase.train_test_pairs(scv, rows, y)
    @test pairs == [
        ([5, 4, 6, 2], [8, 7, 3]),
        ([8, 7, 3, 6, 2], [5, 4]),
        ([8, 7, 3, 5, 4], [6, 2])
    ]

    # test shuffle:
    scv_random = StratifiedCV(nfolds=3, shuffle=true, rng=1)
    pairs_random = MLJBase.train_test_pairs(scv_random, rows, y)
    @test pairs != pairs_random

    # wrong target type throws error:
    @test_throws(ArgumentError,
                 MLJBase.train_test_pairs(scv,
                                          rows,
                                          CategoricalArrays.unwrap.(y)))

    # check class distribution is preserved in a larger randomized example:
    N = 30
    y = shuffle(vcat(fill('a', N), fill('b', 2N),
                        fill('c', 3N), fill('d', 4N))) |> categorical;
    d = Distributions.fit(MLJBase.UnivariateFinite, y)
    pairs = MLJBase.train_test_pairs(scv, 1:10N, nothing, y)
    folds = vcat(first.(pairs), last.(pairs))
    @test all([Distributions.fit(MLJBase.UnivariateFinite, y[fold]) ≈
               d for fold in folds])
end

@testset_accelerated "sample weights in evaluation" accel begin
    # cv:
    x1 = ones(4)
    x2 = ones(4)
    X = (x1=x1, x2=x2)
    y = [1.0, 2.0, 3.0, 1.0]
    w = 1:4
    cv=CV(nfolds=2)
    model = Models.DeterministicConstantRegressor()
    mach = machine(model, X, y)
    e = evaluate!(mach, resampling=cv, measure=l1,
                  weights=w, verbosity=verb, acceleration=accel).measurement[1]
    efold1 = mean([1*1, 1*0])
    efold2 = mean([3*3/2, 4*1/2])
    @test e ≈ mean([efold1, efold2])

    # if I don't specify weights in `evaluate!`, then uniform should
    # be used:
    e = evaluate!(mach, resampling=cv, measure=l1,
                  verbosity=verb, acceleration=accel).measurement[1]
    efold1 = mean([1*1, 1*0])
    efold2 = mean([1*3/2, 1*1/2])
    @test e ≈ mean([efold1, efold2])
end

@testset_accelerated "class weights in evaluation" accel begin
    x = [1,2,3,4,5,6,7]
    X, y = Tables.table([x x x x x x]), coerce([1,2,1,3,1,2,2], Multiclass)
    model = Models.DeterministicConstantClassifier()
    mach = machine(model, X, y)
    cv=CV(nfolds = 2)
    class_w = Dict(1=>1, 2=>2, 3=>3)
    e = evaluate!(mach,
                  resampling=cv,
                  measure=MLJBase.MulticlassFScore(return_type=Vector),
                  class_weights=class_w,
                  verbosity=verb,
                  acceleration=accel).measurement[1]
    @test round(e, digits=3) ≈ 0.217

    # if class weights in `evaluate!` isn't specified:
    e = evaluate!(mach, resampling=cv, measure=multiclass_f1score,
                  verbosity=verb, acceleration=accel).measurement[1]
    @test e ≈ 0.15
end

@testset_accelerated "resampler as machine" accel begin
    N = 50
    X = (x1=rand(rng, N), x2=rand(rng, N), x3=rand(rng, N))
    y = X.x1 -2X.x2 + 0.05*rand(rng, N)
    ridge_model = FooBarRegressor(lambda=20.0)
    holdout = Holdout(fraction_train=0.75)
    resampler = Resampler(resampling=holdout, model=ridge_model, measure=mae,
                          acceleration=accel)
    resampling_machine = machine(resampler, X, y)
    @test_logs((:info, r"^Training"), fit!(resampling_machine))
    e1=evaluate(resampling_machine).measurement[1]
    mach = machine(ridge_model, X, y)
    @test e1 ≈  evaluate!(mach, resampling=holdout,
                          measure=mae, verbosity=verb,
                          acceleration=CPU1()).measurement[1]
    ridge_model.lambda=1.0
    fit!(resampling_machine, verbosity=verb)
    e2=evaluate(resampling_machine).measurement[1]
    @test e1 != e2
    resampler.weights = rand(rng,N)
    fit!(resampling_machine, verbosity=verb)
    e3=evaluate(resampling_machine).measurement[1]
    @test e3 != e2

    @test MLJBase.package_name(Resampler) == "MLJBase"
    @test MLJBase.is_wrapper(Resampler)

    # when only `model` changes, the folds shouldn't change, even in
    # shuffled case:
    cv = CV(rng=StableRNGs.StableRNG(123))
    resampler=Resampler(model=ridge_model,
                        resampling=cv,
                        repeats=3,
                        measure=mae, acceleration=accel)
    mach = machine(resampler, X, y)
    fit!(mach, verbosity=verb)
    ev1 = evaluate(mach)
    rows1 = ev1.train_test_rows
    resampler.model.lambda *= 0.5
    fit!(mach, verbosity=verb)
    ev2 = evaluate(mach)
    @test rows1 == ev2.train_test_rows
    resampler.model.lambda *= 2
    fit!(mach, verbosity=verb)
    @test ev1.measurement[1] ≈ evaluate(mach).measurement[1]

    # but if `resampling` or `repeats`, then new
    # folds should be generated:
    resampler.resampling = CV(rng=cv.rng, nfolds=2)
    fit!(mach, verbosity=verb)
    rows2 = evaluate(mach).train_test_rows
    @test length(rows2) == 2 * resampler.repeats
    @test rows2 != rows1
    resampler.repeats += 1
    fit!(mach, verbosity=verb)
    rows3 = evaluate(mach).train_test_rows
    @test length(rows3) == 2 * resampler.repeats
end

struct DummyResamplingStrategy <: MLJBase.ResamplingStrategy end

@testset_accelerated "custom strategy depending on X, y" accel begin
    function MLJBase.train_test_pairs(resampling::DummyResamplingStrategy,
                              rows, X, y)
        train = filter(rows) do j
            y[j] == y[1]
        end
        test = setdiff(rows, train)
        return [(train, test),]
    end

    X = (x = rand(rng,8), )
    y = categorical(string.([:x, :y, :x, :x, :y, :x, :x, :y]))
    @test MLJBase.train_test_pairs(DummyResamplingStrategy(), 2:6, X, y) ==
        [([3, 4, 6], [2, 5]),]

    e = evaluate(ConstantClassifier(), X, y,
                 measure=misclassification_rate,
                 resampling=DummyResamplingStrategy(),
                 operation=predict_mode,
                 acceleration=accel,
                 verbosity=verb)
    @test e.measurement[1] ≈ 1.0
end

@testset_accelerated "sample weights in training and evaluation" accel begin
    yraw = ["Perry", "Antonia", "Perry", "Antonia", "Skater"]
    X = (x=rand(rng,5),)
    y = categorical(yraw)
    w = [1, 10, 1, 10, 5]

    for cache in [true, false]
        # without weights:
        mach = machine(ConstantClassifier(), X, y, cache=cache)
        e = evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                      operation=predict_mode, measure=misclassification_rate,
                      acceleration=accel, verbosity=verb)
        @test e.measurement[1] ≈ 1.0

        # with weights in training and evaluation:
        mach = machine(ConstantClassifier(), X, y, w, cache=cache)
        e = evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                      operation=predict_mode, measure=misclassification_rate,
                      acceleration=accel, verbosity=verb, weights=w)
        @test e.measurement[1] ≈ mean([10*0, 5*1])

        # with different weights in training and evaluation:
        e = evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                      operation=predict_mode, measure=misclassification_rate,
                      weights = fill(1, 5), acceleration=accel, verbosity=verb)
        @test e.measurement[1] ≈ 1/2

        @test_throws(DimensionMismatch,
                     evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                               operation=predict_mode,
                               measure=misclassification_rate,
                               weights = fill(1, 100), acceleration=accel,
                               verbosity=verb))

        @test_throws(ArgumentError,
                     evaluate!(mach, resampling=Holdout(fraction_train=0.6),
                               operation=predict_mode,
                               measure=misclassification_rate,
                               weights = fill('a', 5), acceleration=accel,
                               verbosity=verb))
    end

    # resampling on a subset of all rows:
    model = KNNClassifier()

    N = 200
    X = (x = rand(rng,3N), );
    y = categorical(rand(rng,"abcd", 3N));
    w = rand(rng,3N);
    class_w = Dict(zip(levels(y), rand(length(levels(y)))));
    rows = StatsBase.sample(1:3N, 2N, replace=false);
    Xsmall = selectrows(X, rows);
    ysmall = selectrows(y, rows);
    wsmall = selectrows(w, rows);

    for cache in [true, false]
        mach1 = machine(model, Xsmall, ysmall, wsmall, cache=cache)
        e1 = evaluate!(mach1,
                       resampling=CV(),
                       measure=misclassification_rate,
                       weights=wsmall,
                       operation=predict_mode,
                       acceleration=accel,
                       verbosity=verb)

        mach2 = machine(model, X, y, w, cache=cache)
        e2 = evaluate!(mach2,
                       resampling=CV(),
                       measure=misclassification_rate,
                       weights=w,
                       operation=predict_mode,
                       rows=rows,
                       acceleration=accel,
                       verbosity=verb)

        @test e1.per_fold ≈ e2.per_fold
    end

    for cache in [true, false]

        # resampler as machine with evaluation weights not specified:
        resampler = Resampler(model=model, resampling=CV();
                              measure=misclassification_rate,
                              operation=predict_mode,
                              cache=cache)
        resampling_machine = machine(resampler, X, y, w, cache=false)
        fit!(resampling_machine, verbosity=verb)
        e1 = evaluate(resampling_machine).measurement[1]
        mach = machine(model, X, y, w, cache=cache)
        e2 = evaluate!(mach, resampling=CV();
                       measure=misclassification_rate,
                       operation=predict_mode,
                       acceleration=accel, verbosity=verb).measurement[1]

        @test e1 ≈ e2

        # resampler as machine with evaluation weights specified:
        weval = rand(rng,3N);
        resampler = Resampler(model=model, resampling=CV();
                              measure=misclassification_rate,
                              operation=predict_mode,
                              weights=weval, acceleration=accel,
                              cache=cache)
        resampling_machine = machine(resampler, X, y, w, cache=false)
        fit!(resampling_machine, verbosity=verb)
        e1   = evaluate(resampling_machine).measurement[1]
        mach = machine(model, X, y, w, cache=cache)
        e2   = evaluate!(mach, resampling=CV();
                         measure=misclassification_rate,
                         operation=predict_mode,
                         weights=weval,
                         acceleration=accel,
                         verbosity=verb).measurement[1]

        @test e1 ≈ e2

    end

    x = [1,2,3,4,5,6,7]
    X, y = Tables.table([x x x x x x]), coerce([1,2,1,3,1,2,2], Multiclass)
    model = Models.DeterministicConstantClassifier()
    class_w = Dict(zip(levels(y), rand(length(levels(y)))))

    for cache in [true, false]

        #resampler as a machine with class weights specified
        cweval = Dict(zip(levels(y), rand(length(levels(y)))));
        resampler = Resampler(model=model, resampling=CV(nfolds=2);
                              measure=MulticlassFScore(return_type=Vector),
                              class_weights=cweval, acceleration=accel)
        resampling_machine = machine(resampler, X, y, cache=false)
        fit!(resampling_machine, verbosity=verb)
        e1   = evaluate(resampling_machine).measurement[1]
        mach = machine(model, X, y, cache=cache)
        e2   = evaluate!(mach, resampling=CV(nfolds=2);
                         measure=MulticlassFScore(return_type=Vector),
                         class_weights=cweval,
                         acceleration=accel,
                         verbosity=verb).measurement[1]

        @test e1 ≈ e2

    end

    @testset "warnings about measures not supporting weights" begin
        model = ConstantClassifier()
        N = 100
        X, y = make_moons(N)
        class_weights = Dict(0=>0.4, 1=>0.6)
        @test_logs((:warn, r"Sample weights"),
                   evaluate(model, X, y,
                            resampling=Holdout(fraction_train=0.5),
                            measure=auc, weights=ones(N)))
        @test_logs((:warn, r"Class weights"),
                   evaluate(model, X, y,
                            resampling=Holdout(fraction_train=0.5),
                            measure=auc, class_weights=class_weights))
    end

end

@testset_accelerated "automatic operations - integration" accel begin
    clf = ConstantClassifier()
    X, y = make_moons(100)
    e1 = evaluate(clf, X, y, resampling=CV(),
                  measures=[LogLoss(), Accuracy()], verbosity=1)
    e2 = evaluate(clf, X, y, resampling=CV(),
                  operation=[predict, predict_mode],
                  measures=[LogLoss(), Accuracy()], verbosity=1)
    @test e1.measurement ≈ e2.measurement
    evaluate(clf, X, y, resampling=CV(),
             operation=predict,
             measures=[LogLoss(), BrierScore()], verbosity=0)
end

#end
true


text:
# MathOptInterface wrapper of Aspasia solver

const DiscreteSet = Union{MOI.Integer, MOI.ZeroOne, MOI.Semiinteger{Float64}}

function _oa_opt(opt::Optimizer)
    if isnothing(opt.oa_opt)
        if isnothing(opt.oa_solver)
            error("No MIP solver specified (set `oa_solver` attribute)")
        end
        opt.oa_opt = MOI.instantiate(opt.oa_solver, with_bridge_type = Float64)
        if !MOI.supports(opt.oa_opt, MOI.LazyConstraintCallback())
            error("MIP solver does not support lazy constraint callbacks")
        end
    end
    return opt.oa_opt
end

MOI.is_empty(opt::Optimizer) = (isnothing(opt.oa_opt) || MOI.is_empty(opt.oa_opt))

MOI.empty!(opt::Optimizer) = _empty(opt)

MOI.get(::Optimizer, ::MOI.SolverName) = "Aspasia"

MOI.supports_incremental_interface(::Optimizer) = true

function MOI.copy_to(opt::Optimizer, src::MOI.ModelLike)
    return MOI.Utilities.default_copy_to(opt, src)
end

# TODO is this OK?
# function MOI.get(opt::Optimizer, attr::MOI.AbstractModelAttribute)
#     return MOI.get(_oa_opt(opt), attr)
# end
function MOI.get(opt::Optimizer, attr::MOI.AnyAttribute, args...)
    return MOI.get(_oa_opt(opt), attr, args...)
end
function MOI.get(model::Optimizer, attr::MOI.AnyAttribute, idxs::Vector)
    return MOI.get.(model, attr, idxs)
end

function MOI.supports(opt::Optimizer, attr::MOI.AnyAttribute, args...)
    return MOI.supports(_oa_opt(opt), attr, args...)
end

function MOI.set(opt::Optimizer, attr::MOI.AnyAttribute, args...)
    return MOI.set(_oa_opt(opt), attr, args...)
end

function MOI.supports_constraint(
    opt::Optimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return MOI.supports_constraint(_oa_opt(opt), F, S)
end

# NOTE slack bridge should handle VAF in cone
# MOI.supports_constraint(::Optimizer, ::Type{<:Union{VV, VAF}}, ::Type{<:ApproxSet}) = true
MOI.supports_constraint(::Optimizer, ::Type{VV}, ::Type{<:ApproxSet}) = true

MOI.add_variable(opt::Optimizer) = MOI.add_variable(_oa_opt(opt))

function MOI.add_constraint(
    opt::Optimizer,
    func::F,
    set::S,
) where {F <: MathOptInterface.AbstractFunction, S <: MathOptInterface.AbstractSet}
    if S <: DiscreteSet
        opt.has_discrete_constraint = true
    end
    if MOI.supports_constraint(_oa_opt(opt), F, S)
        return MOI.add_constraint(_oa_opt(opt), func, set)
    end
    # store for outer approximation
    push!(opt.con_vars, func.variables)
    push!(opt.con_sets, set)
    return MOI.ConstraintIndex{F, S}(length(opt.con_vars))
end

function MOI.get(opt::Optimizer, param::MOI.RawOptimizerAttribute)
    return getproperty(opt, Symbol(param.name))
end

function MOI.supports(::Optimizer, param::MOI.RawOptimizerAttribute)
    return (Symbol(param.name) in fieldnames(Optimizer))
end

function MOI.set(opt::Optimizer, param::MOI.RawOptimizerAttribute, value)
    return setproperty!(opt, Symbol(param.name), value)
end

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(opt::Optimizer, ::MOI.Silent, value::Bool)
    # MOI.set(_oa_opt(opt), MOI.Silent, value)
    opt.verbose = !value
    return
end

MOI.get(opt::Optimizer, ::MOI.Silent) = !opt.verbose

MOI.get(::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION


text:
module TestQuery
using AbstractTables
using Base.Test

include("mytable.jl")

a = [1, 2, 3]
b = [4, 5, 6]
c = ["a", "b", "c"]
d = [:a, :b, :c]
tbl1 = MyTable(
    a = NullableArray(a),
    b = NullableArray(b),
    c = NullableArray(c),
    d = NullableArray(d)
)
_tbl1 = copy(tbl1)

# test select
for (fld, col) in eachcol(tbl1)
    @eval qrya = @query select(tbl1, $fld)
    @eval qryb = @query tbl1 |> select($fld)
    tbl2a = collect(qrya)
    tbl2b = collect(qryb)
    @test isequal(tbl1, _tbl1)
    @test isequal(
        tbl2a,
        @eval MyTable($fld = $col)
    )
    @test isequal(tbl2a, tbl2b)
end

f = [ a[i] * b[i] for i in eachindex(a) ]
tbl2 = MyTable(f = NullableArray(f))
qrya = @query select(tbl1, f = a * b)
qryb = @query tbl1 |> select(f = a * b)
tbl2a = collect(qrya)
tbl2b = collect(qryb)
@test isequal(tbl1, _tbl1)
@test isequal(tbl2, tbl2a)
@test isequal(tbl2a, tbl2b)

# test filter
tbl2 = MyTable(
    a = NullableArray([3]),
    b = NullableArray([6]),
    c = NullableArray(["c"]),
    d = NullableArray([:c])
)
qrya = @query filter(tbl1, a > 2)
qryb = @query tbl1 |> filter(a > 2)
tbl2a = collect(qrya)
tbl2b = collect(qryb)
@test isequal(tbl1, _tbl1)
@test isequal(tbl2, tbl2a)
@test isequal(tbl2a, tbl2b)

# test summarize
tbl2 = MyTable(
    b_avg = NullableArray([mean([4, 5, 6])])
)
qrya = @query summarize(tbl1, b_avg = mean(b))
qryb = @query tbl1 |> summarize(b_avg = mean(b))
tbl2a = collect(qrya)
tbl2b = collect(qryb)
@test isequal(tbl1, _tbl1)
@test isequal(tbl2, tbl2a)
@test isequal(tbl2a, tbl2b)

# test combinations
tbl2 = MyTable(
    b = NullableArray([4])
)
qrya = @query select(filter(tbl1, a == 1), b)
qryb = @query tbl1 |>
    filter(a == 1) |>
    select(b)
tbl2a = collect(qrya)
tbl2b = collect(qryb)
@test isequal(tbl1, _tbl1)
@test isequal(tbl2a, tbl2)
@test isequal(tbl2a, tbl2b)

end


text:
export expected_policy_values

import Base: run

run(agent, env::AbstractEnv, args...) =
    run(DynamicStyle(env), NumAgentStyle(env), agent, env, args...)

function run(
    ::Sequential,
    ::SingleAgent,
    agent::AbstractAgent,
    env::AbstractEnv,
    stop_condition,
    hook::AbstractHook = EmptyHook(),
)

    reset!(env)
    agent(PRE_EPISODE_STAGE, env)
    hook(PRE_EPISODE_STAGE, agent, env)
    action = agent(PRE_ACT_STAGE, env)
    hook(PRE_ACT_STAGE, agent, env, action)

    while true
        env(action)
        agent(POST_ACT_STAGE, env)
        hook(POST_ACT_STAGE, agent, env)

        if get_terminal(env)
            agent(POST_EPISODE_STAGE, env)  # let the agent see the last observation
            hook(POST_EPISODE_STAGE, agent, env)

            stop_condition(agent, env) && break

            reset!(env)
            agent(PRE_EPISODE_STAGE, env)
            hook(PRE_EPISODE_STAGE, agent, env)
            action = agent(PRE_ACT_STAGE, env)
            hook(PRE_ACT_STAGE, agent, env, action)
        else
            stop_condition(agent, env) && break
            action = agent(PRE_ACT_STAGE, env)
            hook(PRE_ACT_STAGE, agent, env, action)
        end
    end
    hook
end

function run(
    ::Sequential,
    ::SingleAgent,
    agent::AbstractAgent,
    env::MultiThreadEnv,
    stop_condition,
    hook::AbstractHook = EmptyHook(),
)

    while true
        reset!(env)
        action = agent(PRE_ACT_STAGE, env)
        hook(PRE_ACT_STAGE, agent, env, action)

        env(action)
        agent(POST_ACT_STAGE, env)
        hook(POST_ACT_STAGE, agent, env)

        if stop_condition(agent, env)
            agent(PRE_ACT_STAGE, env)  # let the agent see the last observation
            break
        end
    end
    hook
end

function run(
    ::Sequential,
    ::MultiAgent,
    agents::Tuple{Vararg{<:AbstractAgent}},
    env::AbstractEnv,
    stop_condition,
    hooks = [EmptyHook() for _ in agents],
)
    @assert length(agents) == get_num_players(env)

    hooks = Dict(get_role(agent) => hook for (agent, hook) in zip(agents, hooks))
    agents = Dict(get_role(agent) => agent for agent in agents)
    reset!(env)

    agent = agents[get_current_player(env)]
    hook = hooks[get_current_player(env)]

    for p in get_players(env)
        agents[p](PRE_EPISODE_STAGE, env)
        hooks[p](PRE_EPISODE_STAGE, agents[p], env)
    end

    action = agent(PRE_ACT_STAGE, env)
    hook(PRE_ACT_STAGE, agent, env, action)

    while true
        env(action)
        agent(POST_ACT_STAGE, env)
        hook(POST_ACT_STAGE, agent, env)

        if get_terminal(env)
            for p in get_players(env)
                agents[p](POST_EPISODE_STAGE, env)
                hooks[p](POST_EPISODE_STAGE, agents[p], env)
            end

            stop_condition(agent, env) && break

            reset!(env)

            for p in get_players(env)
                agents[p](PRE_EPISODE_STAGE, env)
                hooks[p](PRE_EPISODE_STAGE, agents[p], env)
            end

            agent = agents[get_current_player(env)]
            hook = hooks[get_current_player(env)]
            action = agent(PRE_ACT_STAGE, env)
            hook(PRE_ACT_STAGE, agent, env, action)
        else
            stop_condition(agent, env) && break

            agent = agents[get_current_player(env)]
            hook = hooks[get_current_player(env)]
            action = agent(PRE_ACT_STAGE, env)
            hook(PRE_ACT_STAGE, agent, env, action)
        end
    end

    hooks
end

"""
    expected_policy_values(agents, env)

Calculate the expected return of each agent.
"""
function expected_policy_values(agents::Tuple{Vararg{<:AbstractAgent}}, env::AbstractEnv)
    agents = Dict(get_role(agent) => agent for agent in agents)
    values = expected_policy_values(agents, env)
    Dict(zip(get_players(env), values))
end

expected_policy_values(agents::Dict, env::AbstractEnv) = expected_policy_values(
    agents,
    env,
    RewardStyle(env),
    ChanceStyle(env),
    DynamicStyle(env),
)

function expected_policy_values(
    agents::Dict,
    env::AbstractEnv,
    ::TerminalReward,
    ::Union{ExplicitStochastic,Deterministic},
    ::Sequential,
)
    if get_terminal(env)
        [get_reward(env, get_role(agents[p])) for p in get_players(env)]
    elseif get_current_player(env) == get_chance_player(env)
        vals = zeros(length(agents))
        for a::ActionProbPair in get_legal_actions(env)
            vals .+= a.prob .* expected_policy_values(agents, child(env, a))
        end
        vals
    else
        vals = zeros(length(agents))
        probs = get_prob(agents[get_current_player(env)].policy, env)
        actions = get_actions(env)
        for (a, p) in zip(actions, probs)
            if p > 0 #= ignore illegal action =#
                vals .+= p .* expected_policy_values(agents, child(env, a))
            end
        end
        vals
    end
end


text:
using InfrastructureModels
using PowerModels
using PowerSystems
const PM = PowerModels

# required for reducing logging during tests
using Memento

# Suppress warnings during testing.
setlevel!(getlogger(InfrastructureModels), "error")
setlevel!(getlogger(PowerModels), "error")

# required for "with_optimizer" function
using JuMP

# needed for model building (MOI does not currently suppot adding solvers after model creation)
using Ipopt
ipopt_optimizer = with_optimizer(Ipopt.Optimizer, tol=1e-6, print_level=0)

# is this the best way to find a file in a package?
base_dir = dirname(dirname(pathof(PowerSystems)))
case5_data = PM.parse_file(joinpath(base_dir,"data/matpower/case5.m"))
case5_data = InfrastructureModels.replicate(case5_data, 2)

case5_dc_data = PM.parse_file(joinpath(base_dir,"data/matpower/case5_dc.m"))
case5_dc_data = InfrastructureModels.replicate(case5_dc_data, 2)

# TODO: currently JuMP.num_variables is the best we can do to introspect the JuMP model.
#  Ideally this would also test the number of constraints generated

@test try
    pm = PowerSimulations.build_nip_model(case5_data, PM.DCPPowerModel, optimizer=ipopt_optimizer)
    JuMP.num_variables(pm.model) == 34
true finally end

@test try
    pm = PowerSimulations.build_nip_model(case5_data, PM.ACPPowerModel, optimizer=ipopt_optimizer)
    JuMP.num_variables(pm.model) == 96
true finally end

@test try
    pm = PowerSimulations.build_nip_model(case5_data, PM.SOCWRPowerModel, optimizer=ipopt_optimizer)
    JuMP.num_variables(pm.model) == 110
true finally end

# test PowerSimulations type extentions
DCAngleModel = (data::Dict{String,Any}; kwargs...) -> PM.GenericPowerModel(data, PowerSimulations.DCAngleForm; kwargs...)

StandardACModel = (data::Dict{String,Any}; kwargs...) -> PM.GenericPowerModel(data, PowerSimulations.StandardAC; kwargs...)

@test try
    pm = PowerSimulations.build_nip_model(case5_data, DCAngleModel, optimizer=ipopt_optimizer)
    JuMP.num_variables(pm.model) == 34
true finally end

# test PowerSimulations type extentions
@test try
    pm = PowerSimulations.build_nip_model(case5_data, StandardACModel, optimizer=ipopt_optimizer)
    JuMP.num_variables(pm.model) == 96
true finally end

# test models with HVDC line
@test try
    pm = PowerSimulations.build_nip_model(case5_dc_data, PM.DCPPowerModel, optimizer=ipopt_optimizer)
    JuMP.num_variables(pm.model) == 36
true finally end

@test try
    pm = PowerSimulations.build_nip_model(case5_dc_data, DCAngleModel, optimizer=ipopt_optimizer)
    JuMP.num_variables(pm.model) == 48
true finally end

@test try
    base_dir = dirname(dirname(pathof(PowerSystems)))
    include(joinpath(base_dir,"data/data_5bus_pu.jl"))
    PS_struct = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing,  100.0);
    netinjection = PS.instantiate_network(PS.DCAngleForm, PS_struct);
    PM_dict = PowerSimulations.pass_to_pm(PS_struct, netinjection)
    PM_object = PowerSimulations.build_nip_model(PM_dict, PM.DCPPowerModel, optimizer=ipopt_optimizer);
    JuMP.num_variables(PM_object.model) == 384
true finally end


text:
module FV12_plate_examples
using Statistics
using FinEtools
using FinEtools.MeshImportModule
using FinEtools.MeshExportModule
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using PGFPlotsX
using Test

E = 200e3*phun("MPa");
nu = 0.3;
rho = 8000*phun("KG/M^3");
L = 10.00*phun("M"); t = 0.05*phun("M");
nL = 8; nt =2;
neigvs = 14                   # how many eigenvalues
OmegaShift = (10.0*2*pi)^2;
# Fundamental frequency
f_analytical =[0 0 0  0 0 0 1.622 2.360 2.922 4.190 4.190 7.356 7.356 7.668];

function FV12_plate_esnice()
    global a,b,h,na,nb,nh
    MR = DeforModelRed3D
    fens, fes = T4block(L,L,t,nL,nL,nt);

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
    associategeometry!(femm,  geom)

    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")
    println("f/f_analytical = $(fs[7:10]  ./ f_analytical[7:10] .* 100) %")

    vectors = []
    for i = 7:length(fs)
        scattersysvec!(u, v[:, i])
        push!(vectors, ("Mode_$i", deepcopy(u.values)))
    end
    File  =   "rectangular_plate_esnice.vtk"
    vtkexportmesh(File, connasarray(fes), fens.xyz, FinEtools.MeshExportModule.T4; vectors = vectors)
    @async run(`"paraview.exe" $File`)

    true
end # function

function allrun()
    println("#####################################################")
    println("# FV12_plate_esnice ")
    FV12_plate_esnice()
    return true
end # function allrun

end # module FV12_plate_examples