text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module Test

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

using Test

# Be wary of adding new fields to this Config struct. Always think: can it be
# achieved a different way?
mutable struct Config{T<:Real}
    atol::T
    rtol::T
    optimal_status::MOI.TerminationStatusCode
    infeasible_status::MOI.TerminationStatusCode
    exclude::Vector{Any}
end

"""
    Config(
        ::Type{T} = Float64;
        atol::Real = Base.rtoldefault(T),
        rtol::Real = Base.rtoldefault(T),
        optimal_status::MOI.TerminationStatusCode = MOI.OPTIMAL,
        infeasible_status::MOI.TerminationStatusCode = MOI.INFEASIBLE,
        exclude::Vector{Any} = Any[],
    ) where {T}

Return an object that is used to configure various tests.

## Configuration arguments

  * `atol::Real = Base.rtoldefault(T)`: Control the absolute tolerance used
    when comparing solutions.
  * `rtol::Real = Base.rtoldefault(T)`: Control the relative tolerance used
    when comparing solutions.
  * `optimal_status = MOI.OPTIMAL`: Set to `MOI.LOCALLY_SOLVED` if the solver
    cannot prove global optimality.
  * `infeasible_status = MOI.INFEASIBLE`: Set to `MOI.LOCALLY_INFEASIBLE` if the
    solver cannot prove global infeasibility.
  * `exclude = Vector{Any}`: Pass attributes or functions to `exclude` to skip
    parts of tests that require certain functionality. Common arguments include:
     - `MOI.delete` to skip deletion-related tests
     - `MOI.optimize!` to skip optimize-related tests
     - `MOI.ConstraintDual` to skip dual-related tests
     - `MOI.VariableName` to skip setting variable names
     - `MOI.ConstraintName` to skip setting constraint names

## Examples

For a nonlinear solver that finds local optima and does not support finding
dual variables or constraint names:
```julia
Config(
    Float64;
    optimal_status = MOI.LOCALLY_SOLVED,
    exclude = Any[
        MOI.ConstraintDual,
        MOI.VariableName,
        MOI.ConstraintName,
        MOI.delete,
    ],
)
```
"""
function Config(
    ::Type{T} = Float64;
    atol::Real = Base.rtoldefault(T),
    rtol::Real = Base.rtoldefault(T),
    optimal_status::MOI.TerminationStatusCode = MOI.OPTIMAL,
    infeasible_status::MOI.TerminationStatusCode = MOI.INFEASIBLE,
    exclude::Vector{Any} = Any[],
) where {T<:Real}
    return Config{T}(atol, rtol, optimal_status, infeasible_status, exclude)
end

function Base.copy(config::Config{T}) where {T}
    return Config{T}(
        config.atol,
        config.rtol,
        config.optimal_status,
        config.infeasible_status,
        copy(config.exclude),
    )
end

"""
    setup_test(::typeof(f), model::MOI.ModelLike, config::Config)

Overload this method to modify `model` before running the test function `f` on
`model` with `config`. You can also modify the fields in `config` (e.g., to
loosen the default tolerances).

This function should either return `nothing`, or return a function which, when
called with zero arguments, undoes the setup to return the model to its
previous state. You do not need to undo any modifications to `config`.

This function is most useful when writing new tests of the tests for MOI, but it
can also be used to set test-specific tolerances, etc.

See also: [`runtests`](@ref)

## Example

```julia
function MOI.Test.setup_test(
    ::typeof(MOI.Test.test_linear_VariablePrimalStart_partial),
    mock::MOIU.MockOptimizer,
    ::MOI.Test.Config,
)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 0.0]),
    )
    mock.eval_variable_constraint_dual = false

    function reset_function()
        mock.eval_variable_constraint_dual = true
        return
    end
    return reset_function
end
```
"""
setup_test(::Any, ::MOI.ModelLike, ::Config{T}) where {T} = nothing

"""
    version_added(::typeof(function_name))

Returns the version of MOI in which the test `function_name` was added.

This method should be implemented for all new tests.

See the `exclude_tests_after` keyword of [`runtests`](@ref) for more details.
"""
version_added(::F) where {F} = v"0.10.5"  # The default for any unlabeled tests.

"""
    runtests(
        model::MOI.ModelLike,
        config::Config;
        include::Vector{String} = String[],
        exclude::Vector{String} = String[],
        warn_unsupported::Bool = false,
        exclude_tests_after::VersionNumber = v"999.0.0",
    )

Run all tests in `MathOptInterface.Test` on `model`.

## Configuration arguments

 * `config` is a [`Test.Config`](@ref) object that can be used to modify the
   behavior of tests.
 * If `include` is not empty, only run tests that contain an element from
   `include` in their name.
 * If `exclude` is not empty, skip tests that contain an element from `exclude`
   in their name.
 * `exclude` takes priority over `include`.
 * If `warn_unsupported` is `false`, `runtests` will silently skip tests that
   fail with `UnsupportedConstraint` or `UnsupportedAttribute`. When
   `warn_unsupported` is `true`, a warning will be printed. For most cases the
   default behavior (`false`) is what you want, since these tests likely test
   functionality that is not supported by `model`. However, it can be useful to
   run  `warn_unsupported = true` to check you are not skipping tests due to a
   missing `supports_constraint` method or equivalent.
 * `exclude_tests_after` is a version number that excludes any tests to MOI
   added after that version number. This is useful for solvers who can declare a
   fixed set of tests, and not cause their tests to break if a new patch of MOI
   is released with a new test.

See also: [`setup_test`](@ref).

## Example

```julia
config = MathOptInterface.Test.Config()
MathOptInterface.Test.runtests(
    model,
    config;
    include = ["test_linear_"],
    exclude = ["VariablePrimalStart"],
    warn_unsupported = true,
    exclude_tests_after = v"0.10.5",
)
```
"""
function runtests(
    model::MOI.ModelLike,
    config::Config;
    include::Vector{String} = String[],
    exclude::Vector{String} = String[],
    warn_unsupported::Bool = false,
    exclude_tests_after::VersionNumber = v"999.0.0",
)
    for name_sym in names(@__MODULE__; all = true)
        name = string(name_sym)
        if !startswith(name, "test_")
            continue  # All test functions start with test_
        elseif !isempty(include) && !any(s -> occursin(s, name), include)
            continue
        elseif !isempty(exclude) && any(s -> occursin(s, name), exclude)
            continue
        end
        test_function = getfield(@__MODULE__, name_sym)
        if version_added(test_function) > exclude_tests_after
            continue
        end
        @testset "$(name)" begin
            c = copy(config)
            tear_down = setup_test(test_function, model, c)
            # Make sure to empty the model before every test!
            MOI.empty!(model)
            try
                test_function(model, c)
            catch err
                _error_handler(err, name, warn_unsupported)
            end
            if tear_down !== nothing
                tear_down()
            end
        end
    end
    return
end

"""
    RequirementUnmet(msg::String) <: Exception

An error for throwing in tests to indicate that the model does not support some
requirement expected by the test function.
"""
struct RequirementUnmet <: Exception
    msg::String
end

function Base.show(io::IO, err::RequirementUnmet)
    print(io, "RequirementUnmet: $(err.msg)")
    return
end

"""
    @requires(x)

Check that the condition `x` is `true`. Otherwise, throw an [`RequirementUnmet`](@ref)
error to indicate that the model does not support something required by the test
function.

## Examples

```julia
@requires MOI.supports(model, MOI.Silent())
@test MOI.get(model, MOI.Silent())
```
"""
macro requires(x)
    msg = string(x)
    return quote
        if !$(esc(x))
            throw(RequirementUnmet($msg))
        end
    end
end

function _error_handler(
    err::Union{MOI.NotAllowedError,MOI.UnsupportedError,RequirementUnmet},
    name::String,
    warn_unsupported::Bool,
)
    if warn_unsupported
        @warn("Skipping $(name): $(err)")
    end
    return
end

_error_handler(err, ::String, ::Bool) = rethrow(err)

###
### The following are helpful utilities for writing tests in MOI.Test.
###

"""
    Base.isapprox(x, y, config::Config)

A three argument version of `isapprox` for use in MOI.Test.
"""
function Base.isapprox(x, y, config::Config{T}) where {T}
    return Base.isapprox(x, y; atol = config.atol, rtol = config.rtol)
end

"""
    _supports(config::Config, attribute::MOI.AnyAttribute)

Return `true` if the `attribute` is supported by the `config`.

This is helpful when writing tests.

## Example

```julia
if MOI.Test._supports(config, MOI.Silent)
    @test MOI.get(model, MOI.Silent()) == true
end
```
"""
_supports(config::Config, T::Any)::Bool = !(T in config.exclude)

"""
    _test_model_solution(
        model::MOI.ModelLike,
        config::Config;
        objective_value = nothing,
        variable_primal = nothing,
        constraint_primal = nothing,
        constraint_dual = nothing,
    )

Solve, and then test, various aspects of a model.

First, check that `TerminationStatus == config.optimal_status`.

If `objective_value` is not nothing, check that the attribute `ObjectiveValue()`
is approximately `objective_value`.

If `variable_primal` is not nothing, check that the attribute  `PrimalStatus` is
`MOI.FEASIBLE_POINT`. Then for each `(index, value)` in `variable_primal`, check
that the primal value of the variable `index` is approximately `value`.

If `constraint_primal` is not nothing, check that the attribute  `PrimalStatus`
is `MOI.FEASIBLE_POINT`. Then for each `(index, value)` in `constraint_primal`,
check that the primal value of the constraint `index` is approximately `value`.

Finally, if `config.duals = true`, and if `constraint_dual` is not nothing,
check that the attribute  `DualStatus` is `MOI.FEASIBLE_POINT`. Then for each
`(index, value)` in `constraint_dual`, check that the dual of the constraint
`index` is approximately `value`.

### Example

```julia
MOIU.loadfromstring!(model, \"\"\"
    variables: x
    minobjective: 2.0x + 1.0
    c: x >= 1.0
\"\"\")
x = MOI.get(model, MOI.VariableIndex, "x")
c = MOI.get(
    model,
    MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
    "c",
)
_test_model_solution(
    model,
    config;
    objective_value = 3.0,
    variable_primal = [(x, 1.0)],
    constraint_primal = [(c, 1.0)],
    constraint_dual = [(c, 2.0)],
)
```
"""
function _test_model_solution(
    model::MOI.ModelLike,
    config::Config{T};
    objective_value = nothing,
    variable_primal = nothing,
    constraint_primal = nothing,
    constraint_dual = nothing,
) where {T}
    if !_supports(config, MOI.optimize!)
        return
    end
    MOI.optimize!(model)
    # No need to check supports. Everyone _must_ implement ObjectiveValue.
    @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
    if objective_value !== nothing && _supports(config, MOI.ObjectiveValue())
        @test isapprox(
            MOI.get(model, MOI.ObjectiveValue()),
            objective_value,
            config,
        )
    end
    # No need to check supports. Everyone _must_ implement VariablePrimal.
    if variable_primal !== nothing
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        for (index, solution_value) in variable_primal
            @test isapprox(
                MOI.get(model, MOI.VariablePrimal(), index),
                solution_value,
                config,
            )
        end
    end
    if constraint_primal !== nothing && _supports(config, MOI.ConstraintPrimal)
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        for (index, solution_value) in constraint_primal
            @test isapprox(
                MOI.get(model, MOI.ConstraintPrimal(), index),
                solution_value,
                config,
            )
        end
    end
    if constraint_dual !== nothing && _supports(config, MOI.ConstraintDual)
        @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
        for (index, solution_value) in constraint_dual
            @test isapprox(
                MOI.get(model, MOI.ConstraintDual(), index),
                solution_value,
                config,
            )
        end
    end
    return
end

# TODO(odow): The following are helper functions for testing the value types of
#     different attributes. The following attributes are not tested:
# BarrierIterations()
# CallbackNodeStatus()
# ConflictStatus()
# ConstraintBridgingCost()
# ConstraintConflictStatus()
# LazyConstraintCallback()
# NodeCount
# RelativeGap
# SimplexIterations
# VariableBridgingCost()

function _test_attribute_value_type(
    model::MOI.ModelLike,
    attribute::Union{MOI.AbstractModelAttribute,MOI.AbstractOptimizerAttribute},
)
    T = MOI.attribute_value_type(attribute)
    @test MOI.get(model, attribute) isa T
    return
end

function _test_attribute_value_type(
    model::MOI.ModelLike,
    attribute::MOI.AbstractConstraintAttribute,
    ci::MOI.ConstraintIndex,
)
    T = MOI.attribute_value_type(attribute)
    @test @inferred(T, MOI.get(model, attribute, ci)) isa T
    return
end

function _test_attribute_value_type(
    model::MOI.ModelLike,
    attribute::MOI.AbstractVariableAttribute,
    x::MOI.VariableIndex,
)
    T = MOI.attribute_value_type(attribute)
    @test @inferred(T, MOI.get(model, attribute, x)) isa T
    return
end

function _test_variablenames_equal(model, variable_names)
    seen_name = Dict(name => false for name in variable_names)
    for index in MOI.get(model, MOI.ListOfVariableIndices())
        name = MOI.get(model, MOI.VariableName(), index)
        @test haskey(seen_name, name)
        @test seen_name[name] == false
        seen_name[name] = true
    end
    @test all(values(seen_name))
    return
end

function _test_constraintnames_equal(model, constraint_names)
    seen_name = Dict(name => false for name in constraint_names)
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        if F == MOI.VariableIndex
            continue
        end
        for index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            name = MOI.get(model, MOI.ConstraintName(), index)
            @test haskey(seen_name, name)
            @test seen_name[name] == false
            seen_name[name] = true
        end
    end
    @test all(values(seen_name))
    return
end

"""
    util_test_models_equal(
        model1::ModelLike,
        model2::ModelLike,
        variable_names::Vector{String},
        constraint_names::Vector{String},
        single_variable_constraints::Vector{Tuple{String,<:MOI.AbstractScalarSet}}
    )

Test that `model1` and `model2` are identical using `variable_names` as keys for
the variable names and `constraint_names` as keys for the constraint names.

In addition, it checks that there is a VariableIndex-in-Set constraint for each
`(name, set)` tuple in `single_variable_constraints`, where `name` is the name
of the corresponding variable.

!!! warning
    This is not a generic function that works in all cases. It is mainly
    intended for writing tests in which all variables and constraints have
    unique names.
"""
function util_test_models_equal(
    model1::MOI.ModelLike,
    model2::MOI.ModelLike,
    variable_names::Vector{String},
    constraint_names::Vector{String},
    single_variable_constraints::Vector{<:Tuple} = Tuple[],
)
    for (name, set) in single_variable_constraints
        x1 = MOI.get(model1, MOI.VariableIndex, name)
        ci1 = MOI.ConstraintIndex{MOI.VariableIndex,typeof(set)}(x1.value)
        @test MOI.is_valid(model1, ci1)
        @test MOI.get(model1, MOI.ConstraintSet(), ci1) == set
        x2 = MOI.get(model2, MOI.VariableIndex, name)
        ci2 = MOI.ConstraintIndex{MOI.VariableIndex,typeof(set)}(x2.value)
        @test MOI.is_valid(model2, ci2)
        @test MOI.get(model2, MOI.ConstraintSet(), ci2) == set
    end
    _test_variablenames_equal(model1, variable_names)
    _test_variablenames_equal(model2, variable_names)
    _test_constraintnames_equal(model1, constraint_names)
    _test_constraintnames_equal(model2, constraint_names)
    map_2to1 = Dict{MOI.VariableIndex,MOI.VariableIndex}()
    for name in variable_names
        index2 = MOI.get(model2, MOI.VariableIndex, name)
        map_2to1[index2] = MOI.get(model1, MOI.VariableIndex, name)
    end
    for name in constraint_names
        c1 = MOI.get(model1, MOI.ConstraintIndex, name)
        c2 = MOI.get(model2, MOI.ConstraintIndex, name)
        f1 = MOI.get(model1, MOI.ConstraintFunction(), c1)
        f2 = MOI.get(model2, MOI.ConstraintFunction(), c2)
        @test isapprox(f1, MOI.Utilities.map_indices(map_2to1, f2))
        @test MOI.get(model1, MOI.ConstraintSet(), c1) ==
              MOI.get(model2, MOI.ConstraintSet(), c2)
    end
    attrs1 = MOI.get(model1, MOI.ListOfModelAttributesSet())
    attrs2 = MOI.get(model2, MOI.ListOfModelAttributesSet())
    for attr in union(attrs1, attrs2)
        value1 = MOI.get(model1, attr)
        value2 = MOI.get(model2, attr)
        if value1 isa MOI.AbstractFunction
            @test value2 isa MOI.AbstractFunction
            @test isapprox(value1, MOI.Utilities.map_indices(map_2to1, value2))
        else
            @test !(value2 isa MOI.AbstractFunction)
            @test value1 == value2
        end
    end
    return
end

###
### Include all the test files!
###

for file in readdir(@__DIR__)
    if startswith(file, "test_") && endswith(file, ".jl")
        include(file)
    end
end

end # module

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
# Plot recipes

using RecipesBase, ComplexPhasePortrait, PlotUtils

@recipe function f(d::Sampling{2}, u::Matrix{Complex{Float64}})
    x, y = daxes(d)
    let xs = default(:xlims), ys = default(:ylims)
        xlims --> (isnothing(xs) ? (x[1], x[end]) : xs)
        ylims --> (isnothing(ys) ? (y[1], y[end]) : ys)
    end
    yflip := false
    aspect_ratio --> 1
    tick_direction --> :out
    x[:], y[:], transpose(saneportrait(u) .* abs2.(u) / maximum(abs2, u))
end

@recipe function f(d::Sampling{2}, u::Matrix{Float64})
    x, y = daxes(d)
    let xs = default(:xlims), ys = default(:ylims)
        xlims --> (isnothing(xs) ? (x[1], x[end]) : xs)
        ylims --> (isnothing(ys) ? (y[1], y[end]) : ys)
    end
    yflip := false
    aspect_ratio --> 1
    tick_direction --> :out
    if isnothing(default(:clim))
        x[:], y[:], transpose(sense_portrait(u))
    else
        x[:], y[:], transpose(sense_portrait(u, default(:clim)))
    end
end

@recipe f(d::Discretisation{2}, u::Matrix{Float64}) = (d, Complex.(u))

saneportrait(u) = reverse(portrait(u), dims = 1)

"ComplexPhasePortrait, but with real sign instead of phase"
function sense_portrait(xs, mag=maximum(abs, xs))
    C = cgrad([:cyan, :white, :red])
    # TODO stability at |x| ≈ mag
    [C[iszero(mag) ? 1/2 : (x+mag)/2mag] for x in xs]
end

# TODO plot recipes for animations

# BdG spectra

@userplot BdGSpectrum

@recipe function f(uvs::BdGSpectrum)
    # TODO don't require s for J operator
    s, d, ws, us, vs = uvs.args
    J = operators(s, d, :J) |> only
    js = [dot(u, J(u)) + dot(v, J(v)) |> real for (u, v) in zip(us, vs)]
    seriestype := :scatter
    js, ws
end

ext: jl


text:
import DataStructures: OrderedDict

"""
```
HetDSGE{T} <: AbstractHeterogeneousModel{T}
```

### 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.

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

#### Model Specifications and Settings

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

* `subspec::String`: The model subspecification number, indicating that some
  parameters from the original model spec (\"ss0\") 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.

"""
mutable struct HetDSGE{T} <: AbstractHetModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                       # model steady-state values

    # Temporary to get it to work. Need to
    # figure out a more flexible way to define
    # "grids" that are not necessarily quadrature
    # grids within the model
    grids::OrderedDict{Symbol,Union{Grid, Array}}

    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
                                                           # parameters and steady-states

    state_variables::Vector{Symbol}                        # Vector of symbols of the state variables
    jump_variables::Vector{Symbol}                         # Vector of symbols of the jump variables
    normalized_model_states::Vector{Symbol}                # All of the distributional model
                                                           # state variables that need to be normalized
    endogenous_states_unnormalized::OrderedDict{Symbol,UnitRange} # Vector of unnormalized
                                                           # ranges of indices
    endogenous_states::OrderedDict{Symbol,UnitRange}       # Vector of ranges corresponding
                                                           # to normalized (post Klein solution) indices
    exogenous_shocks::OrderedDict{Symbol,Int}              #
    expected_shocks::OrderedDict{Symbol,Int}               #
    equilibrium_conditions::OrderedDict{Symbol,UnitRange}  #
    endogenous_states_augmented::OrderedDict{Symbol, Int}
    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}
end

description(m::HetDSGE) = "HetDSGE, $(m.subspec)"

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

Arguments:
`m:: HetDSGE`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::HetDSGE, states::Vector{Symbol}, jumps::Vector{Symbol})

    endogenous_states = collect(vcat(states, jumps))

    # Exogenous shocks
    exogenous_shocks = collect([:b_sh,:g_sh,:z_sh,:μ_sh,:λ_w_sh, :λ_f_sh,:rm_sh])

    # Equilibrium conditions
    equilibrium_conditions = collect([:eqeuler,:eq_kolmogorov_fwd,
                                      :eq_market_clearing,:eq_lambda, :eq_transfers,
                                      :eq_investment,:eq_tobin_q,:eq_capital_accumulation,
                                      :eq_wage_phillips,:eq_price_phillips,:eq_marginal_cost,
                                      :eq_gdp,:eq_optimal_kl,:eq_taylor, :eq_fisher,
                                      :eq_nominal_wage_inflation,:LR,:LI,:LY,
                                      :LW,:LX,:eq_b,:eq_g,:eq_z,:eq_μ,:eq_λ_w,
                                      :eq_λ_f, :eq_rm]) #, :eq_consumption])

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = [:i_t1, :c_t, :c_t1]

    # Observables
    observables = keys(m.observable_mappings)

    ########################################################################################
    # Setting indices of endogenous_states and equilibrium conditions manually for now

    setup_indices!(m)
    endo = m.endogenous_states_unnormalized
    eqcond = equilibrium_conditions

    ########################################################################################

    m.normalized_model_states = [:kf′_t]

    for (i,k) in enumerate(exogenous_shocks);            m.exogenous_shocks[k]            = i end
    #for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i+length(get_setting(m, :state_indices)) end #length(endogenous_states) end
    for (i,k) in enumerate(observables);                 m.observables[k]                 = i end
end

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

    # Model-specific specifications
    spec               = "het_dsge"
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # initialize empty model
    m = HetDSGE{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(),
            # grids and keys
            OrderedDict{Symbol,Union{Grid, Array}}(), OrderedDict{Symbol,Int}(),

            # normalized_model_states, state_inds, jump_inds
            Vector{Symbol}(), Vector{Symbol}(), Vector{Symbol}(),

            # model indices
            # endogenous states unnormalized, endogenous states normalized
            OrderedDict{Symbol,UnitRange}(), OrderedDict{Symbol,UnitRange}(),
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),
            OrderedDict{Symbol,UnitRange}(),
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),

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

    default_settings!(m)

    # # Set observable transformations
    init_observable_mappings!(m)

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

    # Endogenous states
    states = collect([:kf′_t,:k′_t, :R′_t1,:i′_t1, :y′_t1,:w′_t1,:I′_t1,
                      :b′_t,:g′_t,:z′_t,:μ′_t,:λ_w′_t, :λ_f′_t,:rm′_t]) #, :c′_t1])

    jumps = collect([:l′_t,:R′_t,:i′_t,:t′_t,:w′_t, :L′_t,:π′_t,:π_w′_t,:margutil′_t,:y′_t, :I′_t,
                      :mc′_t,:Q′_t,:capreturn′_t]) #, :c′_t])

    # Initialize model indices
    init_model_indices!(m, states, jumps)

    # Init stuff that keeps track of number of states and jumps and states/jump indices (requires model indices first)
    init_states_and_jumps!(m, states, jumps)

    # Initialize parameters
    init_parameters!(m)

    # Initialize aggregate steady state parameters (necessary for grid construction)
    aggregate_steadystate!(m)

    # Initialize grids
    init_grids!(m)

    # Solve for the steady state
    steadystate!(m)

    # So that the indices of m.endogenous_states reflect the normalization
    normalize_model_state_indices!(m)

    endogenous_states_augmented = [:i_t1, :c_t, :c_t1]
    for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i + first(collect(values(m.endogenous_states))[end]) end
    m <= Setting(:n_model_states_augmented, get_setting(m, :n_model_states) +
                 length(m.endogenous_states_augmented))

    init_subspec!(m)

    return m
end

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

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::HetDSGE)
    ######################################
    # Parameters that affect steady-state
    ######################################
    m <= parameter(:r, 0.01, fixed = true,
                   description = "r: Steady-state real interest rate.", tex_label = "r")
    m <= parameter(:α, 0.3, fixed = true, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), Normal(0.30, 0.05),
                   description = "α: Capital elasticity in the intermediate goods sector's production function (also known as the capital share).",
                   tex_label = "\\alpha")
    # Check this: Previously the description was "Aggregate hours worked"
    m <= parameter(:H, 1.0, (-1000., 1000.), (-1e3, 1e3), Untransformed(),
                   Normal(-45., 5.), fixed = true,
                   description = "Lmean: Mean level of hours", tex_label = "\\bar{L}")
    m <= parameter(:δ, 0.03, fixed = true,
                   description = "δ: The capital depreciation rate", tex_label = "\\delta")
    m <= parameter(:μ_sp, 0.0, fixed = true, description = "μ_sp: The trend in the skill process",
                   tex_label = "\\mu_{sp}")
    m <= parameter(:ρ_sp, 0.7, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = true,
                   description="ρ_sp: AR(1) coefficient in the skill process.",
                   tex_label="\\rho_{sp}")
    m <= parameter(:σ_sp, 0.01, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = true,
                   description="σ_sp: The standard deviation of the skill process.",
                   tex_label="\\sigma_{sp}")

    # Exogenous processes - level
    # Uncomment scaling once adjusted properly in the code
    m <= parameter(:γ, 0.0, (-5.0, 5.0), (-5., 5.), Untransformed(),
                   Normal(0.4, 0.1), fixed = true, # scaling = x -> x/100,
                   description = "γ: The log of the steady-state growth rate of technology",
                   tex_label="100\\gamma")
    # Check this
    m <= parameter(:g, 1/(1-0.01), fixed = true,
                   description = "g_star: 1 - (c_star + i_star)/y_star",
                   tex_label = "g_*")
    # Not in m1002
    m <= parameter(:η, 0.1, description = "η: Borrowing constraint (normalized by TFP)")

    ####################################################
    # Parameters that affect dynamics (not steady-state)
    ####################################################
    m <= parameter(:spp, 4., (-15., 15.), (-15., 15.), Untransformed(),
                   Normal(4., 1.5), fixed = false,
                   description = "S'': The second derivative of households' cost of adjusting investment.",
                   tex_label = "S''")

    m <= parameter(:ϕh, 2., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(),
                   Normal(2, 0.75), fixed = false, description = "inverse frisch elasticity",
                   tex_label = "\\phi_h")

    m <= parameter(:κ_p, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3), fixed = false, description = "κ_p : The slope of the Price Phillips curve", tex_label = "\\kappa_p")
    m <= parameter(:κ_w, 0.5, (1e-5, 5.), (1e-5, 5.), SquareRoot(), GammaAlt(0.5, 0.3), fixed = false, description = "κ_w: The slope of the Wage Phillips curve", tex_label = "\\kappa_w")

    m <= parameter(:ρR , 0.75, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.75, 0.10), fixed = false,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho_R")
    m <= parameter(:ψπ , 1.5, (1e-5, 10.), (1e-5, 10.0), ModelConstructors.Exponential(),
                   Normal(1.5, 0.25), fixed = false,
                   description = "ψ1: Weight on inflation gap in monetary policy rule.",
                   tex_label = "\\psi_\\pi")
    m <= parameter(:ψy , 0.5, (-0.5, 0.5), (-0.5, 0.5), Untransformed(),
                   Normal(0.12, 0.05), fixed = false,
                   description = "ψy: Weight on output growth in monetary policy rule",
                   tex_label = "\\psi_y")

    # Exogenous processes - autocorrelation
    m <= parameter(:ρ_G, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")
    m <= parameter(:ρ_B, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_B")
    m <= parameter(:ρ_μ, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")
    m <= parameter(:ρ_z, 0.5,(1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_z")
    m <= parameter(:ρ_lamf, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")
    m <= parameter(:ρ_lamw, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")
    m <= parameter(:ρ_mon, 0.5, (1e-5, 1 - 1e-5), (1e-5, 1-1e-5), SquareRoot(),
                   BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:σ_g, 0.15, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = false,
                   description = "σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}")
    m <= parameter(:σ_b, 0.15, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = false,
                   description = "σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}")
    m <= parameter(:σ_μ, 0.15, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = false,
                   description = "σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label = "\\sigma_{\\mu}")
    m <= parameter(:σ_z, 0.15,(1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = false,
                   description = "σ_z: The standard deviation of the process describing the stationary component of productivity.",
                   tex_label = "\\sigma_z")
    m <= parameter(:σ_λ_f, 0.15, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = false,
                   description = "σ_λ_f: The mean of the process that generates the price elasticity of the composite good. Specifically, the elasticity is (1+λ_{f,t})/(λ_{f,t}).",
                   tex_label = "\\sigma_{\\lambda_f}")
    m <= parameter(:σ_λ_w, 0.15, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = false,
                   description = "σ_λ_w",
                   tex_label = "\\sigma_{\\lambda_w}")
    m <= parameter(:σ_rm, 0.15, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(),
                   RootInverseGamma(2, 0.10), fixed = false,
                   description = "σ_r_m: The standard deviation of the monetary policy shock.", tex_label = "\\sigma_{r^m}")

    m <= parameter(:π_star, 0.7000, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.62, 0.1), fixed=false, scaling = x -> 1 + x/100,
               description="π_star: The steady-state rate of inflation.",
               tex_label="\\pi_*")

    m <= parameter(:Lmean, -45.9364, (-1000., 1000.), (-1e3, 1e3), Untransformed(), Normal(-45., 5.), fixed=false,
                   description="Lmean: Mean level of hours.",
                   tex_label="\\bar{L}")

    m <= parameter(:e_y, 0.0, fixed = true, description = "e_y: Measurement error on GDP", tex_label = "e_y")
    m <= parameter(:e_L, 0.0, fixed = true, description = "e_L: Measurement error on hours worked", tex_label = "e_L")
    m <= parameter(:e_w, 0.0, fixed = true, description = "e_w: Measurement error on wages", tex_label ="e_w")
    m <= parameter(:e_π, 0.0, fixed = true, description = "e_π: Measurement error on GDP deflator", tex_label = "e_π")
    m <= parameter(:e_R, 0.0, fixed = true, description = "e_R: Measurement error on nominal rate of interest", tex_label = "e_R")
    m <= parameter(:e_c, 0.0, fixed = true, description = "e_c: Measurement error on consumption", tex_label = "e_c")
    m <= parameter(:e_i, 0.0, fixed = true, description = "e_i: Measurement error on investment", tex_label = "e_i")

    m <= parameter(:β_save, NaN, fixed = true, description = "saving the betas per particle")

    # Setting steady-state parameters
    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)

    # Steady state grids for functional/distributional variables and market-clearing discount rate
    m <= SteadyStateParameterGrid(:lstar, fill(NaN, nx*ns), description = "Steady-state expected discounted
                                  marginal utility of consumption", tex_label = "l_*")
    m <= SteadyStateParameterGrid(:cstar, fill(NaN, nx*ns), description = "Steady-state consumption",
                                  tex_label = "c_*")
    m <= SteadyStateParameterGrid(:μstar, fill(NaN, nx*ns), description = "Steady-state cross-sectional density of cash on hand",
                                  tex_label = "\\mu_*")
    m <= SteadyStateParameter(:βstar, NaN, description = "Steady-state discount factor",
                              tex_label = "\\beta_*")
end

"""
```
aggregate_steadystate!(m::HetDSGE)
```
Steady state of aggregate scalar variables (Break out the "analytic" steady-state solution here instead of in steadystate!() since some of these parameters are required to construct the grids that are used in the functional/distributional steady state variable calculation (where init_grids! is called before steadystate! in model construction).
"""
function aggregate_steadystate!(m::HetDSGE)
    m <= SteadyStateParameter(:Rkstar, m[:r] + m[:δ], description = "Rental rate on capital", tex_label = "Rk_*")
    m <= SteadyStateParameter(:ωstar, (m[:α]^(m[:α]/(1-m[:α])))*(1-m[:α])*m[:Rkstar]^(-m[:α]/(1-m[:α])), description = "Real wage", tex_label = "\\omega_*")
    m <= SteadyStateParameter(:klstar, (m[:α]/(1-m[:α]))*(m[:ωstar]/m[:Rkstar])*exp(m[:γ]), description = "Capital/Labor ratio", tex_label = "kl_*")
    m <= SteadyStateParameter(:kstar, m[:klstar]*m[:H], description = "Capital", tex_label = "k_*")
    m <= SteadyStateParameter(:xstar, (1-(1-m[:δ])*exp(-m[:γ]))*m[:kstar], description = "Investment", tex_label = "x_*")
    m <= SteadyStateParameter(:ystar, (exp(-m[:α]*m[:γ])*(m[:kstar]^m[:α])*m[:H]^(1-m[:α])), description = "GDP", tex_label = "y_*")
    m <= SteadyStateParameter(:Tstar, m[:Rkstar]*m[:kstar]*exp(-m[:γ]) - m[:xstar] - (1-(1/m[:g]))*m[:ystar], description = "Net transfer to households", tex_label = "T_*")

    return m
end

"""
```
init_grids!(m::HetDSGE)
```
"""
function init_grids!(m::HetDSGE)

    nx      = get_setting(m, :nx)
    ns      = get_setting(m, :ns)
    λ       = get_setting(m, :λ)

    grids = OrderedDict()

    # Skill grid
    #lsgrid, sprob, sscale = tauchen86(m[:μ_sp].value, m[:ρ_sp].value, m[:σ_sp].value, ns, λ)
    sprob = [[0.9 0.1];[0.1 0.9]]
    sgrid = [0.8;1.2]
    (λs, vs) = eigen(Matrix{Float64}(sprob'))
    order_λs = sortperm(λs, rev = true)
    vs = vs[:,order_λs]
    ss_skill_distr = vs[:,1]/sum(vs[:,1])
    means = ss_skill_distr'*sgrid
    meanz = (get_setting(m, :zhi)+get_setting(m, :zlo))/2.
    sgrid = sgrid/(meanz*means) # so that skills integrate to 1
    sscale = sgrid[2] - sgrid[1]
    swts = (sscale/ns)*ones(ns)
    #sgrid = exp.(lsgrid)
    grids[:sgrid] = Grid(sgrid, swts, sscale)

    # Calculate endogenous grid bounds (bounds are analytic functions of model parameters,
    # which ensure there is no mass in the distributions across x where there shouldn't be)
    # Calculate the lowest possible cash on hand in steady state
    zlo = get_setting(m, :zlo)
    smin = minimum(sgrid)*zlo #exp(m[:μ_sp]/(1-m[:ρ_sp]) - get_setting(m, :λ)*sqrt(m[:σ_sp]^2/(1-m[:ρ_sp])^2))*get_setting(m, :zlo)
    m <= Setting(:xlo, m[:ωstar]*smin*m[:H] - (1+m[:r])*m[:η]*exp(-m[:γ]) + m[:Tstar] + sgrid[1]*m[:ωstar]*m[:H]*0.05, "Lower bound on cash on hand")
    m <= Setting(:xhi, max(get_setting(m, :xlo)*2, get_setting(m, :xlo)+5.), "Upper Bound on cash on hand")
    m <= Setting(:xscale, get_setting(m, :xhi) - get_setting(m, :xlo), "Size of the xgrid")
    xlo     = get_setting(m, :xlo)
    xhi     = get_setting(m, :xhi)
    xscale  = get_setting(m, :xscale)

    # Cash on hand grid
    grids[:xgrid] = Grid(uniform_quadrature(xscale), xlo, xhi, nx, scale = xscale)

    # Markov transition matrix for skill
    grids[:fgrid] = sprob./swts

    # Total grid vectorized across both dimensions
    grids[:sgrid_total] = kron(sgrid, ones(nx))
    grids[:xgrid_total] = kron(ones(ns), grids[:xgrid].points)
    grids[:weights_total] = kron(swts, grids[:xgrid].weights)

    m.grids = grids
end

function model_settings!(m::HetDSGE)
    default_settings!(m)

    # Defaults
    # Data settings for released and conditional data. Default behavior is to set vintage
    # of data to today's date.
    vint = Dates.format(now(), DSGE_DATE_FORMAT)
    m <= Setting(:data_vintage, vint, true, "vint", "Data vintage")

    saveroot = normpath(joinpath(dirname(@__FILE__), "../../../../","save"))
    datapath = normpath(joinpath(dirname(@__FILE__), "../../../../","save","input_data"))

    m <= Setting(:saveroot, saveroot, "Root of data directory structure")
    m <= Setting(:dataroot, datapath, "Input data directory path")

    # Solution method
    m <= Setting(:solution_method, :klein)

    # Likelihood method
    m <= Setting(:use_chand_recursion, true)

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

    # Whether one wishes to re-compute βstar
    m <= Setting(:use_last_βstar, false, "Flag to avoid recomputing steady-state.")

    # Number of states and jumps
    m <= Setting(:normalize_distr_variables, true, "Whether or not to perform the
                 normalization of the distributional states in the Klein solution step")
    m <= Setting(:n_predetermined_variables, 0, "Number of predetermined variables after
                 removing the densities. Calculated with the Jacobian normalization.
                 This set to 0 is just the default setting, since it will always be
                 overwritten once the Jacobian is calculated.")

    # Mollifier setting parameters
    m <= Setting(:In, 0.443993816237631, "Normalizing constant for the mollifier")
    m <= Setting(:zlo, 1.0, "Lower bound on second income shock to mollify actual income")
    m <= Setting(:zhi, 5.0, "Upper bound on second income shock to mollify actual income")

    # s: Skill Distribution/ "Units of effective labor" Grid Setup
    m <= Setting(:ns, 2, "Skill distribution grid points")
    m <= Setting(:λ, 2.0, "The λ parameter in the Tauchen distribution calculation")

    # x: Cash on Hand Grid Setup
    m <= Setting(:nx, 300, "Cash on hand distribution grid points")

    m <= Setting(:trunc_distr, true)
    m <= Setting(:rescale_weights, true)
    m <= Setting(:mindens, 1e-8)

    # Function-valued variables include distributional variables
    m <= Setting(:n_function_valued_backward_looking_states, 1, "Number of function-valued
                 backward looking state variables")
    m <= Setting(:n_backward_looking_distributional_vars, 1, "Number of state variables that are
                 distributional variables.")
    m <= Setting(:n_function_valued_jumps, 1, "Number of function-valued jump variables")
    m <= Setting(:n_jump_distributional_vars, 1, "Number of jump variables that are distributional
                 variables.")

    # Note, these settings assume normalization.
    # The n degrees of freedom removed depends on the distributions/dimensions
    # of heterogeneity that we have discretized over, in this case,
    # cash on hand and the skill distribution. In general the rule of
    # thumb is, remove one degree of freedom for the first endogenous distribution (cash on
    # hand), then one additional degree of freedom for each exogenous distribution (skill
    # distribution). Multiple endogenous distributions only permit removing a single degree
    # of freedom since it is then non-trivial to obtain the marginal distributions.
    m <= Setting(:n_degrees_of_freedom_removed_state, 2, "Number of degrees of freedom from the distributional variables to remove.")
    m <= Setting(:n_degrees_of_freedom_removed_jump, 0, "Number of degrees of freedom from the distributional variables to remove.")

end

function setup_indices!(m::HetDSGE)
    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)
    endo = m.endogenous_states_unnormalized
    eqconds = m.equilibrium_conditions
    nxns = nx*ns

    # Endogenous function-valued states
    endo[:kf′_t]    = 1:nxns    #  combination of lagged ell function and lagged m function that predicts m

    #endogenous scalar-valued states
    endo[:k′_t]  = nxns+1:nxns+1             # capital –dont get confused with steadystate object K
    endo[:R′_t1] = nxns+2:nxns+2              # lagged real interest rate
    endo[:i′_t1] = nxns+3:nxns+3             # lagged nominal interest rate
    endo[:y′_t1] = nxns+4:nxns+4             # lagged gdp
    endo[:w′_t1] = nxns+5:nxns+5           # lag real wages
    endo[:I′_t1] = nxns+6:nxns+6           # lag investment–don't get this confused with i, the nominal interest rate
    # exogenous scalar-valued states:
    endo[:b′_t]    = nxns+7:nxns+7        # discount factor shock
    endo[:g′_t]    = nxns+8:nxns+8        # govt spending
    endo[:z′_t]  = nxns+9:nxns+9        # tfp growth
    endo[:μ′_t]   = nxns+10:nxns+10        # investment shock
    endo[:λ_w′_t] = nxns+11:nxns+11        # wage markup
    endo[:λ_f′_t] = nxns+12:nxns+12        # price markup
    endo[:rm′_t]  = nxns+13:nxns+13        # monetary policy shock
    #endo[:c′_t1]  = nxns+14:nxns+14        # lagged consumption

    # function-valued jumps
    endo[:l′_t]  = nxns+14:2*nxns+13 # ell function

    #scalar-valued jumps
    endo[:R′_t]  = 2*nxns+14:2*nxns+14        # real interest rate
    endo[:i′_t]  = 2*nxns+15:2*nxns+15        # nominal interest rate
    endo[:t′_t]  = 2*nxns+16:2*nxns+16        # transfers + dividends
    endo[:w′_t]  = 2*nxns+17:2*nxns+17        # real wage
    endo[:L′_t]  = 2*nxns+18:2*nxns+18        # hours worked
    endo[:π′_t]  = 2*nxns+19:2*nxns+19        # inflation
    endo[:π_w′_t] = 2*nxns+20:2*nxns+20        # nominal wage inflation
    endo[:margutil′_t] =  2*nxns+21:2*nxns+21        # average marginal utility
    endo[:y′_t]  = 2*nxns+22:2*nxns+22       # gdp
    endo[:I′_t]  = 2*nxns+23:2*nxns+23        # investment
    endo[:mc′_t] = 2*nxns+24:2*nxns+24        # marginal cost - this is ζ in HetDSGEₖd.pdf
    endo[:Q′_t]  = 2*nxns+25:2*nxns+25        # Tobin's qfunction
    endo[:capreturn′_t] = 2*nxns+26:2*nxns+26        # return on capital
    #endo[:c′_t] = 2*nxns+27:2*nxns+27        # consumption

    nvars = 2*nxns+26
    nscalars = 26 # num eqs which output scalars
    nyscalars = 13 # num scalar jumps
    nxscalars = nscalars - nyscalars # num scalar states

# create objects needed for solve.jl
# we will order function blocks as follows:
# 1. all function blocks which output a function (first real eqs, then lags)
# 2. all function blocks which map functions to scalars
# 3. all scalar blocks involving endogenous vbls (first real eqs, then lags)
# 4. shock processes
funops = 1:2 # which operators output a function

    # function blocks which output a function
    eqconds[:eq_euler]              = 1:nxns
    eqconds[:eq_kolmogorov_fwd]     = nx*ns+1:2*nx*ns

    # function blocks which map functions to scalars
    eqconds[:eq_market_clearing]    = 2*nxns+1:2*nxns+1
    eqconds[:eq_lambda]             = 2*nxns+2:2*nxns+2
    #scalar blocks involving endogenous variables
    eqconds[:eq_transfers]            = 2*nxns+3:2*nxns+3 # transfers
    eqconds[:eq_investment]          = 2*nxns+4:2*nxns+4 # investment
    eqconds[:eq_tobin_q]              = 2*nxns+5:2*nxns+5 # tobin's q
    eqconds[:eq_capital_accumulation] = 2*nxns+6:2*nxns+6 # capital accumulation
    eqconds[:eq_wage_phillips] = 2*nxns+7:2*nxns+7 # wage phillips curve
    eqconds[:eq_price_phillips] = 2*nxns+8:2*nxns+8 # price phillips curve
    eqconds[:eq_marginal_cost]  = 2*nxns+9:2*nxns+9 # marginal cost
    eqconds[:eq_gdp]  = 2*nxns+10:2*nxns+10 # gdp
    eqconds[:eq_optimal_kl] = 2*nxns+11:2*nxns+11 # optimal K/L ratio
    eqconds[:eq_taylor] = 2*nxns+12:2*nxns+12 # taylor rule
    eqconds[:eq_fisher] = 2*nxns+13:2*nxns+13 # fisher eqn
    eqconds[:eq_nominal_wage_inflation] = 2*nxns+14:2*nxns+14 # nominal wage inflation
    # lagged variables
    eqconds[:LR] = 2*nxns+15:2*nxns+15 # LR
    eqconds[:LI] = 2*nxns+16:2*nxns+16 # LI
    eqconds[:LY] = 2*nxns+17:2*nxns+17 # LY
    eqconds[:LW]  = 2*nxns+18:2*nxns+18 # LW
    eqconds[:LX] = 2*nxns+19:2*nxns+19 # LX
    # shocks
    eqconds[:eq_b] = 2*nxns+20:2*nxns+20 # discount factor B
    eqconds[:eq_g] = 2*nxns+21:2*nxns+21 # govt spending G
    eqconds[:eq_z] = 2*nxns+22:2*nxns+22 # tfp growth Z
    eqconds[:eq_μ] = 2*nxns+23:2*nxns+23 # investment MU
    eqconds[:eq_λ_w] = 2*nxns+24:2*nxns+24 # wage mkup LAMW
    eqconds[:eq_λ_f] = 2*nxns+25:2*nxns+25 # price mkup LAMF
    eqconds[:eq_rm] = 2*nxns+26:2*nxns+26 # monetary policy MON
    #eqconds[:eq_consumption] = 2*nxns+27:2*nxns+27 # monetary policy MON

    # Total grid x*s
    m <= Setting(:n, get_setting(m, :nx) * get_setting(m, :ns), "Total grid size, multiplying
                         across grid dimensions.")

    m <= Setting(:nvars, 2*get_setting(m, :n) + 26, "num variables")
    m <= Setting(:nscalars, 26, " # num eqs which output scalars")
    m <= Setting(:nyscalars, 13, "num scalar jumps")
    m <= Setting(:nxscalars, 26 - 13, "num scalar states")

    m.endogenous_states = deepcopy(endo)
end

function init_states_and_jumps!(m::AbstractModel, states::Vector{Symbol}, jumps::Vector{Symbol})
    endo = m.endogenous_states_unnormalized

    m <= Setting(:states, states)
    m <= Setting(:jumps, jumps)

    state_indices = stack_indices(endo, states)
    jump_indices = stack_indices(endo, jumps)

    m <= Setting(:state_indices, state_indices, "Which indices of m.endogenous_states correspond to backward looking state variables")
    m <= Setting(:jump_indices, jump_indices, "Which indices of m.endogenous_states correspond to jump
                 variables")

    n_dof_removed_state = get_setting(m, :n_degrees_of_freedom_removed_state)
    n_dof_removed_jump = get_setting(m, :n_degrees_of_freedom_removed_jump)

    ####################################################
    # Calculating the number of backward-looking states
    ####################################################
    n_backward_looking_distr_vars = get_setting(m, :n_backward_looking_distributional_vars)
    m <= Setting(:backward_looking_states_normalization_factor,
                 n_dof_removed_state*n_backward_looking_distr_vars,
                 "The number of dimensions removed from the backward looking state variables
                    for the normalization.")

    nxns = get_setting(m, :n) #nx)*get_setting(m, :ns)
    n_backward_looking_vars = length(get_setting(m, :state_indices))
    n_backward_looking_function_valued_vars = get_setting(m, :n_function_valued_backward_looking_states)
    n_backward_looking_scalar_vars = n_backward_looking_vars - nxns*n_backward_looking_function_valued_vars  #n_backward_looking_function_valued_vars

    m <= Setting(:n_backward_looking_states, nxns*n_backward_looking_distr_vars +
                 n_backward_looking_scalar_vars - get_setting(m, :backward_looking_states_normalization_factor),
                 "Number of state variables, in the true sense (fully
                  backward looking) accounting for the discretization across the grid
                  of function-valued variables and the normalization of
                  distributional variables.")

    ##################################
    # Calculating the number of jumps
    ##################################
    n_jump_distr_vars = get_setting(m, :n_jump_distributional_vars)
    m <= Setting(:jumps_normalization_factor,
                 n_dof_removed_jump*n_jump_distr_vars, "The number of dimensions removed from the
                 jump variables for the normalization.")

    n_jump_vars = length(get_setting(m, :jump_indices))
    n_jump_function_valued_vars = get_setting(m, :n_function_valued_jumps)
    n_jump_scalar_vars = n_jump_vars - nxns*n_jump_function_valued_vars

    m <= Setting(:n_jumps, nxns*n_jump_function_valued_vars +
                 n_jump_scalar_vars - get_setting(m, :jumps_normalization_factor),
                 "Number of jump variables (forward looking) accounting for
                  the discretization across the grid of function-valued variables and the
                  normalization of distributional variables.")

    m <= Setting(:n_model_states, get_setting(m, :n_backward_looking_states) + get_setting(m, :n_jumps),
                 "Number of 'states' in the state space model. Because backward and forward
                 looking variables need to be explicitly tracked for the Klein solution
                 method, we have n_states and n_jumps")
    m <= Setting(:n_model_states_augmented, get_setting(m, :n_model_states) +
                 length(m.endogenous_states_augmented))

end


text:
module Result

using ..FEStructure
using ..LoadCase

include("../util/mmp.jl")

include("./nodal_result.jl")
include("./beam_result.jl")
include("./quad_result.jl")
include("./tria_result.jl")

include("./modal_result.jl")

end


text:
using TeleMap
using Test


text:
struct Complex
    r::Real
    i::Real
end

macro len(x)
    return :( sqrt( $x.r * $x.r + $x.i*$x.i ) )
end

print(@len(Complex(3,4)))


text:
import CSV
import ExcelFiles
using DataFrames
using Query

mit = CSV.read("data/1976-2016-president.csv")
usep1 = DataFrame(ExcelFiles.load("data/1980-2014 November General Election.xlsx", "Turnout Rates", skipstartrows=1))
usep2 = DataFrame(ExcelFiles.load("data/2016 November General Election.xlsx", "Turnout Rates", skipstartrows=1))[1:52, :]
# usep3 = DataFrame(ExcelFiles.load("data/2018 November General Election.xlsx", "Turnout Rates", skipstartrows=1))[1:52, :]

for df in [usep1, usep2#=, usep3=#]
    rename!(df, (replace.(string.(names(df)), Ref(r"\s"=>""))))
    rename!(df, "Voting-EligiblePopulation(VEP)" => "VEP")
end

mit_clean = @from i in mit begin
    @where i.year >= 1980
    @select{i.year, state=i.state, i.candidate, i.party, i.candidatevotes}
    # @dropna(:candidate)
    @collect DataFrame
end

mit_total = @from i in mit begin
    @where i.year >= 1980
    @select{i.year, state=i.state, i.totalvotes}
    @collect DataFrame
end
unique!(mit_total)

usep1_clean = @from i in usep1 begin
    @where i.x1 != "United States" && i.Year % 4 == 0
    @select{year=i.Year, state=i.x1, i.VEP}
    @collect DataFrame
end

usep2_clean = @from i in usep2 begin
    @where i.x1 != "United States"
    @select{year=2016, state=i.x1, i.VEP}
    @collect DataFrame
end

#=
usep3_clean = @from i in usep3 begin
    @where i.x1 != "United States"
    @select{state=i.x1, i.VEP}
    @collect DataFrame
end
=#

usep_clean = vcat(usep1_clean, usep2_clean)

df = outerjoin(usep_clean, mit_total, on=[:year, :state], makeunique=true)

CSV.write("data/candidates.csv", mit_clean)
CSV.write("data/totalvotes.csv", df)

max_stars_repo_name: TnTo/Electoral-Flows

max_issues_repo_name: TnTo/Electoral-Flows

max_forks_repo_name: TnTo/Electoral-Flows


text:
# This example solves the linear 1D rotating shallow water equations
# for the depth-integrated velocities u and v and for the
# surface elevation η in a fluid with constant depth h,
# gravitational acceleration g and Coriolis parameter f:
#
# ∂u/∂t - fv = -gh∂η/∂x - D,   (1)
# ∂v/∂t + fu = - D,            (2)
# ∂η/∂t + ∂u/∂x = - D.         (3)
#
# where D indicates a hyperviscous operator of
# the form ν*(-1)^nν*Laplacian^nν, where ν is the viscosity
# coefficient and nν is the order of the operator.
# The code time-steps Equations (1-3) simultaneously.
#
# The initial conditions prescribe a Gaussian bump larger
# than the deformation radius in a fluid at rest, with
# some random noise smaller than the deformation radius
# added on top. The system develops geostrophically-balanced
# jets around the Gaussian bump, while the smaller-scale noise
# propagates around as inertia-gravity waves.

using FourierFlows, PyPlot
using FFTW: rfft, irfft, fft, ifft
using LinearAlgebra: mul!, ldiv!
using Printf
using Random

function calcN!(N, sol, t, s, v, p, g)
  mul!(v.uh, g.rfftplan, v.u)
  mul!(v.vh, g.rfftplan, v.v)
  mul!(v.ηh, g.rfftplan, v.η)
  rhsu = p.f*v.vh - p.g*p.h*im*g.kr.*v.ηh
  rhsv = - p.f*v.uh
  rhsη = - im*g.kr.*v.uh
  N[:, 1] .= rhsu
  N[:, 2] .= rhsv
  N[:, 3] .= rhsη
  dealias!(N, g, g.kralias)
  nothing
end

struct Vars{Aphys, Atrans} <: AbstractVars
        u :: Aphys
        v :: Aphys
        η :: Aphys
       uh :: Atrans
       vh :: Atrans
       ηh :: Atrans
end

function Vars(::Dev, g::AbstractGrid{T}) where {Dev, T}
  @devzeros Dev T g.nx u v η
  @devzeros Dev Complex{T} (g.nkr) uh vh ηh
  Vars(u, v, η, uh, vh, ηh)
end

struct Params{T} <: AbstractParams
       ν :: T         # Hyperviscosity coefficient
      nν :: Int       # Order of the hyperviscous operator
      g  :: T         # Gravitational acceleration
      h  :: T         # Fluid depth
      f  :: T         # Coriolis parameter
end

function Equation(p, g::AbstractGrid{T}) where T
  L = zeros(T, g.nkr, 3)
  D = - p.ν*g.kr.^(2*p.nν)
  L[:, 1] .= D # u
  L[:, 2] .= D # v
  L[:, 3] .= D # η
  FourierFlows.Equation(L, calcN!, g)
end

"""
    updatevars!(prob)
Update the vars in v on the grid g with the solution in sol.
"""
function updatevars!(prob)
  v, g, sol, p = prob.vars, prob.grid, prob.sol, prob.params
  ldiv!(v.u, g.rfftplan, deepcopy(sol[:, 1]))
  ldiv!(v.v, g.rfftplan, deepcopy(sol[:, 2]))
  ldiv!(v.η, g.rfftplan, deepcopy(sol[:, 3]))
  nothing
end

function set_uvη!(u0, v0, η0, prob)
    g, sol, p, v = prob.grid, prob.sol, prob.params, prob.vars
    mul!(v.uh, g.rfftplan, u0)
    mul!(v.vh, g.rfftplan, v0)
    mul!(v.ηh, g.rfftplan, η0)
    sol[:, 1] = v.uh
    sol[:, 2] = v.vh
    sol[:, 3] = v.ηh
    updatevars!(prob)
end

function updateplot!(axs, prob)
  sca(axs[1])
  cla()
  plot(prob.grid.x, prob.vars.η)
  xlabel(L"$x$ [m]")
  ylabel(L"$\eta$ [m]")
  sca(axs[2])
  cla()
  plot(prob.grid.x, prob.vars.u/prob.params.h) # Plot depth-averaged u.
  xlabel(L"$x$ [m]")
  ylabel(L"$\bar{u}$ [m/s]")
  sca(axs[3])
  cla()
  plot(prob.grid.x, prob.vars.v/prob.params.h) # Plot depth-averaged v.
  xlabel(L"$x$ [m]")
  ylabel(L"$\bar{v}$ [m/s]")
  xlim(prob.grid.x[1], prob.grid.x[end])
  subplots_adjust(hspace=0)
  nothing
end

#---
nx = 256     # Resolution
Lx = 20000π  # Domain length
ν = 1e-32    # Viscosity
nν = 1       # Viscosity order (nν = 1 means Laplacian)
dt = 0.01
g = 10.0     # Gravitational acceleration
h = 200.0    # Fluid depth
f = 1e-2     # Coriolis parameter
dev = CPU()

nt = 100000

grd = OneDGrid(nx, Lx)
vars = Vars(dev, grd)
params = Params(ν, nν, g, h, f)
eqn = Equation(params, grd)
stepper = "FilteredRK4"
iplt = 100

# Set up problem.
prob = FourierFlows.Problem(eqn, stepper, dt, grd, vars, params, dev)

# Set up the initial conditions, in this case a large Gaussian bump
# with small-scale perturbations superimposed.
u0 = zeros(grd.nx)
v0 = zeros(grd.nx)
noise = deepcopy(u0)
noise = Random.randn!(noise)/10
s = 1e4
η0 = @. 3*exp(-grd.x.^2/s^2) + noise
set_uvη!(u0, v0, η0, prob)

# Set up figure.
fig, axs = subplots(nrows=3, figsize=(8, 4), sharex=true)
Ld = @sprintf "%d" sqrt(g*h)/f
suptitle(L"Deformation radius $L_d = \sqrt{gh}/f = $"*string(Ld)*" m")

sca(axs[1])
plot(prob.grid.x, prob.vars.η)
xlabel(L"$x$ [m]")
ylabel(L"$\eta$ [m]")

sca(axs[2])
plot(prob.grid.x, prob.vars.u/prob.params.h) # Plot depth-averaged u.
xlabel(L"$x$ [m]")
ylabel(L"$\bar{u}$ [m/s]")

sca(axs[3])
plot(prob.grid.x, prob.vars.v/prob.params.h) # Plot depth-averaged v.
xlabel(L"$x$ [m]")
ylabel(L"$\bar{v}$ [m/s]")
xlim(prob.grid.x[1], prob.grid.x[end])

subplots_adjust(hspace=0)
show()

# Time-step the problem.
for i=1:nt
    println(i)
    if rem(i, iplt)==0
        updateplot!(axs, prob)
    end
    stepforward!(prob)
    updatevars!(prob)
end


text:
## qValues ##

function qValues{T<:AbstractFloat}(pValues::Vector{T}, pi0::T, pfdr::Bool = false)
    validPValues(pValues)
    validPValues([pi0])
    n = length(pValues)
    u = sortperm(pValues)
    v = competerank(pValues) ## ties with 'min'
    if pfdr
        qvals = (pi0 .* n .* pValues) ./ (v .* (1 - (1 - pValues) .^ n))
    else
        qvals = (pi0 .* n .* pValues) ./ v
    end
    qvals[u[n]] = min(qvals[u[n]], 1)
    for i in (n - 1):-1:1
        qvals[u[i]] = min(qvals[u[i]], qvals[u[i + 1]])
    end
    return qvals
end


text:
abstract type State{V} end

#TODO: support call syntax in Julia 1.3
#(s::State)() = value(s)

value(v) = v
value(s::State) = getfield(s, :value)
value(S::Vector{<:State}) = value.(S)

export value

Base.getindex(s::State) = value(s)
Base.adjoint(s::State) = value(s)

store!(s::State, v) = setfield!(s, :value, convert(valuetype(s), unitfy(v, unittype(s))))

Base.getindex(s::State, i) = value(s)[i]
Base.length(s::State) = 1
Base.iterate(s::State) = (s, nothing)
Base.iterate(s::State, i) = nothing
Base.eltype(::Type{<:State{V}}) where V = V

unittype(::Type{V}) where V = nothing
unittype(::Type{V}) where {V<:Number} = Unitful.unit(V)
unittype(::Type{<:State{V}}) where V = unittype(V)
unittype(::Type{<:Vector{V}}) where V = unittype(V)
unittype(::Type{<:Vector{Union{V,Missing}}}) where V = unittype(V)
unittype(::T) where T = unittype(T)

valuetype(::State{V}) where V = V
valuetype(T, ::Nothing) = T
valuetype(T, U::Units) = Unitful.isunitless(U) ? T : Quantity{T, Unitful.dimension(U), typeof(U)}
valuetype(::Type{Array{T,N}}, U::Units) where {T,N} = Array{valuetype(T, U), N}

rateunittype(U::Nothing, T::Units) = T^-1
rateunittype(U::Units, T::Units) = (R = U/T; Unitful.isunitless(R) ? nothing : R)
rateunittype(U::Units, T::Nothing) = U
rateunittype(U::Nothing, T::Nothing) = nothing

timeunittype(U, TU=u"hr") = isnothing(U) ? TU : (Unitful.dimension(U) == Unitful.𝐓) ? U : TU

struct Nounit{S,U}
    state::S
    unit::U
end

nounit(s::State, u::Units) = Nounit(s, u)
nounit(s::State) = Nounit(s, nothing)
value(s::Nounit) = deunitfy(value(s.state), s.unit)

export nounit

struct Not{S}
    state::S
end

Base.:!(s::State) = Not(s)
value(s::Not) = !value(s.state)

struct And{S1,S2}
    state1::S1
    state2::S2
end

struct Or{S1,S2}
    state1::S1
    state2::S2
end

Base.:&(a::Union{State,And,Or,Not}, b::Union{State,And,Or,Not}) = And(a, b)
Base.:|(a::Union{State,And,Or,Not}, b::Union{State,And,Or,Not}) = Or(a, b)
value(s::And) = value(s.state1) && value(s.state2)
value(s::Or) = value(s.state1) || value(s.state2)

mutable struct StateRef{V,S<:State{V}} <: State{V}
    state::S
end

Base.getindex(r::StateRef) = r.state
value(r::StateRef) = value(r[])
setvalue!(r::StateRef, s::State) = (r.state = s)

#HACK: swap out state variable of mutable System after initialization
setvar!(s::System, k::Symbol, v::State) = begin
    r = s[k]
    @assert r isa StateRef
    setvalue!(r, v)
    a = Dict(fieldnamesalias(s))[k]
    !isnothing(a) && setvalue!(s[a], v)
    nothing
end

abstract type Priority end
struct PrePriority <: Priority end
struct PostPriority <: Priority end

priority(::S) where {S<:State} = priority(S)
priority(::Type{<:State}) = PostPriority()

Base.show(io::IO, s::State) = begin
    v = value(s)
    maxlength = get(io, :maxlength, nothing)
    r = labelstring(v; maxlength)
    print(io, r)
end
Base.show(io::IO, ::MIME"text/plain", s::State) = print(io, value(s))

include("state/hold.jl")
include("state/wrap.jl")
include("state/advance.jl")
include("state/preserve.jl")
include("state/tabulate.jl")
include("state/interpolate.jl")
include("state/track.jl")
include("state/remember.jl")
include("state/provide.jl")
include("state/drive.jl")
include("state/call.jl")
include("state/integrate.jl")
include("state/accumulate.jl")
include("state/capture.jl")
include("state/flag.jl")
include("state/produce.jl")
include("state/bisect.jl")
include("state/solve.jl")


text:
"""
# plotbounds

Plot regression line and intervals based on Stan samples of coeffficients.

$(SIGNATURES)

### Required arguments
```julia
* `df::DataFrame`                      : DataFrame with observed variables and scaled variables
* `xvar::Symbol`                       : X variable in df
* `yvar::Symbol`                       : Y variable in df
* `dfs::DataFrame`                     : DataFrame with Stan samples
* `linkvars::Vector{Symbol}`           : Initial 2 Symbols are regression coefficients,
                                         3rd - only for :predict bounds indicates σ.
```
### Optional arguments
```julia
* `fnc = link`                         : Link function
* `fig::AbstractString=""`             : File to store the plot. If "", a plot is returned
* `bounds::Vector{Symbol}`             : Bounds to display, see below
* `title::AbstractString=""`           : Title for plot
* `xlab::AbstractString=String(xvar)`  : X axis variable label
* `ylab::AbstractString=String(yvar)`  : Y axis variable label
* `alpha::Float64=0.11`                : Interval value, defaults to [0.045, 0.945]
* `colors::Vector{Symbol}`             : Colors for regions, defaults to [:lightgrey, :grey]
* `stepsize::Float64=0.01`             : Stepsize for boundary accuracy
* `rescale_axis=true`                  : Display using un-standardized scale
```

This method is primarily intended to display 2 regions, typically predicted values and
the quantile or hpdi region around the mean line. This is specified using the `bounds`
keyword argument, e.g. `bounds = [:predicted, :hpdi]`.

For the prediction interval, a 3rd parameter needs to be present in the Stan sample
DataFrame (dfs) containing the σ value. This symbol needs to be added to `linkvars`, e.g.
```julia
linkvars = [:a, :bM, :sigma]
```

For other options, :quantile and :hpdi, two parameters suffice (typically the itercept
and slope parameters).

"""
function plotbounds(
  df::DataFrame,
  xvar::Symbol,
  yvar::Symbol,
  dfs::DataFrame,
  linkvars::Vector{Symbol};
  fnc = link,
  fig::AbstractString="",
  bounds::Vector{Symbol}=[:predicted, :hpdi],
  title::AbstractString="",
  xlab::AbstractString=String(xvar),
  ylab::AbstractString=String(yvar),
  alpha::Float64=0.11,
  colors::Vector{Symbol}=[:lightgrey, :grey],
  stepsize::Float64=0.01,
  rescale_axis=true,
  kwargs...
)

  xbar = mean(df[:, xvar])
  xstd = std(df[:, xvar])
  ybar = mean(df[:, yvar])
  ystd = std(df[:, yvar])

  xvar_s = Symbol(String(xvar)*"_s")
  yvar_s = Symbol(String(yvar)*"_s")

  x_s = minimum(df[:, xvar_s]):stepsize:maximum(df[:, xvar_s])
  y_s = fnc(dfs, linkvars, x_s);

  if rescale_axis
    x = rescale(x_s, xbar, xstd)
    y = [rescale(y_s[i], ybar, ystd) for i in 1:length(x)]
  else
    x = x_s
    y = y_s
    xbar = 0.0
    xstd = 1.0
    ybar = 0.0
    ystd = 1.0
  end

  mu = [mean(y[i]) for i in 1:length(x)]

  p = plot(;xlab=xlab, ylab=ylab, title=title, kwargs...)

  if :predicted in bounds
    predictions = zeros(length(x), nrow(dfs))
    for i in 1:length(x)
      for j in 1:nrow(dfs)
        predictions[i, j] =
          ybar .+ rand(Normal(dfs[j, linkvars[1]] + ystd/xstd * dfs[j, linkvars[2]] * (x[i] - xbar), ystd*dfs[j, linkvars[3]]))
      end
    end
    pred_hpd = [hpdi(predictions[i, :], alpha=alpha) for i in 1:length(x)]
    for i in 1:length(x)
      plot!([x[i], x[i]], pred_hpd[i], color=colors[1], leg=false)
    end
  end

  if :quantile in bounds
    bnds_quantile = [quantile(y[i], [alpha/2, 1-alpha/2]) for i in 1:length(x)]
    for i in 1:length(x)
      plot!([x[i], x[i]], bnds_quantile[i], color=colors[2], leg=false)
    end
  end

  if :hpdi in bounds
    bnds_hpd = [hpdi(y[i], alpha=alpha) for i in 1:length(x)]
    for i in 1:length(x)
      plot!([x[i], x[i]], bnds_hpd[i], color=colors[2], leg=false)
    end
  end

  plot!(x, mu, color=:black)
  if rescale_axis
    scatter!(df[:, xvar], df[:, yvar], leg=false, color=:darkblue)
  else
    scatter!(df[:, xvar_s], df[:, yvar_s], leg=false, color=:darkblue)
  end

  if fig == ""
    return(p)
  else
    savefig(p, fig)
  end

end

"""
# plotbounds

Plot regression line and PI interval with alpha=0.11.

$(SIGNATURES)

### Required arguments
```julia
* `df::DataFrame`                      : DataFrame with observed variables and scaled variables
* `xvar::Symbol`                       : X variable in df
* `yvar::Symbol`                       : Y variable in df
* `nt::NamedTuple`                     : NamedTuple with Stan samples
* `linkvars::Vector{Symbol}`           : Initial 2 Symbols are regression coefficients,
                                         3rd - only for :predict bounds indicates σ.
```
### Optional arguments
```julia
* `fnc = link`                         : Link function
* `fig::AbstractString=""`             : File to store the plot. If "", a plot is returned
* `stepsize::Float64=0.01`             : Stepsize for boundary accuracy
* `rescale_axis=true`                  : Display using un-standardized scale
* `title::AbstractString=""`           : Title for plot
* `lab::AbstractString=""`             : X axis variable label
```

This method is intended to the PI region around the mean line.

The 2nd symbol in `linkvars` needs to match the matrix
```julia
linkvars = [:a, :bM]
```

For other options, :quantile and :hpdi, two parameters suffice (typically the itercept
and slope parameters).

"""
function plotbounds(
    df::DataFrame,
    xvar::Symbol,
    yvar::Symbol,
    nt::NamedTuple,
    linkvars::Vector{Symbol};
    fnc::Function=link,
    fig::AbstractString="",
    stepsize=0.01,
    rescale_axis=true,
    lab::AbstractString="",
    title::AbstractString="",
    kwargs...)

    xbar = mean(df[:, xvar])
    xstd = std(df[:, xvar])
    ybar = mean(df[:, yvar])
    ystd = std(df[:, yvar])

    xvar_s = Symbol(String(xvar)*"_s")
    yvar_s = Symbol(String(yvar)*"_s")

    minx = minimum(df[!, xvar_s]) - 0.1
    maxx = maximum(df[!, xvar_s]) + 0.1
    x_s = minx:stepsize:maxx

    k = size(nt[Symbol(linkvars[2])], 1)
    m = fnc(collect(x_s), k)
    y_s = m * nt[Symbol(linkvars[2])]
    y_s .+= nt[Symbol(linkvars[1])]'
    mul = meanlowerupper(y_s)

    p = plot(;title, kwargs...)
    if rescale_axis
        xrange = rescale(x_s, xbar, xstd)
        mm = rescale(mul.mean, ybar, ystd)
        lm = rescale(mul.lower, ybar, ystd)
        um = rescale(mul.upper, ybar, ystd)
        plot!(p, xrange, mm;
            xlab=String(xvar), ylab=String(yvar),
            ribbon=(mm-lm, um-mm),
            lab)
    else
        plot!(p, xrange, mul.mean;
            xlab=String(xvar), ylab=String(yvar),
            ribbon=(mul.mean-mul.lower, mul.upper-mul.mean),
            lab)
    end

    if fig == ""
        return(p)
    else
        savefig(p, fig)
    end
end

export
  plotbounds


text: const CDM_VERSION = v"5.3.1"


text:
typealias NSE Union(Number, Symbol, Expr)
typealias SE Union(Symbol, Expr)

type OdeModel
  file::String

  states::Dict{String, Float64}
  parameters::Dict{String, Float64}
  odes::Dict{String, NSE}
end

OdeModel(file::String) = parse_model(file::String)


text:
using Interact
using Statistics

export iplot, ispecgram

"""
$(SIGNATURES)
Plots interactive timeseries of the signal.
"""
function iplot(s; fs=framerate(s), t0=0.0, legend=false, size=(800,400), ylabel=nothing, kwargs...)
  fs = inHz(fs)
  s = samples(s)
  if isanalytic(s)
    @warn "Plotting only real part of complex signal"
    s = real.(s)
  end
  siglen = Base.size(s, 1)/fs
  panback = button("<")
  panfwd = button(">")
  pan = widget((0:(length(s)-1))/fs; value=0.0, label="Pan:") |> onchange
  zoomout = button("-")
  zoomin = button("+")
  zoom = widget((100:length(s))/fs; value=1000/fs, label="Zoom:") |> onchange
  #pooling = widget(OrderedDict("minmax"=>orderedextrema, "sample"=>nothing, "min"=>minimum, "max"=>maximum, "mean"=>mean); label="Pooling:")
  ymin = minimum(s)
  ymax = maximum(s)
  f = (t, z) -> begin
    t = round(Int, t*fs)
    z = round(Int, z*fs)
    t1 = max(t - z÷2, 1)
    t2 = min(t1 + z, Base.size(s,1))
    t1 = max(t2 - z, 1)
    ds = ceil(Int, (t2-t1)/(10*size[1]))
    #s1 = @view s[t1:t2,:]
    s1 = s[t1:t2,:]
    plot(signal(s1, fs); size=size, t0=t0+float(t1-1)/fs, downsample=ds, leg=legend, kwargs...)
    #timeseries(s[t1:t2,:]; size=size, fs=fs, t0=t0+float(t1-1)/fs, downsample=ds, pooling=pooling[], legend=legend, kwargs...)
    #ds == 1 || annotate!(t0+(t1+t2)/2/fs, ymin+(ymax-ymin)/20, text("Downsampled by $ds", :red, 8))
    ylabel == nothing || ylabel!(ylabel)
    ylims!(ymin, ymax)
  end
  on(n -> pan[] = max(pan[]-zoom[]/2, 0.0), panback)
  on(n -> pan[] = min(pan[]+zoom[]/2, siglen), panfwd)
  on(n -> zoom[] = min(zoom[]*1.5, siglen), zoomout)
  on(n -> zoom[] = max(zoom[]/1.5, 100/fs), zoomin)
  #on(n -> zoom[] = zoom[], pooling)
  vbox(
    hbox(panback, hskip(5px), panfwd, pan, zoomout, hskip(5px), zoomin, zoom),
    #hbox(pooling),
    map(f, pan, zoom)
  )
end

"""
$(SIGNATURES)
Plots interactive spectrogram of the signal.
"""
function ispecgram(s; fs=framerate(s), nfft=min(div(length(s),8),256), noverlap=div(nfft,2), window=nothing, t0=0.0, size=(800,400), kwargs...)
  @assert Base.size(s,2) == 1 "ispecgram only works with vectors"
  fs = inHz(fs)
  siglen = length(s)/fs
  panback = button("<")
  panfwd = button(">")
  pan = widget(0:(length(s)-1)/fs; value=0.0, label="Pan:") |> onchange
  zoomout = button("-")
  zoomin = button("+")
  zoom = widget(0.1:0.1:(length(s)-1)/fs; value=1.0, label="Zoom:") |> onchange
  pooling = widget(OrderedDict("mean"=>mean, "median"=>median, "min"=>minimum, "max"=>maximum, "sample"=>nothing); label="Pooling:")
  f = (t, z) -> begin
    t = round(Int, t*fs)
    z = max(round(Int, z*fs), nfft)
    t1 = max(t - z÷2, 1)
    t2 = min(t1 + z, Base.size(s,1))
    t1 = max(t2 - z, 1)
    ds = ceil(Int, (t2-t1)/(10*size[1]))
    s1 = @view s[t1:t2]
    specgram(s1; size=size, fs=fs, t0=t0+float(t1-1)/fs, downsample=ds, pooling=pooling[], nfft=nfft, noverlap=noverlap, window=window, kwargs...)
    #ds == 1 || annotate!(t0+(t1+t2)/2/fs, ymin+(ymax-ymin)/20, text("Downsampled by $ds", :red, 8))
    #ylims!(ymin, ymax)
  end
  on(n -> pan[] = max(pan[]-zoom[]/2, 0.0), panback)
  on(n -> pan[] = min(pan[]+zoom[]/2, siglen), panfwd)
  on(n -> zoom[] = min(zoom[]*1.5, siglen), zoomout)
  on(n -> zoom[] = max(zoom[]/1.5, 100/fs), zoomin)
  on(n -> zoom[] = zoom[], pooling)
  vbox(
    hbox(panback, hskip(5px), panfwd, pan, zoomout, hskip(5px), zoomin, zoom, pooling),
    map(f, pan, zoom)
  )
end


text:
const DEFAULT_TARGETS = Dict(
    :all_different => Dict(
        :domains => [domain(1:4) for i in 1:4],
    ),
    :all_equal => Dict(
        :domains => [domain(1:4) for i in 1:4],
    ),
    :all_equal_param => Dict(
        :domains => [domain(8:12) for i in 1:4],
        :param => 10,
    ),
    :dist_different => Dict(
        :domains => [domain(1:4) for i in 1:4],
    ),
    :eq => Dict(
        :domains => [domain(1:10) for i in 1:2],
    ),
    :less_than_param => Dict(
        :domains => [domain(0:4) for i in 1:1],
        :param => 2,
    ),
    :minus_equal_param => Dict(
        :domains => [domain(0:4) for i in 1:2],
        :param => 3,
    ),
    :ordered => Dict(
        :domains => [domain(1:4) for i in 1:4],
    ),
    :sequential_tasks => Dict(
        :domains => [domain(1:4) for i in 1:4],
    ),
    :sum_equal_param => Dict(
        :domains => [domain(1:9) for i in 1:3],
        :param => 15,
    ),
)

const DEFAULT_CONFIG = Dict(
    :lang => :Julia,
    :local_iter => 100,
    :global_iter => 100,
    :search => :complete,
    :metric => hamming,
    :population => 400,
)

const DEFAULT_PATH = joinpath(dirname(pathof(Constraints)), "compositions")

function learn_from_icn(;
    targets = DEFAULT_TARGETS,
    config = DEFAULT_CONFIG,
    path = DEFAULT_PATH,
    )
    for t in targets
        @info "Starting learning for $(t.first)"
        name = "icn_$(t.first)"
        compose_to_file!(
            concept(usual_constraints[t.first]),
            name,
            joinpath(path, "$name.jl");
            domains=t.second[:domains],
            param=get(t.second, :param, nothing),
            local_iter=config[:local_iter],
            global_iter=config[:global_iter],
            search=config[:search],
            metric=config[:metric],
            popSize=config[:population],
            language=config[:lang],
        )
        println()
    end

end


text:
using AlphaStructures
using ViewerGL
GL =  ViewerGL

"""
Return points from file.
"""
function load_points(filename::String)::Matrix
    io = open(filename, "r")
    point = readlines(io)
    close(io)

    b = [tryparse.(Float64,split(point[i], " ")) for i in 1:length(point)]
    V = hcat(b...)
    return V
end

V = load_points("examples/examples3D/model/cat.txt")

# view points
GL.VIEW([
	GL.GLPoints(permutedims(V))
])

# points by columns
filtration = AlphaStructures.alphaFilter(V);
VV, EV, FV, TV = AlphaStructures.alphaSimplex(V, filtration, 0.9)

# view all simplices
GL.VIEW(
	[
		GL.GLGrid(V, EV, GL.COLORS[1], 0.6) # White
		GL.GLGrid(V, FV, GL.COLORS[2], 0.3) # Red
		GL.GLGrid(V, TV, GL.COLORS[3], 0.3) # Green
	]
);

#
# Arlecchino's Cat
#

filter_key = sort(unique(values(filtration)))
granular = 7
reduced_filter = filter_key[sort(abs.(rand(Int, granular).%length(filter_key)))]
reduced_filter = [reduced_filter; max(filter_key...)]
points = [[p] for p in VV]

for α in reduced_filter
	@show α
	VV,EV,FV,TV = AlphaStructures.alphaSimplex(V, filtration, α)
	GL.VIEW(
		GL.GLExplode(
			V,
			[[[t] for t in TV]; [[f] for f in FV]; [[e] for e in EV]],
			1., 1., 1.,	# Explode Ratio
			99, 1		# Colors
		)
	)
end

#
# Appearing Colors
#

#=
for i = 2 : length(reduced_filter)
	VV0, EV0, FV0, TV0 = AlphaStructures.alphaSimplex(V, filtration, reduced_filter[i-1])
	VV,  EV,  FV,  TV = AlphaStructures.alphaSimplex(V, filtration, reduced_filter[i])
	EV0mesh = GL.GLGrid(V, EV0)
	FV0mesh = GL.GLGrid(V, FV0)
	TV0mesh = GL.GLGrid(V, TV0)
	EVmesh = GL.GLGrid(V, setdiff(EV, EV0), GL.COLORS[2], 1)
	FVmesh = GL.GLGrid(V, setdiff(FV, FV0), GL.COLORS[7], 1)
	TVmesh = GL.GLGrid(V, setdiff(TV, TV0), GL.COLORS[12], 1)
	GL.VIEW([EV0mesh; FV0mesh; TV0mesh; EVmesh; FVmesh; TVmesh])
end
=#


text:
module TestPadded

include("preamble.jl")
using UnionArrays.Impl: ofsamesize, unpad, Padded

@testset "ofsamesize" begin
    @testset "bigger = $bigger" for (bigger, value) in [
        (Float64, UInt8(0)),
        (NTuple{7, UInt8}, UInt8(0)),
    ]
        @testset for smaller in [value, typeof(value)]
            padded = ofsamesize(bigger, smaller)
            @test sizeof(padded) == sizeof(bigger)
            @test unpad(padded) === smaller
        end
    end
end

@testset "convert" begin
    xs = ofsamesize.(Float64, UInt8.(1:10)) :: Vector
    @test eltype(xs) <: Padded
    @test unpad(xs[1]) === UInt8(1)
    xs[1] = UInt8(2)
    @test unpad(xs[1]) === UInt8(2)
    xs[1] = ofsamesize(Float64, UInt8(3))
    @test unpad(xs[1]) === UInt8(3)
end

end  # module


text:
#=

  The BUGS Book

  Example 2.7.2 Heart transplant cost-effectiveness (continued)

  Cf ~/jags/bugs_book_2_7_2.R
  """
  model {
     thetaN <- 2 # expected lifetime without transplant
     thetaT ~ dbeta(8,2) # probability of surviving operation
     thetaP ~ dnorm(5,1) # expected survival post-transplant (mean 5, sd 1)
     thetaC ~ dnorm(3000,0.000001)
     # expected cost per year (mean 3000, sd 1000)
     E.c <- (20000 + thetaC*thetaT*thetaP)/1000
     # expected additional cost of transplant
     # in thousands of pounds
     E.e <- thetaT*thetaP - thetaN
     # expected total additional survival
     ICER <- E.c/E.e # incremental cost-effectiveness ratio
     for (i in 1:21) {
       K[i] <- (i-1)*5 # constant!
       INB[i] <- E.e*K[i] - E.c
       Q[i] <- step(INB[i])
     }
  }
  """

  Summary Statistics
  parameters        mean         std   naive_se      mcse          ess      rhat   ess_per_sec
      Symbol     Float64     Float64    Float64   Float64      Float64   Float64       Float64

      thetaT      0.8011      0.1208     0.0012    0.0013   10080.2389    1.0000     2794.6324
      thetaP      5.0099      1.0052     0.0101    0.0107    9915.4268    0.9999     2748.9401
      thetaC   3002.9781    317.8311     3.1783    3.2697    9704.9708    1.0000     2690.5935
         E_c     32.0568      3.3304     0.0333    0.0371   10054.3021    0.9999     2787.4417
         E_e      2.0150      1.0209     0.0102    0.0113   10325.3788    0.9999     2862.5946
        ICER     36.5224   1105.9294    11.0593   11.2014   10007.5129    1.0000     2774.4699

  From The WebPPL model:
  expectation:
  [ 'thetaT', 0.798603971520839 ],
  [ 'thetaP', 5.0233592313165465 ],
  [ 'thetaC', 3002.4139520566914 ],
  [ 'E_c', 32.04848460131377 ],
  [ 'E_e', 2.0127388752617894 ],
  [ 'ICER', 23.878283218113094 ] ]

  Cf ~/webppl/bugs_book_2_7_2.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_2_7_2()

    thetaN = 2.0                          #  expected lifetime without transplant
    thetaT ~ Beta(8,2)                    #  probability of surviving operation
    thetaP ~ Normal(5,1)                  #  expected survival post-transplant (mean 5, sd 1)
    thetaC ~ Normal(3000.0,sqrt(100000))  #  expected cost per year (mean 3000, sd 1000)
    E_c    ~  Dirac((20000.0 + thetaC*thetaT*thetaP)/1000.0) #  expected additional cost of transplant
    # in thousands of pounds
    E_e   ~ Dirac(thetaT*thetaP - thetaN) #  expected total additional survival
    ICER  ~ Dirac(E_c/E_e)                #  incremental cost-effectiveness ratio

    n = 21
    INB = Vector{Real}(undef, n)
    Q = Vector{Real}(undef, n)
    for i in 1:n
       K = (i-1)*5 #  constant!
       INB[i] ~ Dirac(E_e*K - E_c)
       Q[i] ~ Dirac(INB[i] > 0 ? 1 : 0)
    end

end

model = bugs_book_2_7_2()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(10_000), 10_000)
chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)
# chns = sample(model, HMC(0.1,5), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
###
#   Script for launching many Himmelblaus
###

numDims = 1:30

script = ""
for num in numDims

    script =
"#!/bin/bash -l
#SBATCH -p phi
#SBATCH -e errors_$num
#SBATCH -o output_$num
#SBATCH -D ./
#SBATCH --export=ALL
#SBATCH -n 10
#SBATCH -t 24:00:00

module purge
module load julia

#
# Should not need to edit below this line
#
echo ========================================================= echo SLURM job: submitted date = `date`
date_start=`date +%s`
echo ========================================================= echo Job output begins
echo -----------------
echo
hostname

julia Nd_gaussians.jl $num
echo
echo ---------------
echo Job output ends
echo =========================================================
echo SLURM job: finished date = `date`
    "

    #println(script)
    open("runSim_$num.sh","w") do io
        print(io, script)
    end
    run(`sbatch runSim_$num.sh`);
end


text:
using PowerSimulationsDynamics
using Sundials

"""
Case 23:
This case study a 15-state droop grid forming inverter against an infinite bus located at bus 1, with inverter located at bus 2.
The perturbation increase the reference power (analogy for mechanical power) from 0.5 to 0.7.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(TEST_FILES_DIR, "data_tests/test23.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

# time span
tspan = (0.0, 4.0);

# PSCAD benchmark data
csv_file = joinpath(TEST_FILES_DIR, "benchmarks/pscad/Test23/Test23_theta.csv")
t_offset = 9.0

#Define Fault using Callbacks
Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)

@testset "Test 23 Droop Inverter ResidualModel" begin
    path = (joinpath(pwd(), "test-23"))
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            omib_sys, # system
            path,
            tspan,
            Pref_change,
        )

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test23_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test23_eigvals) < 1e-3

        #Solve problem
        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        #Obtain data for angles
        series = get_state_series(results, ("generator-102-1", :θ_oc))
        t = series[1]
        θ = series[2]

        #Obtain PSCAD benchmark data
        M = get_csv_data(csv_file)
        t_pscad = M[:, 1] .- t_offset
        θ_pscad = M[:, 2]

        # Test Transient Simulation Results
        @test LinearAlgebra.norm(θ - θ_pscad) <= 3e-2
        @test LinearAlgebra.norm(t - round.(t_pscad, digits = 3)) == 0.0

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 23 Droop Inverter MassMatrixModel" begin
    path = (joinpath(pwd(), "test-23"))
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            omib_sys, # system
            path,
            tspan,
            Pref_change,
        )

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test23_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test23_eigvals) < 1e-3

        #Solve problem
        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        #Obtain data for angles
        series = get_state_series(results, ("generator-102-1", :θ_oc))
        t = series[1]
        θ = series[2]

        #Obtain PSCAD benchmark data
        M = get_csv_data(csv_file)
        t_pscad = M[:, 1] .- t_offset
        θ_pscad = M[:, 2]

        # Test Transient Simulation Results
        @test LinearAlgebra.norm(θ - θ_pscad) <= 3e-2
        @test LinearAlgebra.norm(t - round.(t_pscad, digits = 3)) == 0.0

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end