text:
const SKIP_PM_VALIDATION = false

"""
    System

    A power system defined by fields for basepower, components, and forecasts.

    # Constructor
    ```julia
    System(basepower)
    System(components, basepower)
    System(buses, generators, loads, branches, storage, basepower, services, annex; kwargs...)
    System(buses, generators, loads, basepower; kwargs...)
    System(file; kwargs...)
    System(; buses, generators, loads, branches, storage, basepower, services, annex, kwargs...)
    System(; kwargs...)
    ```

    # Arguments

    * `buses`::Vector{Bus} : an array of buses
    * `generators`::Vector{Generator} : an array of generators of (possibly) different types
    * `loads`::Vector{ElectricLoad} : an array of load specifications that includes timing of the loads
    * `branches`::Union{Nothing, Vector{Branch}} : an array of branches; may be `nothing`
    * `storage`::Union{Nothing, Vector{Storage}} : an array of storage devices; may be `nothing`
    * `basepower`::Float64 : the base power value for the system
    * `services`::Union{Nothing, Vector{ <: Service}} : an array of services; may be `nothing`

    # Keyword arguments

    * `runchecks`::Bool : run available checks on input fields. If an error is found in a field, that component will not be added to the system and InvalidRange is thrown.
    * `configpath`::String : specify path to validation config file
    DOCTODO: any other keyword arguments? genmap_file, REGEX_FILE
"""
struct System <: PowerSystemType
    data::IS.SystemData
    basepower::Float64             # [MVA]
    bus_numbers::Set{Int}
    runchecks::Bool
    internal::InfrastructureSystemsInternal

    function System(data, basepower, internal; kwargs...)
        bus_numbers = Set{Int}()
        runchecks = get(kwargs, :runchecks, true)
        sys = new(data, basepower, bus_numbers, runchecks, internal)
    end
end

"""Construct an empty System. Useful for building a System while parsing raw data."""
function System(basepower; kwargs...)
    return System(_create_system_data_from_kwargs(; kwargs...), basepower)
end

function System(data, basepower; kwargs...)
    return System(data, basepower, InfrastructureSystemsInternal(); kwargs...)
end

"""System constructor when components are constructed externally."""
function System(buses::Vector{Bus},
                generators::Vector{<:Generator},
                loads::Vector{<:ElectricLoad},
                branches::Union{Nothing, Vector{<:Branch}},
                storage::Union{Nothing, Vector{<:Storage}},
                basepower::Float64,
                services::Union{Nothing, Vector{ <: Service}},
                annex::Union{Nothing,Dict}; kwargs...)

    data = _create_system_data_from_kwargs(; kwargs...)

    sys = System(data, basepower; kwargs...)

    arrays = [buses, generators, loads]
    if !isnothing(branches)
        push!(arrays, branches)
    end
    if !isnothing(storage)
        push!(arrays, storage)
    end
    if !isnothing(services)
        push!(arrays, services)
    end

    error_detected = false

    for component in Iterators.flatten(arrays)
        try
            add_component!(sys, component)
        catch e
            if isa(e, InvalidRange)
                error_detected = true
            else
                rethrow()
            end
        end
    end

    load_zones = isnothing(annex) ? nothing : get(annex, :LoadZones, nothing)
    if !isnothing(load_zones)
        for lz in load_zones
            try
                add_component!(sys, lz)
            catch e
                if isa(e, InvalidRange)
                    error_detected = true
                else
                    rethrow()
                end
            end
        end
    end

    runchecks = get(kwargs, :runchecks, true)

    if error_detected
        throw(InvalidRange("Invalid value(s) detected"))
    end

    if runchecks
        check!(sys)
    end

    return sys
end

"""System constructor without nothing-able arguments."""
function System(buses::Vector{Bus},
                generators::Vector{<:Generator},
                loads::Vector{<:ElectricLoad},
                basepower::Float64; kwargs...)
    return System(buses, generators, loads, nothing, nothing, basepower, nothing, nothing, nothing; kwargs...)
end

"""System constructor with keyword arguments."""
function System(; basepower=100.0,
                buses,
                generators,
                loads,
                branches,
                storage,
                services,
                annex,
                kwargs...)
    return System(buses, generators, loads, branches, storage, basepower, services, annex; kwargs...)
end

"""Constructs a non-functional System for demo purposes."""
function System(::Nothing; buses=[Bus(nothing)],
                generators=[ThermalStandard(nothing), RenewableFix(nothing)],
                loads=[PowerLoad(nothing)],
                branches=nothing,
                storage=nothing,
                basepower=100.0,
                services = nothing,
                annex = nothing,
                kwargs...)
    return System(buses, generators, loads, branches, storage, basepower, services, annex; kwargs...)
end

"""
    to_json(sys::System, filename::AbstractString)

Serializes a system to a JSON string.
"""
function to_json(sys::System, filename::AbstractString)
    IS.prepare_for_serialization!(sys.data, filename)
    return IS.to_json(sys, filename)
end

"""
    to_json(io::IO, sys::System)

Serializes a system an IO stream in JSON.
"""
function to_json(io::IO, sys::System)
    return IS.to_json(io, sys)
end

"""Constructs a System from a JSON file."""
function System(filename::String)
    sys = IS.from_json(System, filename)
    check!(sys)
    return sys
end

"""
    add_component!(sys::System, component::T; kwargs...) where T <: Component

Add a component to the system.

Throws ArgumentError if the component's name is already stored for its concrete type.
Throws ArgumentError if any Component-specific rule is violated.
Throws InvalidRange if any of the component's field values are outside of defined valid
range.

# Examples
```julia
sys = System(100.0)

# Add a single component.
add_component!(sys, bus)

# Add many at once.
buses = [bus1, bus2, bus3]
generators = [gen1, gen2, gen3]
foreach(x -> add_component!(sys, x), Iterators.flatten((buses, generators)))
```
"""
function add_component!(sys::System, component::T; kwargs...) where T <: Component
    check_component_addition(sys, component)

    if Bus in fieldtypes(T)
        check_bus(sys, get_bus(component), component)
    end

    if sys.runchecks && !validate_struct(sys, component)
        throw(InvalidValue("Invalid value for $(component)"))
    end

    IS.add_component!(sys.data, component; kwargs...)

    # Whatever this may change should have been validated above in check_component_addition,
    # so this should not fail.
    handle_component_addition!(sys, component)
end

"""
    add_forecasts!(
                   sys::System,
                   metadata_file::AbstractString,
                   label_mapping::Dict{Tuple{String, String}, String};
                   resolution=nothing,
                  )

Adds forecasts from a metadata file or metadata descriptors.

# Arguments
- `sys::System`: system
- `metadata_file::AbstractString`: metadata file for timeseries
  that includes an array of IS.TimeseriesFileMetadata instances or a vector.
- `label_mapping::Dict{Tuple{String, String}, String}`: maps customized component field names to
  PowerSystem field names
- `resolution::DateTime.Period=nothing`: skip forecast that don't match this resolution.
"""
function add_forecasts!(
                        sys::System,
                        metadata_file::AbstractString,
                        label_mapping::Dict{Tuple{String, String}, String};
                        resolution=nothing,
                       )
    return IS.add_forecasts!(Component, sys.data, metadata_file, label_mapping;
                             resolution=resolution)
end

"""
    add_forecasts!(
                   sys::System,
                   timeseries_metadata::Vector{IS.TimeseriesFileMetadata},
                   resolution=nothing,
                  )

Adds forecasts from a metadata file or metadata descriptors.

# Arguments
- `sys::System`: system
- `timeseries_metadata::Vector{IS.TimeseriesFileMetadata}`: metadata for timeseries
- `resolution::DateTime.Period=nothing`: skip forecast that don't match this resolution.
"""
function add_forecasts!(
                        sys::System,
                        timeseries_metadata::Vector{IS.TimeseriesFileMetadata},
                        resolution=nothing
                       )
    return IS.add_forecasts!(Component, sys.data, timeseries_metadata;
                             resolution=resolution)
end

function IS.add_forecast!(
                          ::Type{<:Component},
                          data::IS.SystemData,
                          forecast_cache::IS.ForecastCache,
                          metadata::IS.TimeseriesFileMetadata;
                          resolution=nothing,
                         )
    IS.set_component!(metadata, data, PowerSystems)
    component = metadata.component
    if isnothing(component)
        return
    end

    forecast, ts_data = IS.make_forecast!(forecast_cache, metadata; resolution=resolution)
    if isnothing(forecast)
        return
    end

    if component isa LoadZones
        uuids = Set([IS.get_uuid(x) for x in get_buses(component)])
        for component_ in (load for load in IS.get_components(ElectricLoad, data)
                          if get_bus(load) |> IS.get_uuid in uuids)
            if forecast isa IS.DeterministicInternal
                forecast_ = IS.DeterministicInternal(IS.get_label(forecast), ts_data)
            else
                # TODO: others
                error("forecast type is not supported yet: $(typeof(forecast))")
            end
            IS.add_forecast!(data, component_, forecast, ts_data)
        end
    else
        IS.add_forecast!(data, component, forecast, ts_data)
    end
end

"""
    iterate_components(sys::System)

Iterates over all components.

# Examples
```julia
for component in iterate_components(sys)
    @show component
end
```

See also: [`get_components`](@ref)
"""
function iterate_components(sys::System)
    return IS.iterate_components(sys.data)
end

"""
    clear_components!(sys::System)

Remove all components from the system.
"""
function clear_components!(sys::System)
    return IS.clear_components!(sys.data)
end

"""
    remove_components!(::Type{T}, sys::System) where T <: Component

Remove all components of type T from the system.

Throws ArgumentError if the type is not stored.
"""
function remove_components!(::Type{T}, sys::System) where T <: Component
    for component in IS.remove_components!(T, sys.data)
        handle_component_removal!(sys, component)
    end
end

"""
    remove_component!(sys::System, component::T) where T <: Component

Remove a component from the system by its value.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(sys::System, component::T) where T <: Component
    IS.remove_component!(sys.data, component)
    handle_component_removal!(sys, component)
end

"""
    remove_component!(
                      ::Type{T},
                      sys::System,
                      name::AbstractString,
                      ) where T <: Component

Remove a component from the system by its name.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(
                           ::Type{T},
                           sys::System,
                           name::AbstractString,
                          ) where T <: Component
    component = IS.remove_component!(T, sys.data, name)
    handle_component_removal!(sys, component)
end

"""
    get_component(
                  ::Type{T},
                  sys::System,
                  name::AbstractString
                 )::Union{T, Nothing} where {T <: Component}

Get the component of concrete type T with name. Returns nothing if no component matches.

See [`get_components_by_name`](@ref) if the concrete type is unknown.

Throws ArgumentError if T is not a concrete type.
"""
function get_component(::Type{T}, sys::System, name::AbstractString) where T <: Component
    return IS.get_component(T, sys.data, name)
end

"""
    get_components(
                   ::Type{T},
                   sys::System,
                  )::FlattenIteratorWrapper{T} where {T <: Component}

Returns an iterator of components. T can be concrete or abstract.
Call collect on the result if an array is desired.

# Examples
```julia
iter = PowerSystems.get_components(ThermalStandard, sys)
iter = PowerSystems.get_components(Generator, sys)
generators = PowerSystems.get_components(Generator, sys) |> collect
generators = collect(PowerSystems.get_components(Generator, sys))
```

See also: [`iterate_components`](@ref)
"""
function get_components(
                        ::Type{T},
                        sys::System,
                       )::FlattenIteratorWrapper{T} where {T <: Component}
    return IS.get_components(T, sys.data)
end

"""
    get_components_by_name(
                           ::Type{T},
                           sys::System,
                           name::AbstractString
                          )::Vector{T} where {T <: Component}

Get the components of abstract type T with name. Note that PowerSystems enforces unique
names on each concrete type but not across concrete types.

See [`get_component`](@ref) if the concrete type is known.

Throws ArgumentError if T is not an abstract type.
"""
function get_components_by_name(
                                ::Type{T},
                                sys::System,
                                name::AbstractString
                               )::Vector{T} where {T <: Component}
    return IS.get_components_by_name(T, sys.data, name)
end

"""
    add_forecast!(sys::System, component::Component, forecast::Forecast)

Adds forecast to the system.

Throws ArgumentError if the component is not stored in the system.

"""
function add_forecast!(sys::System, component::Component, forecast::Forecast)
    return IS.add_forecast!(sys.data, component, forecast)
end

"""
    add_forecast!(
                  sys::System,
                  filename::AbstractString,
                  component::Component,
                  label::AbstractString,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast to a system from a CSV file.

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       filename::AbstractString,
                       component::Component,
                       label::AbstractString,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, filename, component, label, scaling_factor)
end

"""
    add_forecast!(
                  sys::System,
                  ta::TimeSeries.TimeArray,
                  component,
                  label,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast to a system from a TimeSeries.TimeArray.

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       ta::TimeSeries.TimeArray,
                       component,
                       label,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, ta, component, label, scaling_factor)
end

"""
    add_forecast!(
                  sys::System,
                  df::DataFrames.DataFrame,
                  component,
                  label,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast to a system from a DataFrames.DataFrame.

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       df::DataFrames.DataFrame,
                       component,
                       label,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, df, component, label, scaling_factor)
end

"""
    make_forecasts(sys::System, metadata_file::AbstractString; resolution=nothing)

Return a vector of forecasts from a metadata file.

# Arguments
- `data::SystemData`: system
- `metadata_file::AbstractString`: path to metadata file
- `resolution::{Nothing, Dates.Period}`: skip any forecasts that don't match this resolution

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of what the file
should contain.
"""
function make_forecasts(sys::System, metadata_file::AbstractString; resolution=nothing)
    return IS.make_forecasts(sys.data, metadata_file, PowerSystems; resolution=resolution)
end

"""
    make_forecasts(data::SystemData, timeseries_metadata::Vector{TimeseriesFileMetadata};
                   resolution=nothing)

Return a vector of forecasts from a vector of TimeseriesFileMetadata values.

# Arguments
- `data::SystemData`: system
- `timeseries_metadata::Vector{TimeseriesFileMetadata}`: metadata values
- `resolution::{Nothing, Dates.Period}`: skip any forecasts that don't match this resolution
"""
function make_forecasts(sys::System, metadata::Vector{IS.TimeseriesFileMetadata};
                        resolution=nothing)
    return IS.make_forecasts(sys.data, metadata, PowerSystems; resolution=resolution)
end

"""
    generate_initial_times(sys::System, interval::Dates.Period, horizon::Int)

Generates all possible initial times for the stored forecasts. This should return the same
result regardless of whether the forecasts have been stored as one contiguous array or
chunks of contiguous arrays, such as one 365-day forecast vs 365 one-day forecasts.

Throws ArgumentError if there are no forecasts stored, interval is not a multiple of the
system's forecast resolution, or if the stored forecasts have overlapping timestamps.
"""
function generate_initial_times(sys::System, interval::Dates.Period, horizon::Int)
    return IS.generate_initial_times(sys.data, interval, horizon)
end

"""
    get_forecast(
                 ::Type{T},
                 component::Component,
                 initial_time::Dates.DateTime,
                 label::AbstractString,
                ) where T <: Forecast

Return a forecast for the entire time series range stored for these parameters.
"""
function get_forecast(
                      ::Type{T},
                      component::Component,
                      initial_time::Dates.DateTime,
                      label::AbstractString,
                      ) where T <: Forecast
    return IS.get_forecast(T, component, initial_time, label)
end

"""
    get_forecast(
                 ::Type{T},
                 component::Component,
                 initial_time::Dates.DateTime,
                 label::AbstractString,
                 horizon::Int,
                ) where T <: Forecast

Return a forecast for a subset of the time series range stored for these parameters.
"""
function get_forecast(
                      ::Type{T},
                      component::InfrastructureSystemsType,
                      initial_time::Dates.DateTime,
                      label::AbstractString,
                      horizon::Int,
                     ) where T <: Forecast
    return IS.get_forecast(T, component, initial_time, label, horizon)
end

function get_forecast_initial_times(
                                    ::Type{T},
                                    component::Component,
                                   ) where T <: Forecast
    return IS.get_forecast_initial_times(T, component)
end

function get_forecast_initial_times(
                                    ::Type{T},
                                    component::Component,
                                    label::AbstractString
                                   ) where T <: Forecast
    return IS.get_forecast_initial_times(T, component, label)
end

function get_forecast_labels(
                             ::Type{T},
                             component::Component,
                             initial_time::Dates.DateTime,
                            ) where T <: Forecast
    return IS.get_forecast_labels(T, component, initial_time)
end

"""
    get_forecast_values(component::Component, forecast::Forecast)

Return a TimeSeries.TimeArray where the forecast data has been multiplied by the forecasted
component field.
"""
function get_forecast_values(component::Component, forecast::Forecast)
    return IS.get_forecast_values(component, forecast)
end

"""
    get_forecast_initial_times(sys::System)::Vector{Dates.DateTime}

Return sorted forecast initial times.

"""
function get_forecast_initial_times(sys::System)::Vector{Dates.DateTime}
    return IS.get_forecast_initial_times(sys.data)
end

"""
    get_forecast_keys(component::InfrastructureSystemsType)

Return an iterable of NamedTuple keys for forecasts stored for this component.
"""
function get_forecast_keys(component::Component)
    return IS.get_forecast_keys(component)
end

"""
    get_forecasts_horizon(sys::System)

Return the horizon for all forecasts.
"""
function get_forecasts_horizon(sys::System)
    return IS.get_forecasts_horizon(sys.data)
end

"""
    get_forecasts_initial_time(sys::System)

Return the earliest initial_time for a forecast.
"""
function get_forecasts_initial_time(sys::System)
    return IS.get_forecasts_initial_time(sys.data)
end

"""
    get_forecasts_interval(sys::System)

Return the interval for all forecasts.
"""
function get_forecasts_interval(sys::System)
    return IS.get_forecasts_interval(sys.data)
end

"""
    get_forecasts_resolution(sys::System)

Return the resolution for all forecasts.
"""
function get_forecasts_resolution(sys::System)
    return IS.get_forecasts_resolution(sys.data)
end

"""
    iterate_forecasts(sys::System)

Iterates over all forecasts in order of initial time.

# Examples
```julia
for forecast in iterate_forecasts(sys)
    @show forecast
end
```
"""
function iterate_forecasts(sys::System)
    return IS.iterate_forecasts(sys.data)
end

"""
    clear_forecasts!(sys::System)

Remove all forecasts from the system.
"""
function clear_forecasts!(sys::System)
    return IS.clear_forecasts!(sys.data)
end

"""
    check_forecast_consistency(sys::System)

Throws DataFormatError if forecasts have inconsistent parameters.
"""
function check_forecast_consistency(sys::System)
    IS.check_forecast_consistency(sys.data)
end

"""
    validate_forecast_consistency(sys::System)

Return true if all forecasts have consistent parameters.
"""
function validate_forecast_consistency(sys::System)
    return IS.validate_forecast_consistency(sys.data)
end

"""
    remove_forecast!(
                     ::Type{T},
                     sys::System,
                     component::Component,
                     initial_time::Dates.DateTime,
                     label::String,
                    )

Remove the time series data for a component.
"""
function remove_forecast!(
                          ::Type{T},
                          sys::System,
                          component::Component,
                          initial_time::Dates.DateTime,
                          label::String,
                         ) where T <: Forecast
    return IS.remove_forecast!(T, sys.data, component, initial_time, label)
end

"""
    validate_struct(sys::System, value::PowerSystemType)

Validates an instance of a PowerSystemType against System data.
Returns true if the instance is valid.

Users implementing this function for custom types should consider implementing
InfrastructureSystems.validate_struct instead if the validation logic only requires data
contained within the instance.
"""
function validate_struct(sys::System, value::PowerSystemType)::Bool
    return true
end

function check!(sys::System)
    buses = get_components(Bus, sys)
    slack_bus_check(buses)
    buscheck(buses)
end

function JSON2.write(io::IO, sys::System)
    return JSON2.write(io, encode_for_json(sys))
end

function JSON2.write(sys::System)
    return JSON2.write(encode_for_json(sys))
end

function encode_for_json(sys::T) where T <: System
    fields = fieldnames(T)
    final_fields = Vector{Symbol}()
    vals = []

    for field in fields
        # Exclude bus_numbers because they will get rebuilt during deserialization.
        if field != :bus_numbers
            push!(vals, getfield(sys, field))
            push!(final_fields, field)
        end
    end

    return NamedTuple{Tuple(final_fields)}(vals)
end

function JSON2.read(io::IO, ::Type{System})
    raw = JSON2.read(io, NamedTuple)
    data = IS.deserialize(IS.SystemData, Component, raw.data)
    sys = System(data, float(raw.basepower); runchecks=raw.runchecks)
    return sys
end

function IS.deserialize_components(
                                   ::Type{Component},
                                   data::IS.SystemData,
                                   raw::NamedTuple,
                                  )
    # TODO: This adds components through IS.SystemData instead of System, which is what
    # should happen. There is a catch-22 between creating System and SystemData.
    component_cache = Dict{Base.UUID, Component}()

    # Buses and Arcs are encoded as UUIDs.
    composite_components = [Bus]
    for composite_component in composite_components
        for component in IS.get_components_raw(IS.SystemData, composite_component, raw)
            comp = IS.convert_type(composite_component, component)
            IS.add_component!(data, comp)
            component_cache[IS.get_uuid(comp)] = comp
        end
    end

    # Skip Services this round because they have Devices.
    for c_type_sym in IS.get_component_types_raw(IS.SystemData, raw)
        c_type = getfield(PowerSystems, Symbol(IS.strip_module_name(string(c_type_sym))))
        (c_type in composite_components || c_type <: Service) && continue
        for component in IS.get_components_raw(IS.SystemData, c_type, raw)
            comp = IS.convert_type(c_type, component, component_cache)
            IS.add_component!(data, comp)
            component_cache[IS.get_uuid(comp)] = comp
        end
    end

    # Now get the Services.
    for c_type_sym in IS.get_component_types_raw(IS.SystemData, raw)
        c_type = getfield(PowerSystems, Symbol(IS.strip_module_name(string(c_type_sym))))
        if c_type <: Service
            for component in IS.get_components_raw(IS.SystemData, c_type, raw)
                comp = IS.convert_type(c_type, component, component_cache)
                IS.add_component!(data, comp)
            end
        end
    end
end

function JSON2.write(io::IO, component::T) where T <: Component
    return JSON2.write(io, encode_for_json(component))
end

function JSON2.write(component::T) where T <: Component
    return JSON2.write(encode_for_json(component))
end

"""
Encode composite buses as UUIDs.
"""
function encode_for_json(component::T) where T <: Component
    fields = fieldnames(T)
    vals = []

    for name in fields
        val = getfield(component, name)
        if val isa Bus
            push!(vals, IS.get_uuid(val))
        else
            push!(vals, val)
        end
    end

    return NamedTuple{fields}(vals)
end

function IS.convert_type(
                         ::Type{T},
                         data::NamedTuple,
                         component_cache::Dict,
                        ) where T <: Component
    @debug T data
    values = []
    for (fieldname, fieldtype)  in zip(fieldnames(T), fieldtypes(T))
        val = getfield(data, fieldname)
        if fieldtype <: Bus
            uuid = Base.UUID(val.value)
            bus = component_cache[uuid]
            push!(values, bus)
        elseif fieldtype <: Component
            # Recurse.
            push!(values, IS.convert_type(fieldtype, val, component_cache))
        else
            obj = IS.convert_type(fieldtype, val)
            push!(values, obj)
        end
    end

    return T(values...)
end

"""
    get_bus(sys::System, name::String)

Return bus with name.
"""
function get_bus(sys::System, name::String)
    return get_component(Bus, sys, name)
end

"""
    get_bus(sys::System, bus_number::Int)

Return bus with bus_number.
"""
function get_bus(sys::System, bus_number::Int)
    for bus in get_components(Bus, sys)
        if bus.number == bus_number
            return bus
        end
    end

    return nothing
end

"""
    get_buses(sys::System, bus_numbers::Set{Int})

Return all buses values with bus_numbers.
"""
function get_buses(sys::System, bus_numbers::Set{Int})
    buses = Vector{Bus}()
    for bus in get_components(Bus, sys)
        if bus.number in bus_numbers
            push!(buses, bus)
        end
    end

    return buses
end

"""
Throws ArgumentError if a PowerSystems rule blocks addition to the system.

This method is tied with handle_component_addition!. If the methods are re-implemented for
a subtype then whatever is added in handle_component_addition! must be checked here.
"""
function check_component_addition(sys::System, component::Component)
    # no-op
end

"""
Refer to docstring for check_component_addition!
"""
function handle_component_addition!(sys::System, component::Component)
    # no-op
end

function handle_component_removal!(sys::System, component::Component)
    # no-op
end

function check_component_addition(sys::System, branch::Branch)
    arc = get_arc(branch)
    check_bus(sys, get_from(arc), arc)
    check_bus(sys, get_to(arc), arc)
end

function check_component_addition(sys::System, bus::Bus)
    number = get_number(bus)
    if number in sys.bus_numbers
        throw(ArgumentError("bus number $number is already stored in the system"))
    end
end

function handle_component_addition!(sys::System, bus::Bus)
    number = get_number(bus)
    @assert !(number in sys.bus_numbers) "bus number $number is already stored"
    push!(sys.bus_numbers, number)
end

"""
Throws ArgumentError if the bus number is not stored in the system.
"""
function handle_component_removal!(sys::System, bus::Bus)
    number = get_number(bus)
    @assert number in sys.bus_numbers "bus number $number is not stored"
    pop!(sys.bus_numbers, number)
end

"""
    get_bus_numbers(sys::System)

Return a sorted vector of bus numbers in the system.
"""
function get_bus_numbers(sys::System)
    return sort(collect(sys.bus_numbers))
end

"""
Throws ArgumentError if the bus is not stored in the system.
"""
function check_bus(sys::System, bus::Bus, component::Component)
    name = get_name(bus)
    if isnothing(get_component(Bus, sys, name))
        throw(ArgumentError("$component has bus $name that is not stored in the system"))
    end
end

function IS.compare_values(x::System, y::System)::Bool
    match = true

    if !IS.compare_values(x.data, y.data)
        @debug "SystemData values do not match"
        match = false
    end

    if x.basepower != y.basepower
        @debug "basepower does not match" x.basepower y.basepower
        match = false
    end

    return match
end

function _create_system_data_from_kwargs(; kwargs...)
    validation_descriptor_file = nothing
    runchecks = get(kwargs, :runchecks, true)
    if runchecks
        validation_descriptor_file = get(kwargs, :configpath,
                                         POWER_SYSTEM_STRUCT_DESCRIPTOR_FILE)
    end

    return IS.SystemData(; validation_descriptor_file=validation_descriptor_file)
end

function parse_types(mod)
    for name in names(mod)
        mod_type = getfield(mod, name)
        try
            !isstructtype(mod_type) && continue
        catch(e)
            continue
        end
        !isconcretetype(mod_type) && continue
        println("object $mod_type")
    end
    for name in names(mod)
        mod_type = getfield(mod, name)
        !isstructtype(mod_type) && continue
        !isconcretetype(mod_type) && continue
        for (fname, ftype) in zip(fieldnames(mod_type), fieldtypes(mod_type))
            if ftype in names(mod)
                println("$mod_type o-- $ftype")
            end
        end
    end
end


text:
mutable struct CPU
    a::Int64
    b::Int64
    c::Int64
    d::Int64
end

import Base.==
function ==(x::CPU, y::CPU)::Bool
    if x.a == y.a && x.b == y.b && x.c == y.c && x.d == y.d
        return true
    else
        return false
    end
end

struct Op
    instruction::Int64
    a::Int64
    b::Int64
    c::Int64
end

include("util.jl")

function parse_registers(line::String)::CPU
    matches = match(r".*\[(\d+), (\d+), (\d+), (\d+)\]", line)
    return CPU(parse(Int64, matches[1]), parse(Int64, matches[2]), parse(Int64, matches[3]), parse(Int64, matches[4]))
end

function parse_op(line::String)::Op
    matches = match(r"(\d+) (\d+) (\d+) (\d+)", line)
    return Op(parse(Int64, matches[1]), parse(Int64, matches[2]), parse(Int64, matches[3]), parse(Int64, matches[4]))
end

function chunk_input()
    lines = getlines("inputs/16a.txt")
    chunks = []
    while length(lines) > 2
        before = lines[1]
        op = lines[2]
        after = lines[3]
        push!(chunks, (before, op, after))
        deleteat!(lines, 1)
        deleteat!(lines, 1)
        deleteat!(lines, 1)
        try
            deleteat!(lines, 1)
        catch
        end
    end
    parsed_chunks = []
    for chunk in chunks
        before = parse_registers(chunk[1])
        op = parse_op(chunk[2])
        after = parse_registers(chunk[3])
        push!(parsed_chunks, (before, op, after))
    end
    return parsed_chunks
end

function run_instruction(before::CPU, op::Op, i::String)::CPU
    r = [before.a, before.b, before.c, before.d]
    if i == "addr"
        c = CPU(before.a, before.b, before.c, before.d)
        val = r[op.a+1] + r[op.b+1]
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "addi"
        c = CPU(before.a, before.b, before.c, before.d)
        val = r[op.a+1] + op.b
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "mulr"
        c = CPU(before.a, before.b, before.c, before.d)
        val = r[op.a+1] * r[op.b+1]
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "muli"
        c = CPU(before.a, before.b, before.c, before.d)
        val = r[op.a+1] * op.b
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "banr"
        c = CPU(before.a, before.b, before.c, before.d)
        val = r[op.a+1] & r[op.b+1]
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "bani"
        c = CPU(before.a, before.b, before.c, before.d)
        val = r[op.a+1] & op.b
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "borr"
        c = CPU(before.a, before.b, before.c, before.d)
        val = r[op.a+1] | r[op.b+1]
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "bori"
        c = CPU(before.a, before.b, before.c, before.d)
        val = r[op.a+1] | op.b
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "setr"
        c = CPU(before.a, before.b, before.c, before.d)
        val = 0
        if op.a == 0
            val = before.a
        elseif op.a == 1
            val = before.b
        elseif op.a == 2
            val = before.c
        elseif op.a == 3
            val = before.d
        end
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "seti"
        c = CPU(before.a, before.b, before.c, before.d)
        val = op.a
        if op.c == 0
            c.a = val
        elseif op.c == 1
            c.b = val
        elseif op.c == 2
            c.c = val
        elseif op.c == 3
            c.d = val
        end
        return c
    elseif i == "gtir"
        c = CPU(before.a, before.b, before.c, before.d)
        if op.a > r[op.b+1]
            if op.c == 0
                c.a = 1
            elseif op.c == 1
                c.b = 1
            elseif op.c == 2
                c.c = 1
            elseif op.c == 3
                c.d = 1
            end
        else
            if op.c == 0
                c.a = 0
            elseif op.c == 1
                c.b = 0
            elseif op.c == 2
                c.c = 0
            elseif op.c == 3
                c.d = 0
            end
        end
        return c
    elseif i == "gtri"
        c = CPU(before.a, before.b, before.c, before.d)
        if r[op.a+1] > op.b
            if op.c == 0
                c.a = 1
            elseif op.c == 1
                c.b = 1
            elseif op.c == 2
                c.c = 1
            elseif op.c == 3
                c.d = 1
            end
        else
            if op.c == 0
                c.a = 0
            elseif op.c == 1
                c.b = 0
            elseif op.c == 2
                c.c = 0
            elseif op.c == 3
                c.d = 0
            end
        end
        return c
    elseif i == "gtrr"
        c = CPU(before.a, before.b, before.c, before.d)
        if r[op.a+1] > r[op.b+1]
            if op.c == 0
                c.a = 1
            elseif op.c == 1
                c.b = 1
            elseif op.c == 2
                c.c = 1
            elseif op.c == 3
                c.d = 1
            end
        else
            if op.c == 0
                c.a = 0
            elseif op.c == 1
                c.b = 0
            elseif op.c == 2
                c.c = 0
            elseif op.c == 3
                c.d = 0
            end
        end
        return c
    elseif i == "eqir"
        c = CPU(before.a, before.b, before.c, before.d)
        if op.a == r[op.b+1]
            if op.c == 0
                c.a = 1
            elseif op.c == 1
                c.b = 1
            elseif op.c == 2
                c.c = 1
            elseif op.c == 3
                c.d = 1
            end
        else
            if op.c == 0
                c.a = 0
            elseif op.c == 1
                c.b = 0
            elseif op.c == 2
                c.c = 0
            elseif op.c == 3
                c.d = 0
            end
        end
        return c
    elseif i == "eqri"
        c = CPU(before.a, before.b, before.c, before.d)
        if r[op.a+1] == op.b
            if op.c == 0
                c.a = 1
            elseif op.c == 1
                c.b = 1
            elseif op.c == 2
                c.c = 1
            elseif op.c == 3
                c.d = 1
            end
        else
            if op.c == 0
                c.a = 0
            elseif op.c == 1
                c.b = 0
            elseif op.c == 2
                c.c = 0
            elseif op.c == 3
                c.d = 0
            end
        end
        return c
    elseif i == "eqrr"
        c = CPU(before.a, before.b, before.c, before.d)
        if r[op.a+1] == r[op.b+1]
            if op.c == 0
                c.a = 1
            elseif op.c == 1
                c.b = 1
            elseif op.c == 2
                c.c = 1
            elseif op.c == 3
                c.d = 1
            end
        else
            if op.c == 0
                c.a = 0
            elseif op.c == 1
                c.b = 0
            elseif op.c == 2
                c.c = 0
            elseif op.c == 3
                c.d = 0
            end
        end
        return c
    end
end

function num_possible_instructions(before::CPU, op::Op, after::CPU)::Int64
    count = 0
    for i in ["addr", "addi", "mulr", "muli", "banr", "bani", "borr", "bori", "setr", "seti", "gtir", "gtri", "gtrr", "eqir", "eqri", "eqrr"]
        c = run_instruction(before, op, i)
        if c == after
            count += 1
        end
    end
    return count
end

function is_instruction(before::CPU, op::Op, after::CPU, reject::Array{String})::String
    count = 0
    ins = []
    for i in ["addr", "addi", "mulr", "muli", "banr", "bani", "borr", "bori", "setr", "seti", "gtir", "gtri", "gtrr", "eqir", "eqri", "eqrr"]
        if !(i in reject)
            c = run_instruction(before, op, i)
            if c == after
                count += 1
                push!(ins, i)
            end
        end
    end
    if count == 1
        return string("Opcode " * string(op.instruction) * " is " * string(ins[1]))
    end
    return ""
end

function part_1(chunks)::Int64
    count = 0
    for chunk in chunks
        if num_possible_instructions(chunk[1], chunk[2], chunk[3]) >= 3
            count += 1
        end
    end
    return count
end

function part_2(chunks, _opcodes)
    opcodes::Set{String} = Set()
    reject = [o[2] for o in _opcodes]
    for chunk in chunks
        push!(opcodes, is_instruction(chunk[1], chunk[2], chunk[3], reject))
    end
    return opcodes
end

#c = run_instruction(chunks[1][1], chunks[1][2], "addi")
#n = num_possible_instructions(chunks[1][1], chunks[1][2], chunks[1][3])
#p1 = part_1(chunks)
main() = begin
    chunks = chunk_input()
    opcodes::Dict{Int64,String} = Dict()
    opcodes[0] = "addi"
    opcodes[1] = "bani"
    opcodes[2] = "gtir"
    opcodes[3] = "borr"
    opcodes[4] = "eqrr"
    opcodes[5] = "bori"
    opcodes[6] = "gtrr"
    opcodes[7] = "setr"
    opcodes[8] = "muli"
    opcodes[9] = "seti"
    opcodes[10] = "banr"
    opcodes[11] = "gtri"
    opcodes[12] = "eqir"
    opcodes[13] = "eqri"
    opcodes[14] = "addr"
    opcodes[15] = "mulr"
    part_2(chunks, opcodes)
end
main()

function part_2b()
    opcodes::Dict{Int64,String} = Dict()
    opcodes[0] = "addi"
    opcodes[1] = "bani"
    opcodes[2] = "gtir"
    opcodes[3] = "borr"
    opcodes[4] = "eqrr"
    opcodes[5] = "bori"
    opcodes[6] = "gtrr"
    opcodes[7] = "setr"
    opcodes[8] = "muli"
    opcodes[9] = "seti"
    opcodes[10] = "banr"
    opcodes[11] = "gtri"
    opcodes[12] = "eqir"
    opcodes[13] = "eqri"
    opcodes[14] = "addr"
    opcodes[15] = "mulr"
    lines = getlines("inputs/16b.txt")
    ops = []
    for line in lines
        push!(ops, parse_op(line))
    end
    c = CPU(0, 0, 0, 0)
    for op in ops
        println(op)
        c = run_instruction(c, op, opcodes[op.instruction])
        println(c)
    end
end

part_2b()

ext: jl


text:
using Test
using QuantumLattices.Prerequisites.VectorSpaces
using QuantumLattices.Interfaces: dimension, rank

import QuantumLattices.Prerequisites.VectorSpaces: shape, ndimshape
import QuantumLattices.Prerequisites.Traits: contentnames, getcontent

struct SimpleVectorSpace{B, N} <: EnumerativeVectorSpace{B}
    sorted::Bool
    table::NTuple{N, B}
    SimpleVectorSpace(sorted::Bool, table::NTuple{N, B}) where {B, N} = new{B, N}(sorted, table)
end
@inline contentnames(::Type{<:SimpleVectorSpace}) = (:sorted, :table)
@inline Base.issorted(vs::SimpleVectorSpace) = vs.sorted

@testset "EnumerativeVectorSpace" begin
    @test contentnames(EnumerativeVectorSpace) == (:table,)

    id0, id4 = (1, 0), (1, 4)
    id1, id2, id3 = (1, 1), (1, 2), (1, 3)
    vs = SimpleVectorSpace(true, (id1, id2, id3))
    @test vs == deepcopy(vs)
    @test isequal(vs, deepcopy(vs))
    @test vs|>size == (3,)
    @test vs|>dimension == 3
    @test vs|>collect == [id1, id2, id3]
    @test vs[1]==id1 && vs[2]==id2 && vs[3]==id3
    @test searchsortedfirst(vs, id0)==1 && searchsortedfirst(vs, id4)==4
    @test searchsortedfirst(vs, id1)==1 && searchsortedfirst(vs, id2)==2 && searchsortedfirst(vs, id3)==3
    @test isnothing(findfirst(id0, vs)) && isnothing(findfirst(id4, vs))
    @test findfirst(id1, vs)==1 && findfirst(id2, vs)==2 && findfirst(id3, vs)==3
    @test (id0∉vs) && (id4 ∉ vs)
    @test (id1∈vs) && (id2∈vs) && (id3 ∈ vs)

    vs = SimpleVectorSpace(false, (id1, id2, id3))
    @test searchsortedfirst(vs, id0)==4 && searchsortedfirst(vs, id4)==4
    @test searchsortedfirst(vs, id1)==1 && searchsortedfirst(vs, id2)==2 && searchsortedfirst(vs, id3)==3
    @test isnothing(findfirst(id0, vs)) && isnothing(findfirst(id4, vs))
    @test findfirst(id1, vs)==1 && findfirst(id2, vs)==2 && findfirst(id3, vs)==3
end

struct SimpleIndices{N} <: CartesianVectorSpace{CartesianIndex{N}}
    shape::NTuple{N, UnitRange{Int}}
    SimpleIndices(shape::NTuple{N, UnitRange{Int}}) where N = new{N}(shape)
end
SimpleIndices(shape::UnitRange{Int}...) = SimpleIndices(shape)
@inline shape(vs::SimpleIndices) = vs.shape
@inline ndimshape(::Type{<:SimpleIndices{N}}) where N = N
@inline Base.CartesianIndex(basis::CartesianIndex{N}, ::SimpleIndices{N}) where N = basis

@testset "CartesianVectorSpace" begin
    foi = SimpleIndices(2:3, 2:3, 2:3)
    @test dimension(foi) == 8
    @test ndimshape(foi) == ndimshape(typeof(foi)) == 3
    @test issorted(foi) == true
    @test foi|>collect == CartesianIndex.([(2, 2, 2), (3, 2, 2), (2, 3, 2), (3, 3, 2), (2, 2, 3), (3, 2, 3), (2, 3, 3), (3, 3, 3)])
    for (i, index) in enumerate(CartesianIndices((2:3, 2:3, 2:3)))
        @test foi[i] == foi[index] == index
        @test findfirst(index, foi) == i
        @test searchsortedfirst(foi, index) == i
        @test index∈foi
    end
    i1, i2 = CartesianIndex(1, 1, 1), CartesianIndex(4, 4, 4)
    @test isnothing(findfirst(i1, foi)) && isnothing(findfirst(i2, foi))
    @test searchsortedfirst(foi, i1) == 1 && searchsortedfirst(foi, i2) == 9
    @test i1∉foi && i2∉foi
end

@testset "NamedVectorSpace" begin
    nvs = ZipNamedVectorSpace{(:t, :U)}([1, 2], [8.0, 9.0])
    @test nvs==deepcopy(nvs) && isequal(nvs, deepcopy(nvs))
    @test nvs≠ZipNamedVectorSpace{(:t, :V)}([1, 2], [8.0, 9.0])
    @test !isequal(nvs, ZipNamedVectorSpace{(:t, :V)}([1, 2], [8.0, 9.0]))
    @test nvs|>keys == nvs|>typeof|>keys == (:t, :U)
    @test nvs|>values == ([1, 2], [8.0, 9.0])
    @test nvs|>pairs|>collect == [:t=>[1, 2], :U=>[8.0, 9.0]]
    @test eltype(nvs, 1) == eltype(nvs|>typeof, 1) == Int
    @test eltype(nvs, 2) == eltype(nvs|>typeof, 2) == Float64
    @test nvs|>rank == nvs|>typeof|>rank == 1
    @test shape(nvs) == (1:2,)
    @test ndimshape(nvs) == ndimshape(typeof(nvs)) == 1
    elements = [(t = 1, U = 8.0), (t = 2, U = 9.0)]
    for i = 1:dimension(nvs)
        @test NamedTuple(CartesianIndex(elements[i], nvs), nvs) == elements[i]
    end
    @test nvs|>collect == elements

    nvs = ProductNamedVectorSpace{(:t, :U)}([1.0, 2.0], [8.0, 9.0])
    @test nvs|>rank == nvs|>typeof|>rank == 2
    @test shape(nvs) == (1:2, 1:2)
    @test ndimshape(nvs) == ndimshape(typeof(nvs)) == 2
    elements = [(t = 1.0, U = 8.0), (t = 2.0, U = 8.0), (t = 1.0, U = 9.0), (t = 2.0, U = 9.0)]
    for i = 1:dimension(nvs)
        @test NamedTuple(CartesianIndex(elements[i], nvs), nvs) == elements[i]
    end
    @test nvs|>collect == elements
end


text:
using ObjectPools, Test
using ObjectPools: HomogeneousArrayPool, _new!, _recycle!

####
#### tests for internals
####

###
### HomogeneousArrayPools
###

@testset "HomogeneousArrayPool (internals)" begin
    T = Float64
    D = (2,4,3)
    h = HomogeneousArrayPool{T}(D)
    @test h.used == 0
    A1 = _new!(h)
    @test A1 isa Array{T,3} && size(A1) == D
    A2 = _new!(h)
    @test A1 ≢ A2
    @test A2 isa Array{T,3} && size(A2) == D
    @test h.used == 2
    @test _recycle!(h) ≡ nothing
    @test h.used == 0
    @test length(h.arrays) == 2
    @test _new!(h) ≡ A1
    @test _new!(h) ≡ A2
    @test h.used == 2
    A3 = _new!(h)
    @test A3 isa Array{T,3} && size(A3) == D && A1 ≢ A3 ≢ A2
    @test length(h.arrays) == h.used == 3
end

####
#### tests for exposed API
####

@testset "ArrayPool" begin
    pool = ArrayPool()
    D = (3, 2)
    A1 = @inferred new!(pool, Array{Float64}, D)
    @test A1 isa Matrix{Float64} && size(A1) == D
    A2 = @inferred new!(pool, Vector{Int}, 9)
    @test A2 isa Vector{Int} && length(A2) == 9
    @test recycle!(pool) ≡ nothing
    @test A1 ≡ @inferred new!(pool, Matrix{Float64}, D)
    @test A2 ≡ @inferred new!(pool, Array{Int}, 9)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct CloudError <: SwaggerModel
    error::Any # spec type: Union{ Nothing, CloudErrorBody } # spec name: error

    function CloudError(;error=nothing)
        o = new()
        validate_property(CloudError, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type CloudError

const _property_map_CloudError = Dict{Symbol,Symbol}(Symbol("error")=>Symbol("error"))
const _property_types_CloudError = Dict{Symbol,String}(Symbol("error")=>"CloudErrorBody")
Base.propertynames(::Type{ CloudError }) = collect(keys(_property_map_CloudError))
Swagger.property_type(::Type{ CloudError }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CloudError[name]))}
Swagger.field_name(::Type{ CloudError }, property_name::Symbol) =  _property_map_CloudError[property_name]

function check_required(o::CloudError)
    true
end

function validate_property(::Type{ CloudError }, name::Symbol, val)
end


text:
###############################
#NOTE: Experiments with parsing
##############################
#=function parsestring(T::Type, s::AbstractString)
  v::Union{Nothing,T} = tryparse(T,s)

  return v
end

parsestring(Int, "3")
parsestring(Float64, "3.4")

function parsestringFloat64(s::AbstractString)
  v::Union{Nothing,Float64} = tryparse(Float64,s)

  return v
end

function parsestringInt(s::AbstractString)
  v::Union{Nothing,Int} = tryparse(Int,s)

  return v
end

parsestringInt("3")
parsestringFloat64("3.4")

function parsestring(s::AbstractString, ::Val{T}) where T
  v::Union{Nothing,T} = tryparse(T,s)

  return v
end

parsestring("3", Val(Int))
parsestring("3.4", Val(Float64))=#

###############################
#NOTE: Experiments with conversion
##############################
#=using Revise
using DataFrames, CodecLz4, Serialization

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

function lz4mwe(infile::String, outfile::String)
  @assert isfile(infile)
  df = deserialize(infile)
  println("size after read: $(size(df))")
  outdata = open(outfile, "w")

  s = LZ4CompressorStream(open(outfile, "w"))

  serialize(s, df)
  close(s)

  s = LZ4DecompressorStream(open(outfile))

  df = deserialize(s)

  close(s)
  println("size after decompression: $(size(df))")

end=#

###############################
#NOTE: Experiments with testing values in a df
###############################
#=
using DataStructures, DataFrames, Finometrics, BenchmarkTools

function testDS(NSet::Int = 1000, NSearched::Int = 1_000_000)
  df::DataFrame = DataFrame(cint = Vector{MInt}(rand(1:10^5, NSearched)),
    cfloat = Vector{MFloat64}(rand(1:10^5, NSearched)))

  local found::Int

  basefvec::Vector{Float64} = collect(1:NSet)
  baseivec::Vector{Int} = Vector{Int}(basefvec)

  #sftuple::NTuple{NSet, Float64} = NTuple{NSet, Float64}(basefvec)
  #situple::NTuple{NSet, Int} = NTuple{NSet, Int}(baseivec)

  sfset::Set{Float64} = Set{Float64}(basefvec)
  siset::Set{Int} = Set{Int}(baseivec)

  sfsortedset::SortedSet{Float64} = SortedSet{Float64}(basefvec)
  sisortedset::SortedSet{Int} = SortedSet{Int}(baseivec)

  print("\n")
  println("Testing Set")
  found=0
  @time found += sum((f::MFloat64->f ∈ sfset).(df.cfloat))
  @time found += sum((f::MInt->f ∈ siset).(df.cint))
  println("Set test complete found $found")

  println("Testing SortedSet")
  found=0
  @time found += sum((f::MFloat64->f ∈ sfsortedset).(df.cfloat))
  @time found += sum((f::MInt->f ∈ sisortedset).(df.cint))
  println("Set test complete found $found")

  println("Testing base")
  found=0
  @time found += sum((f::MFloat64->f ∈ basefvec).(df.cfloat))
  @time found += sum((f::MInt->f ∈ baseivec).(df.cint))
  println("Base test complete found $found")
  return nothing
end

testDS()=#

#=function mwe()
  println("first run:")
  for i ∈ 1:100:401
    basevec::Vector{Float64} = collect(1:i)
    print(i, ": ");
    @time NTuple{i, Float64}(basevec)
  end

  println("second run:")
  for i ∈ 1:100:401
    basevec::Vector{Float64} = collect(1:i)
    print(i, ": ");
    @time NTuple{i, Float64}(basevec)
  end

  println("uniontype first run")
  for i ∈ 1:18
    basevec::Vector{Union{Nothing,Float64}} = collect(1:i)
    print(i, ": ");
    @time NTuple{i, Union{Nothing,Float64}}(basevec)
  end

  println("uniontype second run")
  for i ∈ 1:18
    basevec::Vector{Union{Nothing,Float64}} = collect(1:i)
    print(i, ": ");
    @time NTuple{i, Union{Nothing,Float64}}(basevec)
  end
end
mwe()=#

#lz4mwe("data\\CRSP-D.jls", "data\\CRSP-D.csv.lz4")

##############
#NOTE: functions related to multi-threaded dataframes
################
#=using DataFrames
#WARNING: DO NOT RUN THIS

function tsmwecorrupt(N=100_000)
  df = DataFrame(rand(N,100))
  df.grpcol = (i->i%50).(1:N)

  Threads.@threads for sdf ∈ groupby(df, :grpcol)
    sdf.x3 .= -1.
  end

  println(sum(df.x3))
end

tsmwecorrupt()

function tsmwe(N=100_000)
  df = DataFrame(rand(N,100))
  df.grpcol = (i->i%50).(1:N)

  Threads.@threads for r ∈ eachrow(df)
    r.x3 = -1.
  end

  println(sum(df.x3))
end

tsmwe()

function tsmwe2(N=100_000)
  df = DataFrame(rand(N,100))
  df.grpcol = (i->i%50).(1:N)

  Threads.@threads for sdf ∈ collect(groupby(df, :grpcol))
    sdf.x3 .= -1.
  end

  println(sum(df.x3))
end

tsmwe2()=#

######################
#NOTE: Macro experiments
###################
using Random

macro tt(args...)
    println(args)

    na = length(args)
    if na != 1
        throw(ArgumentError("wrong number of arguments in @threads"))
    end
    ex = args[1]
    if !isa(ex, Expr)
        throw(ArgumentError("need an expression argument to @threads"))
    end
    if ex.head === :for
        return "ok"
    else
        throw(ArgumentError("unrecognized argument to @threads"))
    end
end

PAR = false
Random.seed!(11)
condg = true
macro mpar(cond, expr)
  quote
    if $(esc(cond))
        :($(Threads.@threads($expr)))
    else
        :($($expr))
    end
  end

end

function testmacro()

  v::Vector{Float64} = collect(1:100_000_000)
  cond::Bool = false

  @time @mpar true for i ∈ 1:length(v)
    v[i] = log10(i)
    end
  #=@mpar for i ∈ 1:length(v)
    v[i] = 1
end=#

  #println(s)
  println(sum(v))
end

testmacro()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module GeometricExamples

    using GeometricIntegrators

    export generate_html, package_version

    include("generate_html.jl")
    include("package_version.jl")

    export get_tableau_list_erk, get_tableau_list_firk, get_tableau_list_vprk,
           get_tableau_list_vprk_projection, get_tableau_list_vprk_radau,
           get_tableau_list_splitting

    include("tableau_lists.jl")

end


text:
using Algorand
using Test

@testset "Algorand.jl" begin
    # Write your tests here.
end


text:
import CLIMA

# Properties of dry air
CLIMA.Parameters.Planet.molmass_dryair(ps::EarthParameterSet) = 28.97e-3
CLIMA.Parameters.Planet.R_d(ps::EarthParameterSet) =
    gas_constant() / molmass_dryair(ps)
CLIMA.Parameters.Planet.kappa_d(ps::EarthParameterSet) = 2 / 7
CLIMA.Parameters.Planet.cp_d(ps::EarthParameterSet) = R_d(ps) / kappa_d(ps)
CLIMA.Parameters.Planet.cv_d(ps::EarthParameterSet) = cp_d(ps) - R_d(ps)

# Properties of water
CLIMA.Parameters.Planet.ρ_cloud_liq(ps::EarthParameterSet) = 1e3
CLIMA.Parameters.Planet.ρ_cloud_ice(ps::EarthParameterSet) = 916.7
CLIMA.Parameters.Planet.molmass_water(ps::EarthParameterSet) = 18.01528e-3
CLIMA.Parameters.Planet.molmass_ratio(ps::EarthParameterSet) =
    molmass_dryair(ps) / molmass_water(ps)
CLIMA.Parameters.Planet.R_v(ps::EarthParameterSet) =
    gas_constant() / molmass_water(ps)
CLIMA.Parameters.Planet.cp_v(ps::EarthParameterSet) = 1859
CLIMA.Parameters.Planet.cp_l(ps::EarthParameterSet) = 4181
CLIMA.Parameters.Planet.cp_i(ps::EarthParameterSet) = 2100
CLIMA.Parameters.Planet.cv_v(ps::EarthParameterSet) = cp_v(ps) - R_v(ps)
CLIMA.Parameters.Planet.cv_l(ps::EarthParameterSet) = cp_l(ps)
CLIMA.Parameters.Planet.cv_i(ps::EarthParameterSet) = cp_i(ps)
CLIMA.Parameters.Planet.T_freeze(ps::EarthParameterSet) = 273.15
CLIMA.Parameters.Planet.T_min(ps::EarthParameterSet) = 150.0
CLIMA.Parameters.Planet.T_max(ps::EarthParameterSet) = 1000.0
CLIMA.Parameters.Planet.T_icenuc(ps::EarthParameterSet) = 233.00
CLIMA.Parameters.Planet.T_triple(ps::EarthParameterSet) = 273.16
CLIMA.Parameters.Planet.T_0(ps::EarthParameterSet) = T_triple(ps)
CLIMA.Parameters.Planet.LH_v0(ps::EarthParameterSet) = 2.5008e6
CLIMA.Parameters.Planet.LH_s0(ps::EarthParameterSet) = 2.8344e6
CLIMA.Parameters.Planet.LH_f0(ps::EarthParameterSet) = LH_s0(ps) - LH_v0(ps)
CLIMA.Parameters.Planet.e_int_v0(ps::EarthParameterSet) =
    LH_v0(ps) - R_v(ps) * T_0(ps)
CLIMA.Parameters.Planet.e_int_i0(ps::EarthParameterSet) = LH_f0(ps)
CLIMA.Parameters.Planet.press_triple(ps::EarthParameterSet) = 611.657

# Properties of sea water
CLIMA.Parameters.Planet.ρ_ocean(ps::EarthParameterSet) = 1.035e3
CLIMA.Parameters.Planet.cp_ocean(ps::EarthParameterSet) = 3989.25

# Planetary parameters
CLIMA.Parameters.Planet.planet_radius(ps::EarthParameterSet) = 6.371e6
CLIMA.Parameters.Planet.day(ps::EarthParameterSet) = 86400
CLIMA.Parameters.Planet.Omega(ps::EarthParameterSet) = 7.2921159e-5
CLIMA.Parameters.Planet.grav(ps::EarthParameterSet) = 9.81
CLIMA.Parameters.Planet.year_anom(ps::EarthParameterSet) = 365.26 * day(ps)
CLIMA.Parameters.Planet.orbit_semimaj(ps::EarthParameterSet) = 1 * astro_unit()
CLIMA.Parameters.Planet.TSI(ps::EarthParameterSet) = 1362
CLIMA.Parameters.Planet.MSLP(ps::EarthParameterSet) = 1.01325e5

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
s = """
Dict(:a=>1, :b=>2, :a=>3)
"""
msgs = lintstr(s)
@test msgs[1].code == :E334
@test contains(msgs[1].message, "duplicate key in Dict")

s = """
Dict{Symbol,Int}(:a=>1, :b=>"")
"""
msgs = lintstr(s)
@test msgs[1].code == :E532
@test contains(msgs[1].message, "multiple value types detected. Use Dict{K,Any}() for " *
    "mixed type dict")

s = """
Dict{Symbol,Int}(:a=>1, "b"=>2)
"""
msgs = lintstr(s)
@test msgs[1].code == :E531
@test contains(msgs[1].message, "multiple key types detected. Use Dict{Any,V}() for " *
    "mixed type dict")


text:
#=function FR.filter_exp!(u::AbstractMatrix{T}, args...) where {T<:AbstractFloat}
    nx, nz = size(u)
    λx, λξ = args[1:2]

    if length(args) >= 3
        Ncx = args[3]
    else
        Ncx = 0
    end
    if length(args) >= 4
        Ncξ = args[4]
    else
        Ncξ = 0
    end

    σ1 = FR.filter_exp1d(nx-1, λx, Ncx)
    σ2 = FR.filter_exp1d(nz-1, λξ, Ncξ)

    for j in axes(u, 2), i in axes(u, 1)
        u[i, j] *= σ1[i] * σ2[j]
    end

    return nothing
end=#

function KitBase.filter_exp!(u::AbstractMatrix{T}, args...) where {T<:AbstractFloat}
    nx, nz = size(u)
    spx = args[1]

    if length(args) >= 2
        spy = args[2]
    else
        spy = spx
    end
    if length(args) >= 3
	lamdt = args[3]
    else
	lamdt = 1
    end

    σ = filter_exp2d(nx-1, nz-1, spx,spy).^lamdt

    for j in axes(u, 2), i in axes(u, 1)
        u[i, j] *= σ[i, j]
    end

    return nothing
end

function filter_exp2d(Nx, Ny, spx, spy, Nc = 0)
    alpha = -log(eps())

    filterdiag = ones((Nx + 1), (Ny + 1))
    for i = 0:Nx
        for j = 0:Ny
            if i + j >= Nc
     		#filterdiag[i+1, j+1] = exp(-alpha * ((i + j - Nc) / (Nx + Ny - Nc))^sp)
		filterdiag[i+1, j+1] = exp(-alpha * ((i/(Nx+1))^spx + (j/(Ny+1))^spy))
            end
        end
    end

    return filterdiag
end

function KitBase.filter_houli!(u::AbstractMatrix{T}, args...) where {T<:AbstractFloat}
    nx, nz = size(u)
    λx, λξ = args[1:2]

    if length(args) >= 3
        Ncx = args[3]
    else
        Ncx = 0
    end
    if length(args) >= 4
        Ncξ = args[4]
    else
        Ncξ = 0
    end

    σ1 = FR.filter_exp1d(nx-1, λx, Ncx)
    σ2 = FR.filter_exp1d(nz-1, λξ, Ncξ)

    for i in eachindex(σ1)
        if i/length(σ1) <= 2/3
            σ1[i] = 1.0
        end
    end
    for i in eachindex(σ2)
        if i/length(σ2) <= 2/3
            σ2[i] = 1.0
        end
    end

    for j in axes(u, 2), i in axes(u, 1)
        u[i, j] *= σ1[i] * σ2[j]
    end

    return nothing
end

begin # compute L1 norms of basis
    Nx = ps.nx
    nLocal = ps.deg + 1
    nMoments = uq.nm + 1

    NxHat = 100
    xHat = collect(range(-1,stop = 1,length = NxHat))
    dxHat = xHat[2]-xHat[1];
    PhiL1 = zeros(nLocal,nMoments)
    for i = 1:nLocal
        for j = 1:nMoments
            PhiJ = sum(@. abs(uq.op.quad.weights * uq.phiRan[:, j])) / (uq.t2Product[j-1, j-1] + 1.e-7)
            PhiI = dxHat*sum(abs.(JacobiP(xHat, 0, 0, i-1)))
            PhiL1[i,j] = PhiI*PhiJ
        end
    end
end

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
using MiniLoggers
using Documenter

DocMeta.setdocmeta!(MiniLoggers, :DocTestSetup, :(using MiniLoggers); recursive=true)

makedocs(;
    modules=[MiniLoggers],
    authors="Andrey Oskin",
    repo="https://github.com/JuliaLogging/MiniLoggers.jl/blob/{commit}{path}#{line}",
    sitename="MiniLoggers.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaLogging.github.io/MiniLoggers.jl",
        siteurl="https://github.com/JuliaLogging/MiniLoggers.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaLogging/MiniLoggers.jl",
)


text:
"""
An OPF formulation conforming to the ARPA-e GOC Challenge 1 specification.
Power balance and branch flow constraints are strictly enforced.
The primary departure from the PowerModels standard formulation is dispatchable
bus shunts and a slight change in the transformer model.
"""
function run_opf_shunt(file, model_constructor, solver; kwargs...)
    return _PM.run_model(file, model_constructor, solver, build_opf_shunt; ref_extensions=[ref_add_goc!], kwargs...)
end

function build_opf_shunt(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm)
    _PM.variable_gen_power(pm)
    _PM.variable_branch_power(pm)

    variable_shunt_admittance_imaginary(pm)

    _PM.objective_min_fuel_cost(pm)

    _PM.constraint_model_voltage(pm)

    for i in ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        constraint_power_balance_shunt_dispatch(pm, i)
    end

    for (i,branch) in ref(pm, :branch)
        constraint_ohms_yt_from_goc(pm, i)
        _PM.constraint_ohms_yt_to(pm, i)

        _PM.constraint_voltage_angle_difference(pm, i)

        _PM.constraint_thermal_limit_from(pm, i)
        _PM.constraint_thermal_limit_to(pm, i)
    end
end

"""
An OPF formulation conforming to the ARPA-e GOC Challenge 1 specification.
Power balance are strictly enforced and the branch flow violations are
penalized based on a conservative linear approximation of the formulation's
flow violation penalty specification.
"""
function run_opf_cheap(file, model_constructor, solver; kwargs...)
    return _PM.run_model(file, model_constructor, solver, build_opf_cheap; ref_extensions=[ref_add_goc!], kwargs...)
end

function build_opf_cheap(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm)
    _PM.variable_gen_power(pm)
    _PM.variable_branch_power(pm, bounded=false)

    variable_branch_power_slack(pm)
    variable_shunt_admittance_imaginary(pm)

    _PM.constraint_model_voltage(pm)

    for i in ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        constraint_power_balance_shunt_dispatch(pm, i)
    end

    for (i,branch) in ref(pm, :branch)
        constraint_ohms_yt_from_goc(pm, i)
        _PM.constraint_ohms_yt_to(pm, i)

        _PM.constraint_voltage_angle_difference(pm, i)

        constraint_thermal_limit_from_soft(pm, i)
        constraint_thermal_limit_to_soft(pm, i)
    end

    ##### Setup Objective #####
    _PM.objective_variable_pg_cost(pm)
    # explicit network id needed because of conductor-less
    pg_cost = var(pm, :pg_cost)
    sm_slack = var(pm, :sm_slack)

    @objective(pm.model, Min,
        sum( pg_cost[i] for (i,gen) in ref(pm, :gen) ) +
        sum( 5e5*sm_slack[i] for (i,branch) in ref(pm, :branch_sm_active) )
    )
end

"""
A variant of `run_opf_cheap` model, specialized for solving very large
AC Power Flow models in rectangular coordinates for faster derivative
computations.  Support sparse collections of flow constrains for
increased performance.
"""
function run_opf_cheap_lazy_acr(file, solver; kwargs...)
    return _PM.run_model(file, _PM.ACRPowerModel, solver, build_opf_cheap_lazy_acr; ref_extensions=[ref_add_goc!], kwargs...)
end

""
function build_opf_cheap_lazy_acr(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm, bounded=false)
    _PM.variable_gen_power(pm)

    variable_branch_power_slack(pm)
    variable_shunt_admittance_imaginary(pm)

    variable_bus_voltage_magnitude_delta(pm)
    variable_gen_power_real_delta(pm)

    vvm = var(pm)[:vvm] = @variable(pm.model,
        [i in ids(pm, :bus)], base_name="vvm",
        lower_bound = ref(pm, :bus, i, "vmin")^2,
        upper_bound = ref(pm, :bus, i, "vmax")^2,
        start = 1.0
    )

    _PM.constraint_model_voltage(pm)

    for i in ids(pm, :branch)
        expression_branch_power_ohms_yt_from_goc(pm, i)
        _PM.expression_branch_power_ohms_yt_to(pm, i)
    end

    vr = var(pm, :vr)
    vi = var(pm, :vi)
    for (i,bus) in ref(pm, :bus)
        vm_midpoint = (bus["vmax"] + bus["vmin"])/2.0
        vm_target = min(vm_midpoint + 0.04, bus["vmax"])
        #vm_target = bus["vm"]

        @constraint(pm.model, vr[i]^2 + vi[i]^2 == vvm[i])
        @constraint(pm.model, vvm[i] == vm_target^2 + var(pm, :vvm_delta, i))
    end

    for i in ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end
    #Memento.info(_LOGGER, "misc constraints time: $(time() - start_time)")

    start_time = time()
    for (i,branch) in ref(pm, :branch)
        if haskey(branch, "rate_a")
            f_bus_id = branch["f_bus"]
            t_bus_id = branch["t_bus"]
            f_idx = (i, f_bus_id, t_bus_id)
            t_idx = (i, t_bus_id, f_bus_id)

            p_fr = @variable(pm.model, base_name="p_fr", start = 0.0)
            q_fr = @variable(pm.model, base_name="q_fr", start = 0.0)
            p_to = @variable(pm.model, base_name="p_to", start = 0.0)
            q_to = @variable(pm.model, base_name="q_to", start = 0.0)

            @constraint(pm.model, var(pm, :p, f_idx) == p_fr)
            @constraint(pm.model, var(pm, :q, f_idx) == q_fr)
            @constraint(pm.model, var(pm, :p, t_idx) == p_to)
            @constraint(pm.model, var(pm, :q, t_idx) == q_to)

            rating = branch["rate_a"]
            sm_slack = var(pm, :sm_slack, i)
            JuMP.@constraint(pm.model, p_fr^2 + q_fr^2 <= (rating + sm_slack)^2)
            JuMP.@constraint(pm.model, p_to^2 + q_to^2 <= (rating + sm_slack)^2)
        end
    end
    #Memento.info(_LOGGER, "flow expr time: $(time() - start_time)")

    start_time = time()
    for (i,gen) in ref(pm, :gen)
        constraint_gen_power_real_deviation(pm, i)
    end
    #Memento.info(_LOGGER, "gen expr time: $(time() - start_time)")

    p = var(pm, :p)
    q = var(pm, :q)
    pg = var(pm, :pg)
    qg = var(pm, :qg)
    bs = var(pm, :bs)
    for (i,bus) in ref(pm, :bus)
        #_PM.constraint_power_balance(pm, i)

        bus_arcs = ref(pm, :bus_arcs, i)
        bus_gens = ref(pm, :bus_gens, i)
        bus_loads = ref(pm, :bus_loads, i)
        bus_shunts_const = ref(pm, :bus_shunts_const, i)
        bus_shunts_var = ref(pm, :bus_shunts_var, i)

        bus_pd = Dict(k => ref(pm, :load, k, "pd") for k in bus_loads)
        bus_qd = Dict(k => ref(pm, :load, k, "qd") for k in bus_loads)

        bus_gs_const = Dict(k => ref(pm, :shunt, k, "gs") for k in bus_shunts_const)
        bus_bs_const = Dict(k => ref(pm, :shunt, k, "bs") for k in bus_shunts_const)

        @constraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*vvm[i])
        @constraint(pm.model, 0 == - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*vvm[i] + sum(bs[s]*vvm[i] for s in bus_shunts_var))
    end
    #Memento.info(_LOGGER, "power balance constraint time: $(time() - start_time)")

    vvm_delta = var(pm, :vvm_delta)
    sm_slack = var(pm, :sm_slack)
    pg_delta = var(pm, :pg_delta)

    @objective(pm.model, Min,
        sum( 1e7*vvm_delta[i]^2 for (i,bus) in ref(pm, :bus)) +
        sum( 5e5*sm_slack[i] for (i,branch) in ref(pm, :branch_sm_active)) +
        sum( 1e5*pg_delta[i]^2 for (i,gen) in ref(pm, :gen))
    )
end

""
function run_opf_cheap_target_acp(file, solver; kwargs...)
    return _PM.run_model(file, _PM.ACPPowerModel, solver, build_opf_cheap_target_acp; ref_extensions=[ref_add_goc!], kwargs...)
end

""
function build_opf_cheap_target_acp(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm)
    _PM.variable_gen_power(pm)
    _PM.variable_branch_power(pm, bounded=false)

    variable_branch_power_slack(pm)
    variable_shunt_admittance_imaginary(pm)

    variable_bus_voltage_magnitude_delta(pm)
    variable_gen_power_real_delta(pm)

    _PM.constraint_model_voltage(pm)

    vm = var(pm, :vm)
    for (i,bus) in ref(pm, :bus)
        vm_midpoint = (bus["vmax"] + bus["vmin"])/2.0
        vm_target = min(vm_midpoint + 0.04, bus["vmax"])

        @constraint(pm.model, vm[i] == vm_target + var(pm, :vvm_delta, i))
    end

    for i in ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    start_time = time()
    for (i,gen) in ref(pm, :gen)
        constraint_gen_power_real_deviation(pm, i)
    end

    for i in ids(pm, :bus)
        constraint_power_balance_shunt_dispatch(pm, i)
    end

    for (i,branch) in ref(pm, :branch)
        constraint_ohms_yt_from_goc(pm, i)
        _PM.constraint_ohms_yt_to(pm, i)

        _PM.constraint_voltage_angle_difference(pm, i)

        _PM.constraint_thermal_limit_from(pm, i)
        _PM.constraint_thermal_limit_to(pm, i)
    end

    vvm_delta = var(pm, :vvm_delta)
    sm_slack = var(pm, :sm_slack)
    pg_delta = var(pm, :pg_delta)

    @objective(pm.model, Min,
        sum( 1e8*vvm_delta[i]^2 for (i,bus) in ref(pm, :bus)) +
        sum( 5e5*sm_slack[i] for (i,branch) in ref(pm, :branch)) +
        sum( 1e5*pg_delta[i]^2 for (i,gen) in ref(pm, :gen))
    )
end


text:
module ResNet
export resnet, ResNetModel

using Flux

include("model.jl")

end


text:
using AbstractMCMC
using Atom.Progress: JunoProgressLogger
using ConsoleProgressMonitor: ProgressLogger
using IJulia
using LoggingExtras: TeeLogger, EarlyFilteredLogger
using TerminalLoggers: TerminalLogger
using Transducers

using Distributed
import Logging
using Random
using Statistics
using Test
using Test: collect_test_logs

const LOGGERS = Set()
const CURRENT_LOGGER = Logging.current_logger()

include("utils.jl")

@testset "AbstractMCMC" begin
    include("sample.jl")
    include("stepper.jl")
    include("transducer.jl")
    include("deprecations.jl")
end


text:
path = "docs/src/assets/"

using ISA, Plots
fs = 16_000
t = 0.0:1/fs:3.0
𝐶₀ = AMFMtriplet(t->0.2+0.8cos(11t), t->200.0, 0.0)
𝐶₁ = AMFMtriplet(t->exp(-abs(t/3)), t->100t, 0.1)
𝐶₂ = AMFMtriplet(t->𝒩ᵤ(t; μ=1.5, σ=1.0), t->150 + 125sin(5t), π)
𝐶₃ = AMFMtriplet(t->u(t-1.725)-u(t-2.475), t->50, 0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂,𝐶₃])
z = AMFMmodel(𝑆)

plot(𝐶₀,timeaxis=t)
plot!(𝐶₁,timeaxis=t)
plot!(𝐶₂,timeaxis=t)
plot!(𝐶₃,linealpha = sign.(abs.(𝐶₃.a.(t))),timeaxis = t)
png(path*"paperFigureIS.png")

plot(z,timeaxis=t)
png(path*"paperFigureArgand.png")


text:
# works well. no complaints
function solve_hpd(domain_name, problem_name; max_levels=20)
    graph = create_funnel_graph(domain_name, problem_name; max_levels=max_levels)

    model = Model(Gurobi.Optimizer)
    set_silent(model)

    T = graph.num_levels
    max_nf = max([length([prop for prop in graph.props[i][:discrete]]) for i=1:T]...)
    max_na = max([length([act for act in graph.acts[i] ]) for i=1:T-1]...)

    @variable(model, p[i=1:max_nf, j=1:T], Bin)
    @variable(model, a[i=1:max_na, j=1:T-1], Bin)
    @variable(model, xᵣ[i=1:2, j=1:T-1])

    # init constraints
    num_init_props = length(graph.props[1][:discrete])
    [@constraint(model, p[i, 1] == 1) for i = 1: num_init_props]
    [@constraint(model, p[i, 1] == 0) for i = num_init_props+1 : max_nf]
    @constraint(model, xᵣ[1, 1] == graph.initprops[:continuous][1][:xr])
    @constraint(model, xᵣ[2, 1] == graph.initprops[:continuous][1][:yr])

    # goal constraints
    goalinds = get_fluent_terminal_indices(graph; loc=:goals)
    [@constraint(model, p[f, T] == 1) for f in goalinds]
    @constraint(model, xᵣ[1, T-1] == graph.goalprops[:continuous][:xr])
    @constraint(model, xᵣ[2, T-1] == graph.goalprops[:continuous][:yr])

    # precondition constraints
    for t = 1:T-1
        fluents = graph.props[t][:discrete]
        for f = 1:length(fluents)
            fluent = fluents[f]
            act_inds = get_actions_precondition_sat(fluent, t, graph)
            [@constraint(model, a[i,t] ≤ p[f,t]) for i in act_inds]
        end
    end

    # effect constraints
    for t = 1:T-1
        fluents = graph.props[t+1][:discrete]
        for f = 1:length(fluents)
            fluent = fluents[f]
            act_inds = get_actions_effect_sat(fluent, t, graph)
            fnext = get_fluent_index(fluent, t+1, graph)
            if !isnothing(fnext) && !isempty(act_inds)
                @constraint(model, p[fnext, t+1] ≤ sum([a[i,t] for i in act_inds]))
            end
        end
        na = length(graph.acts[t])
        @constraint(model, a[na+1:max_na, t] .== 0)
        # for i = 1:na
        #     if graph.acts[t][i].name == :noop  && i<=max_na
        #         @constraint(model, a[i, t] == 0)
        #     end
        # end
        # @constraint(model, sum(a[:, t])==1)
    end

    # mutex constraints
    for t = 1:T-1
        for f=1:length(graph.props[t][:discrete])
            fluent = graph.props[t][:discrete][f]
            act_pairs = get_mutex_actions(fluent, t, graph)
            [@constraint(model, a[i[1], t] + a[i[2], t] <= 1) for i in act_pairs]
        end
    end

    # continuous constraints
    # does this for all actions. There might be some contradictions.
    # I only want to do it for a single action in every level
    for t = 1:T-1
        for i=1:length(graph.acts[t])
            if graph.acts[t][i].name != :noop
                if isa(graph.acts[t][i].end_region[:xr], Float64)
                    @constraint(model, a[i,t] => {xᵣ[1,t] == graph.acts[t][i].end_region[:xr]})
                    @constraint(model, a[i,t] => {xᵣ[2,t] == graph.acts[t][i].end_region[:yr]})
                else
                    @constraint(model, a[i,t] => {graph.acts[t][i].end_region[:xr][1] ≤ xᵣ[1,t] ≤ graph.acts[t][i].end_region[:xr][2]})
                    @constraint(model, a[i,t] => {graph.acts[t][i].end_region[:yr][1] ≤ xᵣ[2,t] ≤ graph.acts[t][i].end_region[:yr][2] })
                end
            end
        end
    end

    # t=@variable(model)
    # @constraint(model, [t; sum([xᵣ[:,t]-xᵣ[:,t+1] for t=1:T-2])] in SecondOrderCone())
    # @objective(model, Min, t)

    # @NLexpression(model, n, [nm(xᵣ[:,t]-xᵣ[:,t+1]) for t=1:T-1])
    # @NLobjective(model, Min, sum(n))

    # @objective(model, Min, sum(a))

    nm(xs) = sum([x^2 for x in xs])
    @objective(model, Min, sum([nm(xᵣ[:,t]-xᵣ[:,t+1]) for t=1:T-2]))

    optimize!(model)
    optimal_funnel_indices = value.(a)
    plan = render_action(optimal_funnel_indices, graph)
    X = value.(xᵣ)
    robot_poses = [X[:,i] for i in 1:size(X,2)]
    sol = Soln(plan, robot_poses)
    return sol
end

function get_fluent_terminal_indices(graph::Graph; loc=:init)
    indices = Set()
    T = graph.num_levels
    if loc == :init
        for i=1:length(graph.initprops[:discrete])
            ind = findall(x->x == graph.initprops[:discrete], graph.props[1][:discrete])[1]
            push!(indices, ind)
        end
    else
        for i=1:length(graph.goalprops[:discrete])
            ind = findall(x->x == graph.goalprops[:discrete][i], graph.props[T][:discrete])[1]
            push!(indices, ind)
        end
    end
    return collect(indices)
end

function get_actions_precondition_sat(fluent, level::Int, graph::Graph)
    actions = [act for act in graph.acts[level] ]
    act_inds = []
    for i = 1:length(actions)
        act = graph.acts[level][i]
        if fluent in act.pos_prec  push!(act_inds, i) end
    end
    return act_inds
end

function get_actions_effect_sat(fluent, level::Int, graph::Graph)
    actions = [act for act in graph.acts[level] ]
    act_inds = []
    for i = 1:length(actions)
        act = graph.acts[level][i]
        if fluent in act.pos_eff push!(act_inds, i) end
    end
    return act_inds
end

function get_fluent_index(fluent, level::Int, graph::Graph)
    fluents = graph.props[level][:discrete]
    for i = 1:length(fluents)
        if fluents[i] == fluent return i end
    end
    return nothing
end

function get_mutex_actions(fluent, level, graph)
    pref = get_actions_precondition_sat(fluent, level, graph)
    preneff = union(pref, get_actions_effect_sat(fluent, level, graph))
    act_mutexes = []
    [push!(act_mutexes, [i,j]) for i in pref for j in preneff if i!=j && !([j,i] in act_mutexes)]

    return act_mutexes
end

function get_mutex_indices(level, graph)
    act_mutexes = []
    for pair in graph.μacts[level]
        i = findall(x->x == pair[1], graph.acts[level])[1]
        j = findall(x->x == pair[2], graph.acts[level])[1]

        if !([j,i] in act_mutexes) push!(act_mutexes, [i,j]) end
    end
    return act_mutexes
end

function render_action(y, graph::Graph)
    T = graph.num_levels
    actions = []
    for t=1:T-1
        ind = findall(x->x == 1.0, y[1:end,t])
        id = ind[1]
        # for i in ind if graph.acts[t][i].name != :noop id = i; break; end end
        act  = graph.acts[t][id]
        push!(actions, (act.name, act.params))
    end
    return actions

end


text:
@testset "errors" begin

for (text, msg) in [("a 1 ]", "Expected key"),
                    ("a [1 2]", "Expected ]"),
                    ("a [a -w]", "Expected value")]
    try
        println(parse_raw(text))
    catch x
        if isa(x, ParserError)
            @test occursin(msg, x.msg)
        else
            println(x)
        end
    end

    @test_throws ParserError{Int} parse_raw(text)
end

s = open(s -> read(s, String), "gml/error.gml")
try
    parse_raw(s)
    @test false
catch x
    @test isa(x, ParserError)
    @test x.msg == "Expected ] at (2,15)\n  node [ id 1 \"sausage\" ]\n              ^\n"
end

parse_dict("graph [a 0] graph [a 1 a 2]"; lists=[:a,:graph])
for unsafe in (true, false)
    try
        parse_dict("graph [a 0] graph [a 1 a 2]"; unsafe=unsafe)
        @test unsafe == true  # can parse when unsafe
    catch x
        println(x)
        @test unsafe == false # error when safe
        @test isa(x, GMLError)
        @test occursin(x.msg, "a is a list")
    end
end

end


text:
#
# units.jl -
#
# Define constants for physical units.  All physical units (only lengths and
# angles) are assumed to be in SI units (hence lengths are in meters, angles
# are in radians, etc.).
#

"""
Usage:

``julia
using FourierOptics.Units
```

provides short names to be used as units.  For instance `λ = 500nm` or `λ =
0.5µm` can be used to define the wavelength, or `f = 15.2mm` to define a focal
length.

"""
module Units

export
    showlength,
    showposition,
    m, cm, mm, µm, nm,
    rd, deg, arcdeg, arcmin, arcsec, marcsec, µarcsec, mas, µas

# Lengths are given in SI units (meters).
const m  = 1.0
const cm = 1E-2
const mm = 1E-3
const µm = 1E-6
const nm = 1E-9

# Angles are given in SI units (radians).
const rd = 1.0
const arcdeg = π/180
const arcmin = π/(180*60)
const arcsec = π/(180*60*60)
const marcsec = 1E-3arcsec
const µarcsec = 1E-6arcsec
const deg = arcdeg
const mas = marcsec
const µas = µarcsec

"""
```julia
showlength(io, x)
```

prints length `x` (assumed to be in SI units, that is meters) to `io`
in a human readable way.

"""
function showlength(io::IO, x::Real)
    a = abs(x)
    if 1mm ≤ a < 1m
        print(io, @sprintf("%.6gmm", x/1mm))
    elseif 1µm ≤ a < 1mm
        print(io, @sprintf("%.6gµm", x/1µm))
    elseif 1nm ≤ a < 1µm
        print(io, @sprintf("%.6gnm", x/1nm))
    else
        print(io, @sprintf("%.6gm", x/1m))
    end
end

showlength(io::IO, x::Union{Rational,Irrational}) =
    showlength(io, convert(Float64, x))

"""
```julia
showposition(io, (x,y))
```

prints position `(x,y)` (both coordinates are assumed to be in SI units, that
is meters) to `io` in a human readable way.

"""
function showposition(io::IO, c::NTuple{2,Union{Real,Rational,Irrational}})
    print(io, "(")
    showlength(io, c[1])
    print(io, ",")
    showlength(io, c[2])
    print(io, ")")
end

end


text:
module Wave

    include("network.jl")
    include("dataprep.jl")
    include("train.jl")

end


text:
# TODO: add preconditioner to the Problem description??

type WalkerProblem
   V
   E
   dE
   x0
   v0
   id
end

testsets = []
methods = []

maxnumdE = 1_000

# ============ Generate the Benchmark Tests ===========

# TODO: maybe we can assign a "good" initial step too all the methods?

# MullerPotential with good IC
V = MullerPotential()
x0, v0 = ic_dimer(V, :near)
E, dE = objective(V)
push!(testsets, WalkerProblem(V, E, dE, x0, v0, "Muller:near"))

# MüllerPotential with bad IC
V = MullerPotential()
x0, v0 = ic_dimer(V, :far)
E, dE = objective(V)
push!(testsets, WalkerProblem(V, E, dE, x0, v0, "Muller:far"))

# ============ List of Methods to Benchmark ===========

push!(methods,  BBDimer(a0_trans=0.002, a0_rot=0.002, maxnumdE=maxnumdE,
                        verbose=verbose, ls = StaticLineSearch(), id="BBDimer") )

push!(methods,  BBDimer(a0_trans=0.002, a0_rot=0.002, maxnumdE=maxnumdE,
                        verbose=verbose, ls = Backtracking(), id="BBDimer+LS") )

push!(methods,  RotOptimDimer(a_trans=0.001, maxnumdE=maxnumdE, verbose=verbose) )

push!(methods,  SuperlinearDimer(maximum_translation=0.1, max_num_rot=1, maxnumdE=maxnumdE, verbose=verbose) )

# ============= run the all the benchmarks ================

using DataFrames
results = DataFrame()
results[:Methods] = [m.id for m in methods]
for t in testsets
   r = Int[]
   for method in methods
      try
         x, v, log = run!(method, t.E, t.dE, t.x0, t.v0)
         push!(r, maximum(log.numdE))
         @show x, v
      catch e
         println("($(t.id), $(method.id)) : failed")
         @show e
         push!(r, -1)
      end
      # TODO: test the result is correct
   end
   results[Symbol(t.id)] = r
end

display(results)


text:
using TopOpt
# import Makie
# using TopOpt.TopOptProblems.Visualization: visualize

using TimerOutputs
using Suppressor

# function run_topopt()
println("Start running.")
# https://github.com/KristofferC/TimerOutputs.jl
to = TimerOutput()
reset_timer!(to)
Nonconvex.NonconvexCore.show_residuals[] = true

# Define the problem
E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0; # downward force

# Parameter settings
V = 0.3 # volume fraction
# xmin = 0.001 # minimum density
xmin = 1e-6 # minimum density
rmin = 2.0; # density filter radius

nels = (1440, 480) # (1000, 1000) # (500, 500) # (100, 100) # (1440, 480) # | (720, 240)
sizes = (1.0, 1.0)
@timeit to "problem def" problem = PointLoadCantilever(Val{:Linear}, nels, sizes, E, v, f);

# Define a finite element solver
@timeit to "penalty def" penalty = TopOpt.PowerPenalty(3.0)
@timeit to "solver def" solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty);

# Define compliance objective
@timeit to "objective def" begin
    # Define compliance objective
    comp = Compliance(problem, solver)
    filter = DensityFilter(solver; rmin=rmin)
    obj = x -> comp(filter(x))
end

# Define volume constraint
@timeit to "constraint def" begin
    volfrac = TopOpt.Volume(problem, solver)
    constr = x -> volfrac(filter(x)) - V
end

# Define subproblem optimizer
# ! seems to be absolute diff
mma_options =
    options = MMAOptions(;
        maxiter=1000, tol=Tolerance(; x=1e-3, fabs=1e-3, frel=1e-3, kkt=1e-3)
    )
convcriteria = GenericCriteria()
# convcriteria = KKTCriteria()

x0 = fill(V, length(solver.vars))
@timeit to "optimizer def" optimizer = Optimizer(
    obj, constr, x0, MMA87(); options=mma_options, convcriteria=convcriteria
);

# Define SIMP optimizer
@timeit to "simp def" simp = SIMP(optimizer, solver, penalty.p);

# Solve
# initial solution, critical to set it to volfrac! (blame non-convexity :)
@timeit to "simp run" result = simp(x0);

# Print the timings in the default way
println()
show(to)

@show result.convstate
@show result.objval
try
    @show optimizer.workspace.iter
catch
    # IpoptWorkspace has no field iter
end

output = @capture_out begin
    show(to)
    @show result.convstate
    @show result.objval
    try
        @show optimizer.workspace.iter
    catch
        # IpoptWorkspace has no field iter
    end
end;

open("jl-topopt.py_$(nels).txt", "w") do io
    write(io, output)
end

# # # Visualize the result using Makie.jl
# fig = visualize(problem; topology=result.topology,
#     default_exagg_scale=0.07, scale_range=10.0, vector_linewidth=3, vector_arrowsize=0.5)
# Makie.display(fig)

# Makie.save("jl-topopt.py_$(nels).png", fig)

# return problem, result
# end

# problem, result = run_topopt();


text:
# TODO: These tests should be expanded to cover everything

@testset "fixhorizon Grid World" begin
    gw = SimpleGridWorld()
    fhgw = fixhorizon(gw, 2)

    @test length(states(fhgw)) == 3*length(states(gw))
    @test collect(first(state) for state in stage_states(fhgw, 1)) == collect(states(gw))
    @test HorizonLength(fhgw) == FiniteHorizon()
    @test horizon(fhgw) == 2
    @test has_consistent_distributions(fhgw)
    is = (convert(statetype(gw), [1,1]), 1)::statetype(fhgw)
    @show actiontype(fhgw)
    @test length(collect(stepthrough(fhgw, FunctionPolicy(s->:left), is, "s,a,r,sp"))) == 2
end

@testset "fixhorizon Baby" begin
    m = BabyPOMDP()
    fhb = fixhorizon(m, 2)

    @test length(states(fhb)) == 3*length(states(m))
    @test collect(first(state) for state in stage_states(fhb, 1)) == collect(states(m))
    @test HorizonLength(fhb) == FiniteHorizon()
    @test horizon(fhb) == 2
    @test has_consistent_distributions(fhb)
    @test length(collect(stepthrough(fhb, FunctionPolicy(s->true)))) == 2
end

@testset "fixhorizon General" begin
    m = BabyPOMDP()
    fhb = fixhorizon(m, 3)
    state = (false, 2)
    action = true

    states = [(false, 2), (true, 2)]
    @test stateindex(fhb, ordered_states(fhb)[3]) == 3
    @test stage_stateindex(fhb, ordered_states(fhb)[3]) == 1
    @test ordered_states(fhb)[3] == state
    @test ordered_stage_states(fhb, 2) == [(false, 2), (true, 2)]

    obs = [(false, 2), (true, 2)]
    @test collect(first(obs) for obs in stage_observations(fhb, 1)) == collect(observations(m))
    @test obsindex(fhb, obs[2]) == 4
    @test stage_obsindex(fhb, obs[2]) == 2
    @test ordered_stage_observations(fhb, 2) == [(false, 2), (true, 2)]
end

@testset "solver" begin
    fhgw = fixhorizon(SimpleGridWorld(), 3)
    # TODO: Change test to pass without boolean value
    # Dependencies are destroying testing with this (probably until the repos are registered ?)
    # @test test_solver(FiniteHorizonSolver(), fhgw) == 9.025
end


text:
function MOI.get(opt::Optimizer, ::MOI.TerminationStatus)
    if opt.result === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
    stat = opt.result.termination_status
    if stat == OPTIMAL
        return MOI.OPTIMAL
    elseif stat == INFEASIBLE
        return MOI.INFEASIBLE
    elseif stat == INF_OR_UNBOUNDED
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif stat == EARLY_TERMINATION
        # TODO: Plug in a time limit, which will invalidate this.
        return MOI.NODE_LIMIT
    elseif stat == NOT_OPTIMIZED
        return MOI.OPTIMIZE_NOT_CALLED
    else
        @warn "Unhandled termination status $stat."
        return MOI.OTHER_ERROR
    end
end

function _has_feasible_solution(opt::Optimizer)
    if opt.result === nothing
        return false
    end
    x = opt.result.best_solution
    return !isempty(x)
end

function MOI.get(opt::Optimizer, ::MOI.PrimalStatus)
    if opt.result === nothing
        return MOI.NO_SOLUTION
    end
    if _has_feasible_solution(opt)
        stat = opt.result.termination_status
        @assert stat != INFEASIBLE && stat != INF_OR_UNBOUNDED
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function MOI.set(opt::Optimizer, ::MOI.ObjectiveValue, value::Real)
    @assert opt.result === nothing
    opt.primal_bound = value
    return nothing
end

function MOI.get(opt::Optimizer, ::MOI.ObjectiveValue)
    return if opt.result === nothing
        opt.primal_bound
    elseif _is_max_sense(opt)
        -opt.result.primal_bound
    else
        opt.result.primal_bound
    end
end

function MOI.get(opt::Optimizer, ::MOI.ObjectiveBound)
    @assert opt.result !== nothing
    return _is_max_sense(opt) ? -opt.result.dual_bound : opt.result.dual_bound
end

function MOI.get(opt::Optimizer, ::MOI.SolveTime)
    if opt.result === nothing
        return NaN
    end
    return opt.result.total_elapsed_time_sec
end

function MOI.get(opt::Optimizer, ::MOI.SimplexIterations)
    return opt.result.total_simplex_iters
end

function MOI.get(opt::Optimizer, ::MOI.NodeCount)
    return opt.result.total_node_count
end

function MOI.get(opt::Optimizer, ::MOI.RelativeGap)
    return _optimality_gap(primal, dual)
end

function MOI.get(opt::Optimizer, ::MOI.ResultCount)
    return _has_feasible_solution(opt) ? 1 : 0
end


text:
@testset "Maxlease and state_if_never_drilled" begin
    for unitprob ∈ (LeasedProblem, LeasedProblemContsDrill,)
        for wpp in (unitprob(3,4, 5,3,2), unitprob(3,4, 2,3,2), unitprob(3,4, 1,3,5), unitprob(3,16, 5,10,3), unitprob(3,4,1,-1,0))
            SS = ShaleDrillingModel.state_space_vector(wpp)
            @test maximum(ShaleDrillingModel._τrem(s) for s in SS) == ShaleDrillingModel.maxlease(wpp)
            for sidx ∈ 1:length(wpp)
                for d in ShaleDrillingModel.actionspace(wpp,sidx)
                    @test sprime(wpp,sidx,d) ∈ 1:length(wpp)
                end
            end

            # test state if never drilled...
            for state0 = 1:ShaleDrillingModel.end_ex0(wpp)
                statet = state0
                for t = 1:ShaleDrillingModel.end_ex0(wpp)
                    statet = sprime(wpp,statet,0)
                    ifno_drill = ShaleDrillingModel.state_if_never_drilled(wpp,state0,t)
                    @test statet == ifno_drill
                    # println("t,state0 = ($t,$state0): recursive sprime = $statet, state_if_never_drilled = $ifno_drill")
                end
            end

        end
    end
end

@testset "Drilling transitions are ok" begin
    for unitprob ∈ (LeasedProblem, LeasedProblemContsDrill,PerpetualProblem,)
        wp = unitprob(4, 4, 5, 3, 2)
        nS = length(wp)

        # infill: d1 == 0 cascades down
        for i = end_lrn(wp)+ShaleDrillingModel._nstates_per_D(wp) : ShaleDrillingModel._nstates_per_D(wp) : nS-1
            @test sprime(wp,i,0) == i
        end

        # infill: d1 == 1
        if isa(wp,LeasedProblemContsDrill)
            for i = end_lrn(wp)+1 : 2 : nS-2
                @test sprime(wp,i,0) == i+1
                for d in 1:ShaleDrillingModel._dmax(wp,i)
                    @test sprime(wp,i,d) == i+2*d
                end
            end
        end

        # initial lease term
        if isa(wp,PerpetualProblem)
            @test sprime(wp,1,0) == 1
        else
            for i = 1 :  end_ex1(wp)-1
                @test sprime(wp,i,0) == i+1
            end

            # last lease term
            for i = end_ex1(wp)+1 : end_ex0(wp)
                @test sprime(wp,i,0) == i+1
            end

            # any drilling during a primary term
            for i = 0+1 : end_ex0(wp)
                for d = 1:4
                    @test sprime(wp,i,d) == end_ex0(wp)+1+d
                end
            end
        end
    end
end

@testset "sprime is in 1:length(wp)?" begin
    for unitprob ∈ (ShaleDrillingModel.LeasedProblem, ShaleDrillingModel.LeasedProblemContsDrill, ShaleDrillingModel.PerpetualProblem,)
        for dmx in 3:4
            wp = unitprob(dmx, 4, 5, 3, 2)
            SS = ShaleDrillingModel.state_space_vector(wp)

            @test length(ShaleDrillingModel.inf_fm_lrn(wp)) == length(ShaleDrillingModel.exploratory_learning(wp))

            for (i,s) in enumerate(SS)
                for d in  ShaleDrillingModel.actionspace(wp,i)
                    sp = sprime(wp,i,d)
                    s = SS[i]  # this state
                    t = SS[sp] # next state
                    sp ∈ 1:length(wp) || println("dmax = $dmx d = $d, i = $i with si = $(SS[i]), sprime = $(sprime(wp,i,d))")
                    @test sp ∈ 1:length(wp)
                    @test t.D == s.D + d

                    @test SS[i] == state(wp,i)
                    @test SS[sp] == state(wp,sp)

                    if  i < length(wp) && i ∉ ShaleDrillingModel.ind_lrn(wp) && isa(wp,LeasedProblemContsDrill)
                        @test t.d1 == Int(d > 0)
                    end

                    if i ∈ ShaleDrillingModel.ind_ex1(wp)
                        @test t.τ1 == ( d == 0 ? s.τ1 -1 : -1)
                        @test t.τ0 == ( d == 0 ? s.τ0 - (i == end_ex1(wp)) : -1)
                    end

                    if i ∈ ShaleDrillingModel.ind_ex0(wp)
                        @test t.τ1 ==  s.τ1 == -1
                        if isa(wp, PerpetualProblem)
                            @test t.τ0 == -1
                        else
                            @test t.τ0 == ( d == 0 ? s.τ0 - 1 : -1)
                        end
                    end

                end
            end
        end
    end
end

@testset "make State Space" begin

    let dmx = 4,
        Dmx = 4,
        τ0mx = 5,
        τ1mx = 3,
        emx = 2

        for unitprob ∈ (ShaleDrillingModel.LeasedProblem, ShaleDrillingModel.LeasedProblemContsDrill, ShaleDrillingModel.PerpetualProblem,)
            wp = unitprob(dmx, Dmx, τ0mx, τ1mx, emx)

            SS = ShaleDrillingModel.state_space_vector(wp)

            SS[ 1:end_ex1(wp) ]                                # Exploratory 1
            SS[ 1:end_ex0(wp) ]                                # Exploratory 0
            SS[  end_ex0(wp) .+ (1:dmx+1)]                     # Exploratory Terminal + learning update
            # SS[  (τ0mx + τ1mx + dmx+3) .+ (1:2*Dmx-2)]   # Infill drilling
            # SS[  (τ0mx + τ1mx + dmx+3)  +   2*Dmx-2+1]   # Terminal

            nS = length(SS)
            [ShaleDrillingModel.actionspace(wp,i) for i in 1:nS]

            a = SS, [
                (i,
                ShaleDrillingModel._regime(wp,i),
                ShaleDrillingModel._τ1(wp,i),
                ShaleDrillingModel._τ0(wp,i),
                ShaleDrillingModel._D(wp,i),
                ShaleDrillingModel._d1(wp,i),
                collect(sprime(wp,i,d) for d in ShaleDrillingModel.actionspace(wp,i))
                ) for i in 1:nS]

            idxs = [ShaleDrillingModel.ind_exp(wp)..., ShaleDrillingModel.ind_inf(wp)..., end_inf(wp)..., ShaleDrillingModel.ind_lrn(wp)...]

            @test idxs ⊆ 1:length(wp)
            @test 1:length(wp) ⊆ idxs

            # test that we get back the state we want to get.
            for i in [ShaleDrillingModel.ind_exp(wp)..., ShaleDrillingModel.ind_inf(wp)..., end_inf(wp)...,]
                st = SS[i]
                i_of_st = state_idx(wp, st.τ1, st.τ0, st.D, st.d1)
                @test i_of_st == i
                @test i ∈ ShaleDrillingModel.s_of_D(wp, st.D)
            end

            for s in ShaleDrillingModel.ind_inf(wp)
                @test ShaleDrillingModel._horizon(wp,s) ∈ (:Infinite, :Finite)
            end

            for s in ShaleDrillingModel.ind_exp(wp)
                if isa(wp,PerpetualProblem)
                    @test ShaleDrillingModel._horizon(wp,s) == :Infinite
                else
                    @test ShaleDrillingModel._horizon(wp,s) == :Finite
                end
            end

            for s in ShaleDrillingModel.ind_lrn(wp)
                @test ShaleDrillingModel._horizon(wp,s) ∈ (:Terminal, :Learning)
            end
        end
    end
end

@testset "Checking on dims of dEVσ" begin
    for unitprob ∈ (LeasedProblem, LeasedProblemContsDrill,)
        for wpp in (unitprob(3,4, 5,3,2), unitprob(3,4, 2,3,2), unitprob(3,4, 1,3,5), unitprob(3,16, 5,10,3), unitprob(3,4,1,-1,0))
            nsexp = ShaleDrillingModel._nSexp(wpp)
            for sidx ∈ 1:length(wpp)
                if sidx <= end_ex0(wpp)
                    for d in ShaleDrillingModel.actionspace(wpp,sidx)
                        @test sprime(wpp,sidx,d) ∈ 1:nsexp
                    end
                end
            end
        end
    end
end


text:
include("./unit/runtests.jl")
include("./integration/runtests.jl")

#TODO I'm not familiar with the tests yet, so just dumping it here to get us started.
include("./integration/testStandardAPI.jl")