text:
try
# construct absolute path
depsdir  = splitdir(Base.source_path())[1]
builddir = joinpath(depsdir,"builds")
srcdir   = joinpath(depsdir,"src")

println("=== Building Multigrid ===")
println("depsdir  = $depsdir")
println("builddir = $builddir")
println("srcdir   = $srcdir")

if !isdir(builddir)
	println("creating build directory")
	mkdir(builddir)
	if !isdir(builddir)
		error("Could not create build directory")
	end
end

@static if Sys.isunix()
	src1 = joinpath(srcdir,"Vanka.c")
	src2 = joinpath(srcdir,"parRelax.c");
	src3 = joinpath(srcdir,"parLU.cpp");
	outfile1 = joinpath(builddir,"Vanka.so")
	outfile2 = joinpath(builddir,"parRelax.so")
	outfile3 = joinpath(builddir,"parLU.so")
	println("gcc version")
	run(`gcc --version`)
	run(`gcc -O3 -fPIC -cpp -fopenmp -shared  $src1 -o $outfile1`)
	run(`gcc -O3 -fPIC -cpp -fopenmp -shared  $src2 -o $outfile2`)
	run(`g++ -O3 -fPIC -cpp -fopenmp -shared  $src3 -o $outfile3`)
end

@static if Sys.iswindows()
	src1 = joinpath(srcdir,"Vanka.c")
	src2 = joinpath(srcdir,"parRelax.c");
	src3 = joinpath(srcdir,"parLU.cpp");
	outfile1 = joinpath(builddir,"Vanka.dll")
	outfile2 = joinpath(builddir,"parRelax.dll")
	outfile3 = joinpath(builddir,"parLU.dll")
	println("gcc version")
	run(`gcc --version`)
	run(`gcc -O3 -cpp -fopenmp -shared -DBUILD_DLL  $src1 -o $outfile1`)
	run(`gcc -O3 -cpp -fopenmp -shared -DBUILD_DLL  $src2 -o $outfile2`)
	run(`g++ -O3 -cpp -fopenmp -shared -DBUILD_DLL  $src3 -o $outfile3`)
end
catch
	@warn "Multigrid::build: Unable to build Multigrid"
end


text:
const SKIP_PM_VALIDATION = false

"""
System

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

```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 and when add_component! is called.
  Throws InvalidRange if an error is found.
- `time_series_in_memory::Bool=false`: Store time series data in memory instead of HDF5.
- `configpath::String`: specify path to validation config file
"""
struct System <: PowerSystemType
    data::IS.SystemData
    basepower::Float64 # [MVA]
    frequency::Float64 # [Hz]
    bus_numbers::Set{Int}
    runchecks::Bool
    internal::IS.InfrastructureSystemsInternal

    function System(data, basepower, internal; kwargs...)
        bus_numbers = Set{Int}()
        frequency = get(kwargs, :frequency, DEFAULT_SYSTEM_FREQUENCY)
        runchecks = get(kwargs, :runchecks, true)
        sys = new(data, basepower, frequency, bus_numbers, runchecks, internal)
        return sys
    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; kwargs...)
end

"""Construct a `System` from `InfrastructureSystems.SystemData`"""
function System(data, basepower; kwargs...)
    internal = get(kwargs, :internal, IS.InfrastructureSystemsInternal())
    return System(data, basepower, internal; 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, IS.InvalidRange)
                error_detected = true
            else
                rethrow()
            end
        end
    end

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

    runchecks = get(kwargs, :runchecks, true)

    if error_detected
        throw(IS.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

"""
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

"""
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

"""
Return a user-modifiable dictionary to store extra information.
"""
get_ext(sys::System) = IS.get_ext(sys.internal)

"""
Return the system's base power.
"""
get_basepower(sys::System) = sys.basepower

"""
Return the system's frequency.
"""
get_frequency(sys::System) = sys.frequency

"""
Clear any value stored in ext.
"""
clear_ext(sys::System) = IS.clear_ext(sys.internal)

"""
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)
    check_for_services_on_addition(sys, component)

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

    if sys.runchecks && !validate_struct(sys, component)
        throw(IS.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_service!(sys::System, service::Service, contributing_devices; kwargs...)

Similar to [`add_component!`](@ref) but for services.

# Arguments
- `sys::System`: system
- `service::Service`: service to add
- `contributing_devices`: Must be an iterable of type Device
"""
function add_service!(sys::System, service::Service, contributing_devices; kwargs...)
    if sys.runchecks && !validate_struct(sys, service)
        throw(InvalidValue("Invalid value for $service"))
    end

    for device in contributing_devices
        device_type = typeof(device)
        if !(device_type <: Device)
            throw(ArgumentError("contributing_devices must be of type Device"))
        end

        name = get_name(device)
        sys_device = get_component(device_type, sys, name)
        if isnothing(sys_device)
            throw(ArgumentError("device $device_type $name is not part of the system"))
        end
    end

    # Since this isn't atomic, order is important. Add to system before adding to devices.
    IS.add_component!(sys.data, service; kwargs...)

    for device in contributing_devices
        add_service!(device, service)
    end
end

"""
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.
- `resolution::DateTime.Period=nothing`: skip forecast that don't match this resolution.
"""
function add_forecasts!(sys::System, metadata_file::AbstractString; resolution = nothing)
    return IS.add_forecasts!(Component, sys.data, metadata_file; resolution = resolution)
end

"""
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 Area
        uuids = Set{UUIDs.UUID}()
        for bus in _get_buses(data, component)
            push!(uuids, IS.get_uuid(bus))
        end
        for component_ in (
            load for load in IS.get_components(ElectricLoad, data) if IS.get_uuid(get_bus(
                load,
            )) in uuids
        )
            IS.add_forecast!(data, component_, forecast, ts_data)
        end
    else
        IS.add_forecast!(data, component, forecast, ts_data)
    end
end

"""
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

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

"""
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 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 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 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

"""
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 = collect(PowerSystems.get_components(Generator, sys))
```

See also: [`iterate_components`](@ref)
"""
function get_components(::Type{T}, sys::System) where {T <: Component}
    # TODO: Verify that return type annotation is not required
    return IS.get_components(T, sys.data)
end

function _get_components_by_name(abstract_types, data::IS.SystemData, name::AbstractString)
    _components = []
    for subtype in abstract_types
        component = IS.get_component(subtype, data, name)
        if !isnothing(component)
            push!(_components, component)
        end
    end

    return _components
end

"""
    get_components_by_name(
                           ::Type{T},
                           sys::System,
                           name::AbstractString
                          ) 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,
) where {T <: Component}
    # TODO: Verify that return type annotation is not required
    return IS.get_components_by_name(T, sys.data, name)
end

"""
    get_contributing_devices(sys::System, service::Service)

Return a vector of devices contributing to the service.
"""
function get_contributing_devices(sys::System, service::T) where {T <: Service}
    if isnothing(get_component(T, sys, get_name(service)))
        throw(ArgumentError("service $(get_name(service)) is not part of the system"))
    end

    return [x for x in get_components(Device, sys) if has_service(x, service)]
end

struct ServiceContributingDevices
    service::Service
    contributing_devices::Vector{Device}
end

const ServiceContributingDevicesKey = NamedTuple{(:type, :name), Tuple{DataType, String}}
const ServiceContributingDevicesMapping =
    Dict{ServiceContributingDevicesKey, ServiceContributingDevices}

"""
    get_contributing_device_mapping(sys::System)

Return an instance of ServiceContributingDevicesMapping.
"""
function get_contributing_device_mapping(sys::System)
    services = ServiceContributingDevicesMapping()
    for service in get_components(Service, sys)
        uuid = IS.get_uuid(service)
        devices = ServiceContributingDevices(service, Vector{Device}())
        for device in get_components(Device, sys)
            for _service in get_services(device)
                if IS.get_uuid(_service) == uuid
                    push!(devices.contributing_devices, device)
                    break
                end
            end

            key = ServiceContributingDevicesKey((typeof(service), get_name(service)))
            services[key] = devices
        end
    end

    return services
end

"""
    get_components_in_aggregation_topology(
        ::Type{T},
        sys::System,
        aggregator::AggregationTopology,
    ) where {T <: StaticInjection}

Return a vector of components with buses in the AggregationTopology.
"""
function get_components_in_aggregation_topology(
    ::Type{T},
    sys::System,
    aggregator::AggregationTopology,
) where {T <: StaticInjection}
    buses = Set{String}((get_name(x) for x in get_buses(sys, aggregator)))
    components = Vector{T}()
    for component in get_components(T, sys)
        bus = get_bus(component)
        bus_name = get_name(bus)
        if bus_name in buses
            push!(components, component)
        end
    end

    return components
end

"""
    get_aggregation_topology_mapping(::Type{T}, sys::System) where {T<:AggregationTopology}

Return a mapping of AggregationTopology name to vector of buses within it.
"""
function get_aggregation_topology_mapping(
    ::Type{T},
    sys::System,
) where {T <: AggregationTopology}
    mapping = Dict{String, Vector{Bus}}()
    accessor_func = get_aggregation_topology_accessor(T)
    for bus in get_components(Bus, sys)
        aggregator = accessor_func(bus)
        name = get_name(aggregator)
        buses = get(mapping, name, nothing)
        if isnothing(buses)
            mapping[name] = Vector{Bus}([bus])
        else
            push!(buses, bus)
        end
    end

    return mapping
end

"""
    get_buses(sys::System, aggregator::AggregationTopology)

Return a vector of buses contained within the AggregationTopology.
"""
function get_buses(sys::System, aggregator::AggregationTopology)
    return _get_buses(sys.data, aggregator)
end

function _get_buses(data::IS.SystemData, aggregator::T) where {T <: AggregationTopology}
    accessor_func = get_aggregation_topology_accessor(T)
    buses = Vector{Bus}()
    for bus in IS.get_components(Bus, data)
        _aggregator = accessor_func(bus)
        if IS.get_uuid(_aggregator) == IS.get_uuid(aggregator)
            push!(buses, bus)
        end
    end

    return buses
end

"""
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 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 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 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

"""
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

"""
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

"""
Return true if forecasts are stored contiguously.

Throws ArgumentError if there are no forecasts stored.
"""
function are_forecasts_contiguous(sys::System)
    return IS.are_forecasts_contiguous(sys.data)
end

"""
"""
function are_forecasts_contiguous(component::Component)
    return IS.are_forecasts_contiguous(component)
end

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

# Arguments
- `sys::System`: System.
- `interval::Dates.Period`: Amount of time in between each initial time.
- `horizon::Int`: Length of each forecast array.
- `initial_time::Union{Nothing, Dates.DateTime}=nothing`: Start with this time. If nothing,
  use the first initial time.
"""
function generate_initial_times(
    sys::System,
    interval::Dates.Period,
    horizon::Int;
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
)
    return IS.generate_initial_times(
        sys.data,
        interval,
        horizon;
        initial_time = initial_time,
    )
end

"""
Generate initial times for a component.
"""
function generate_initial_times(
    component::IS.InfrastructureSystemsType,
    interval::Dates.Period,
    horizon::Int;
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
)
    return IS.generate_initial_times(
        component,
        interval,
        horizon;
        initial_time = initial_time,
    )
end

"""
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

"""
Return a forecast for a subset of the time series range stored for these parameters.
"""
function get_forecast(
    ::Type{T},
    component::IS.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

"""
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(PowerSystems, component, forecast)
end

"""
Return sorted forecast initial times.
"""
function get_forecast_initial_times(sys::System)
    return IS.get_forecast_initial_times(sys.data)
end

"""
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

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

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

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

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

"""
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

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

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

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

"""
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

"""
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)
    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)
    internal = IS.convert_type(InfrastructureSystemsInternal, raw.internal)
    sys = System(data, float(raw.basepower); internal = internal, 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.
    # This means that the restrictions in add_component! are not applied here.

    # Maintain a lookup of UUID to component because some component types encode
    # composed types as UUIDs instead of actual types.
    component_cache = Dict{Base.UUID, Component}()

    components_as_uuids = [Bus]
    for component_as_uuid in components_as_uuids
        for component in IS.get_components_raw(IS.SystemData, component_as_uuid, raw)
            comp = IS.convert_type(component_as_uuid, component)
            IS.add_component!(data, comp)
            component_cache[IS.get_uuid(comp)] = comp
        end
    end

    # Skip Devices this round because they have Services.
    for c_type_sym in IS.get_component_types_raw(IS.SystemData, raw)
        c_type = _get_component_type(c_type_sym)
        (c_type in components_as_uuids || c_type <: Device) && 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 Devices.
    for c_type_sym in IS.get_component_types_raw(IS.SystemData, raw)
        c_type = _get_component_type(c_type_sym)
        if c_type <: Device
            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 _get_component_type(component_type::Symbol)
    base_name = IS.strip_module_name(string(component_type))
    type_name, parameters = IS.separate_type_and_parameter_types(base_name)
    c_type = getfield(PowerSystems, Symbol(type_name))
    if !isempty(parameters)
        parametric_types = [getfield(PowerSystems, Symbol(x)) for x in parameters]
        c_type = c_type{parametric_types...}
    end

    return c_type
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 composed 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

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

"""
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

"""
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

check_for_services_on_addition(sys::System, component::Component) = nothing

function check_for_services_on_addition(sys::System, component::Device)
    if supports_services(component) && length(get_services(component)) > 0
        throw(ArgumentError("type Device cannot be added with services"))
    end
    return
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.
"""
check_component_addition(sys::System, component::Component) = nothing

"""
Refer to docstring for check_component_addition!
"""
handle_component_addition!(sys::System, component::Component) = nothing

handle_component_removal!(sys::System, component::Component) = nothing

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

function handle_component_removal!(sys::System, device::Device)
    # This may have to be refactored if handle_component_removal! needs to be implemented
    # for a subtype.
    clear_services!(device)
end

function handle_component_removal!(sys::System, service::Service)
    for device in get_components(Device, sys)
        _remove_service!(device, service)
    end
end

"""
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)
    match = true

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

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

    return match
end

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

    return IS.SystemData(;
        validation_descriptor_file = validation_descriptor_file,
        time_series_in_memory = time_series_in_memory,
    )
end

"""
Converts a Line component to a MonitoredLine component and replaces the original in the
system
"""
function convert_component!(
    linetype::Type{MonitoredLine},
    line::Line,
    sys::System;
    kwargs...,
)
    new_line = linetype(
        line.name,
        line.available,
        line.activepower_flow,
        line.reactivepower_flow,
        line.arc,
        line.r,
        line.x,
        line.b,
        (from_to = line.rate, to_from = line.rate),
        line.rate,
        line.anglelimits,
        line.services,
        line.ext,
        line.forecasts,
        line.internal,
    )
    remove_component!(sys, line)
    add_component!(sys, new_line)
end

"""
Converts a MonitoredLine component to a Line component and replaces the original in the
system
"""
function convert_component!(
    linetype::Type{Line},
    line::MonitoredLine,
    sys::System;
    kwargs...,
)
    force = get(kwargs, :force, false)
    if force
        @warn("Possible data loss converting from $(typeof(line)) to $linetype")
    else
        error("Possible data loss converting from $(typeof(line)) to $linetype, add `force = true` to convert anyway.")
    end

    new_line = linetype(
        line.name,
        line.available,
        line.activepower_flow,
        line.reactivepower_flow,
        line.arc,
        line.r,
        line.x,
        line.b,
        line.rate,
        line.anglelimits,
        line.services,
        line.ext,
        line.forecasts,
        line.internal,
    )
    remove_component!(sys, line)
    add_component!(sys, new_line)
end

ext: jl


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("MLJTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# [MLJ.jl]: https://github.com/alan-turing-institute/MLJ.jl# [RDatasets.jl]: https://github.com/JuliaStats/RDatasets.jl# [NearestNeighbors.jl]: https://github.com/KristofferC/NearestNeighbors.jl## ## Tuning a single hyperparameter## In MLJ, tuning is implemented as a model wrapper.# After wrapping a model in a _tuning strategy_ (e.g. cross-validation) and binding the wrapped model to data in a _machine_, fitting the machine initiates a search for optimal model hyperparameters.## Let's use a decision tree classifier and tune the maximum depth of the tree.# As usual, start by loading data and the model
using MLJ, PrettyPrinting

X, y = @load_iris
@load DecisionTreeClassifier

# ### Specifying a range of value## To specify a range of value, you can use the `range` function:
dtc = DecisionTreeClassifier()
r   = range(dtc, :max_depth, lower=1, upper=5)

# As you can see, the range function takes a model (`dtc`), a symbol for the hyperparameter of interest (`:max_depth`) and indication of how to samples values.# For hyperparameters of type `<:Real`, you should specify a range of values as done above.# For hyperparameters of other type (e.g. `Symbol`), you should use the `values=...` keyword.## Once a range of values has been defined, you can then wrap the model in a `TunedModel` specifying the tuning strategy:
tm = TunedModel(model=dtc, ranges=[r, ], measure=cross_entropy)

# ### Fitting and inspecting a tuned model## To fit a tuned model, you can use the usual syntax:
m = machine(tm, X, y)
fit!(m)

# In order to inspect the best model, you can use the function `fitted_params` on the machine and inspect the `best_model` field:
fitted_params(m).best_model.max_depth

# Note that here we have tuned a probabilistic model and consequently used a probabilistic measure for the tuning.# We could also have decided we only cared about the mode and the misclassification rate, to do this, just use `operation=predict_mode` in the tuned model:
tm = TunedModel(model=dtc, ranges=r, operation=predict_mode,
                measure=misclassification_rate)
m = machine(tm, X, y)
fit!(m)
fitted_params(m).best_model.max_depth

# Let's check the misclassification rate for the best model:
r = report(m)
r.best_result

# Anyone wants plots? of course:
using PyPlot
figure(figsize=(8,6))
res = r.plotting # contains all you need for plotting
plot(res.parameter_values, res.measurements, ls="none", marker="o")

xticks(1:5, fontsize=12)
yticks(fontsize=12)
xlabel("Maximum depth", fontsize=14)
ylabel("Misclassification rate", fontsize=14)
ylim([0, 1])

# \figalt{hyperparameter heatmap}{A-model-tuning-hpt}
# ## Tuning nested hyperparameters
# Let's generate simple dummy regression data
X = (x1=rand(100), x2=rand(100), x3=rand(100))
y = 2X.x1 - X.x2 + 0.05 * randn(100);

# Let's then build a simple ensemble model with decision tree regressors:
dtr = @load DecisionTreeRegressor
forest = EnsembleModel(atom=dtr)

# Such a model has *nested* hyperparameters in that the ensemble has hyperparameters (e.g. the `:bagging_fraction`) and the atom has hyperparameters (e.g. `:n_subfeatures` or `:max_depth`).# You can see this by inspecting the parameters using `params`:
params(forest) |> pprint

# Range for nested hyperparameters are specified using dot syntax, the rest is done in much the same way as before:
r1 = range(forest, :(atom.n_subfeatures), lower=1, upper=3)
r2 = range(forest, :bagging_fraction, lower=0.4, upper=1.0)
tm = TunedModel(model=forest, tuning=Grid(resolution=12),
                resampling=CV(nfolds=6), ranges=[r1, r2],
                measure=rms)
m = machine(tm, X, y)
fit!(m);

# A useful function to inspect a model after fitting it is the `report` function which collects information on the model and the tuning, for instance you can use it to recover the best measurement:
r = report(m)
r.best_result

# Let's visualise this
figure(figsize=(8,6))

res = r.plotting

vals_sf = res.parameter_values[:, 1]
vals_bf = res.parameter_values[:, 2]

tricontourf(vals_sf, vals_bf, res.measurements)
xlabel("Number of sub-features", fontsize=14)
ylabel("Bagging fraction", fontsize=14)
xticks([1, 2, 3], fontsize=12)
yticks(fontsize=12)

# \figalt{Hyperparameter heatmap}{A-model-tuning-hm.svg}

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
"""
    FieldBoundaryConditions

An alias for `NamedTuple{(:x, :y, :z)}` that represents a set of three `CoordinateBoundaryCondition`s
applied to a field along x, y, and z.
"""
const FieldBoundaryConditions = NamedTuple{(:x, :y, :z)}

"""
    FieldBoundaryConditions(x, y, z)

Construct a `FieldBoundaryConditions` using a `CoordinateBoundaryCondition` for each of the
`x`, `y`, and `z` coordinates.
"""
FieldBoundaryConditions(x, y, z) = FieldBoundaryConditions((x, y, z))

"""
    DefaultBoundaryCondition(::Union{Type{Periodic}, Type{Flat}}, loc)

Returns [`PeriodicBoundaryCondition`](@ref).
"""
DefaultBoundaryCondition(::Union{Type{Grids.Periodic}, Type{Flat}}, loc) = PeriodicBoundaryCondition()

"""
    DefaultBoundaryCondition(::Type{Bounded}, ::Type{Cell})

Returns [`NoFluxBoundaryCondition`](@ref).
"""
DefaultBoundaryCondition(::Type{Bounded}, ::Type{Cell}) = NoFluxBoundaryCondition()

"""
    DefaultBoundaryCondition(::Type{Bounded}, ::Type{Face})

Returns [`ImpenetrableBoundaryCondition`](@ref).
"""
DefaultBoundaryCondition(::Type{Bounded}, ::Type{Face}) = ImpenetrableBoundaryCondition()

function validate_bcs(topology, left_bc, right_bc, default_bc, left_name, right_name, dir)
    if topology isa Periodic && (left_bc != default_bc || right_bc != default_bc)
        e = "Cannot specify $left_name or $right_name boundary conditions with " *
            "$topology topology in $dir-direction."
        throw(ArgumentError(e))
    end
    return true
end

relocate_function_boundary_condition(bc, args...) = bc # fallback

relocate_function_boundary_condition(bc::BoundaryCondition{C, <:BoundaryFunction}, boundary, loc1, loc2) where C =
    BoundaryCondition(C, BoundaryFunction{boundary, loc1, loc2}(bc.condition.func; parameters=bc.condition.parameters))

"""
    FieldBoundaryConditions(grid, loc;   east = DefaultBoundaryCondition(topology(grid, 1), loc[1]),
                                         west = DefaultBoundaryCondition(topology(grid, 1], loc[1]),
                                        south = DefaultBoundaryCondition(topology(grid, 2), loc[2]),
                                        north = DefaultBoundaryCondition(topology(grid, 2), loc[2]),
                                       bottom = DefaultBoundaryCondition(topology(grid, 3), loc[3]),
                                          top = DefaultBoundaryCondition(topology(grid, 3), loc[3]))

Construct `FieldBoundaryConditions` for a field with location `loc` (a 3-tuple of `Face` or `Cell`)
defined on `grid`.

Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)` and `loc`.
"""
function FieldBoundaryConditions(grid, loc;   east = DefaultBoundaryCondition(topology(grid, 1), loc[1]),
                                              west = DefaultBoundaryCondition(topology(grid, 1), loc[1]),
                                             south = DefaultBoundaryCondition(topology(grid, 2), loc[2]),
                                             north = DefaultBoundaryCondition(topology(grid, 2), loc[2]),
                                            bottom = DefaultBoundaryCondition(topology(grid, 3), loc[3]),
                                               top = DefaultBoundaryCondition(topology(grid, 3), loc[3]))

    east = relocate_function_boundary_condition(east, :x, loc[2], loc[3])
    west = relocate_function_boundary_condition(west, :x, loc[2], loc[3])
    south = relocate_function_boundary_condition(south, :y, loc[1], loc[3])
    north = relocate_function_boundary_condition(north, :y, loc[1], loc[3])
    bottom = relocate_function_boundary_condition(bottom, :z, loc[1], loc[2])
    top = relocate_function_boundary_condition(top, :z, loc[1], loc[2])

    TX, TY, TZ = topology(grid)
    x_default_bc = DefaultBoundaryCondition(topology(grid, 1), loc[1])
    y_default_bc = DefaultBoundaryCondition(topology(grid, 2), loc[2])
    z_default_bc = DefaultBoundaryCondition(topology(grid, 3), loc[3])

    validate_bcs(TX, west,   east, x_default_bc, :west,   :east, :x)
    validate_bcs(TY, south, north, y_default_bc, :south, :north, :y)
    validate_bcs(TZ, bottom,  top, z_default_bc, :bottom,  :top, :z)

    x = CoordinateBoundaryConditions(west, east)
    y = CoordinateBoundaryConditions(south, north)
    z = CoordinateBoundaryConditions(bottom, top)

    return FieldBoundaryConditions(x, y, z)
end

"""
    UVelocityBoundaryConditions(grid;   east = DefaultBoundaryCondition(topology(grid, 1), Face),
                                        west = DefaultBoundaryCondition(topology(grid, 1), Face),
                                       south = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                       north = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                      bottom = DefaultBoundaryCondition(topology(grid, 3), Cell),
                                         top = DefaultBoundaryCondition(topology(grid, 3), Cell))

Construct `FieldBoundaryConditions` for the `u`-velocity field on `grid`.
Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`.
"""
UVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Face, Cell, Cell); user_defined_bcs...)

"""
    VVelocityBoundaryConditions(grid;   east = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                        west = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                       south = DefaultBoundaryCondition(topology(grid, 2), Face),
                                       north = DefaultBoundaryCondition(topology(grid, 2), Face),
                                      bottom = DefaultBoundaryCondition(topology(grid, 3), Cell),
                                         top = DefaultBoundaryCondition(topology(grid, 3), Cell))

Construct `FieldBoundaryConditions` for the `v`-velocity field on `grid`.
Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`.
"""
VVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Face, Cell); user_defined_bcs...)

"""
    WVelocityBoundaryConditions(grid;   east = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                        west = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                       south = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                       north = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                      bottom = DefaultBoundaryCondition(topology(grid, 3), Face),
                                         top = DefaultBoundaryCondition(topology(grid, 3), Face))

Construct `FieldBoundaryConditions` for the `w`-velocity field on `grid`.
Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`.
"""
WVelocityBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Face); user_defined_bcs...)

"""
    TracerBoundaryConditions(grid;   east = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                     west = DefaultBoundaryCondition(topology(grid, 1), Cell),
                                    south = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                    north = DefaultBoundaryCondition(topology(grid, 2), Cell),
                                   bottom = DefaultBoundaryCondition(topology(grid, 3), Cell),
                                      top = DefaultBoundaryCondition(topology(grid, 3), Cell))

Construct `FieldBoundaryConditions` for a tracer field on `grid`.
Boundary conditions on `x`-, `y`-, and `z`-boundaries are specified via
keyword arguments:

    * `west` and `east` for the `-x` and `+x` boundary;
    * `south` and `north` for the `-y` and `+y` boundary;
    * `bottom` and `top` for the `-z` and `+z` boundary.

Default boundary conditions depend on `topology(grid)`. See `DefaultBoundaryCondition`.
"""
TracerBoundaryConditions(grid; user_defined_bcs...) = FieldBoundaryConditions(grid, (Cell, Cell, Cell); user_defined_bcs...)

const PressureBoundaryConditions = TracerBoundaryConditions
const DiffusivityBoundaryConditions = TracerBoundaryConditions


text:
using ParametricSurfaces # for generating some simple point distributions
using LoopVectorization  # for automatically vectorizing some typical kernels

# some point distributions for testing
function cube_uniform_surface_mesh(dx)
    geo = ParametricSurfaces.Cube() # a cube with low corner at (0,0,0) and high corner (2,2,2)
    range = -1+dx/2:dx:1-dx/2
    pts   = Vector{Geometry.Point3D}()
    for ent in boundary(geo)
        for u in range, v in range
            push!(pts,ent((u,v)))
        end
    end
    return pts
end

function spheroid(N,r,center,zstretch)
    pts   = Vector{Geometry.Point3D}(undef,N)
    phi   = π*(3-sqrt(5)) # golden angle in radians
    for i = 1:N
        ytmp   = 1 - ((i-1)/(N-1))*2
        radius = sqrt(1-ytmp^2)
        theta    = phi*i
        x = cos(theta) * radius * r + center[1]
        y = ytmp * r + center[2]
        z = zstretch * sin(theta) * radius * r + center[3]
        pts[i] = SVector(x,y,z)
    end
    # Output
    return pts
end

sphere_uniform(N,r,center=(0,0,0)) = spheroid(N,r,center,1)

function cube_nonuniform_surface_mesh(dx,p=2)
    geo = ParametricSurfaces.Cube() # a cube with low corner at (0,0,0) and high corner (2,2,2)
    range = -1+dx/2:dx:1-dx/2
    pts   = Vector{Geometry.Point3D}()
    # degenerate change of variables mapping [0,1] onto itself with p
    # derivatives vanishing at each endpoint
    cov   = Integration.KressP(order=p)
    # map cov to act on [-1,1]
    χ     = (u) -> begin
        û = 0.5*(u+1)
        ŝ = cov(û)
        s = 2*ŝ - 1
    end
    for ent in boundary(geo)
        for u in range, v in range
            push!(pts,ent((χ(u),χ(v))))
        end
    end
    return pts
end

function cube_uniform_volume_mesh(dx)
    range = dx/2:dx:1-dx/2
    [SVector(a,b,c) for a in range, b in range, c in range] |> vec
end

function sphere_uniform_surface_mesh(dx)
    geo   = ParametricSurfaces.Sphere(radius=1) # a sphere of radius 1 centered at (0,0,0)
    range = -1+dx/2:dx:1-dx/2
    pts   = Vector{Geometry.Point3D}()
    for ent in boundary(geo)
        for u in range, v in range
            push!(pts,ent((u,v)))
        end
    end
    return pts
end

# define some typical kernels
struct HelmholtzKernel
    k::Float64
end

function (K::HelmholtzKernel)(x,y)
    d = norm(x-y)
    v = exp(im*K.k*d)/(4*π*d)
    return (!iszero(d))*v
end

IFGF.wavenumber(K::HelmholtzKernel) = K.k

function IFGF.near_interaction!(C,K::HelmholtzKernel,X,Y,σ,I,J)
    Tx = eltype(X)
    Ty = eltype(Y)
    @assert Tx <: SVector && Ty <: SVector
    Xm = reshape(reinterpret(eltype(Tx),X),3,:)
    Ym = reshape(reinterpret(eltype(Ty),Y),3,:)
    @views helmholtz3d_sl_vec!(C[I],Xm[:,I],Ym[:,J],σ[J],K.k)
end

function IFGF.transfer_factor(K::HelmholtzKernel,x,Y)
    yc  = IFGF.center(Y)
    yp  = IFGF.center(IFGF.parent(Y))
    d   = norm(x-yc)
    dp  = norm(x-yp)
    exp(im*K.k*(d-dp))*dp/d
end

function helmholtz3d_sl_vec!(C,X,Y,σ,k)
    m,n = size(X,2), size(Y,2)
    C_T = reinterpret(Float64, C)
    C_r = @views C_T[1:2:end,:]
    C_i = @views C_T[2:2:end,:]
    σ_T = reinterpret(Float64, σ)
    σ_r = @views σ_T[1:2:end,:]
    σ_i = @views σ_T[2:2:end,:]
    @turbo for j in 1:n
        for i in 1:m
            d2 = (X[1,i] - Y[1,j])^2
            d2 += (X[2,i] - Y[2,j])^2
            d2 += (X[3,i] - Y[3,j])^2
            d  = sqrt(d2)
            s, c = sincos(k * d)
            zr = inv(4π*d) * c
            zi = inv(4π*d) * s
            C_r[i] += (!iszero(d))*(zr*σ_r[j] - zi*σ_i[j])
            C_i[i] += (!iszero(d))*(zi*σ_r[j] + zr*σ_i[j])
        end
    end
    return C
end

struct LaplaceKernel
end

function (K::LaplaceKernel)(x,y)
    d = norm(x-y)
    v = 1/(4*π*d)
    return (!iszero(d))*v
end

function IFGF.near_interaction!(C,K::LaplaceKernel,X,Y,σ,I,J)
    Xm = reshape(reinterpret(Float64,X),3,:)
    Ym = reshape(reinterpret(Float64,Y),3,:)
    @views laplace3d_sl_vec!(C[I],Xm[:,I],Ym[:,J],σ[J])
end

function laplace3d_sl_vec!(C,X,Y,σ)
    m,n = size(X,2), size(Y,2)
    @turbo for j in 1:n
        for i in 1:m
            d2 = (X[1,i] - Y[1,j])^2
            d2 += (X[2,i] - Y[2,j])^2
            d2 += (X[3,i] - Y[3,j])^2
            # fast invsqrt code taken from here
            # https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-julia-8.html
            invd    = @fastmath Float64(1 / sqrt(Float32(d2)))
            invd = 1.5invd - 0.5d2 * invd * (invd * invd)
            C[i] += (!iszero(d2))*(inv(4π)*invd*σ[j])
            # C[i] += inv(4π*sqrt(d2))*σ[j] # significalty slower
        end
    end
    return C
end

struct MaxwellKernel
    k::Float64
end

function (K::MaxwellKernel)(x,y)
    r   = x - y
    d   = norm(r)
    # helmholtz greens function
    g   = exp(im*K.k*d)/(4π*d)
    gp  = im*K.k*g - g/d
    gpp = im*K.k*gp - gp/d + g/d^2
    RRT = r*transpose(r) # rvec ⊗ rvecᵗ
    G   = g*LinearAlgebra.I + 1/K.k^2*(gp/d*LinearAlgebra.I + (gpp/d^2 - gp/d^3)*RRT)
    return (!iszero(d))*G
end

IFGF.wavenumber(K::MaxwellKernel) = K.k

# function IFGF.centered_factor(K::MaxwellKernel,x,Y)
#     yc = IFGF.center(Y)
#     r  = x-yc
#     d = norm(r)
#     g   = exp(im*k*d)/(4π*d)
#     return g    # works fine!
# end

# function IFGF.transfer_factor(K::HelmholtzKernel,x,Y)
#     yc  = IFGF.center(Y)
#     yp  = IFGF.center(IFGF.parent(Y))
#     d   = norm(x-yc)
#     dp  = norm(x-yp)
#     exp(im*K.k*(d-dp))*dp/d
# end

# function helmholtz3d_sl_vec!(C,X,Y,σ,k)
#     m,n = size(X,2), size(Y,2)
#     C_T = reinterpret(Float64, C)
#     C_r = @views C_T[1:2:end,:]
#     C_i = @views C_T[2:2:end,:]
#     σ_T = reinterpret(Float64, σ)
#     σ_r = @views σ_T[1:2:end,:]
#     σ_i = @views σ_T[2:2:end,:]
#     @turbo for j in 1:n
#         for i in 1:m
#             d2 = (X[1,i] - Y[1,j])^2
#             d2 += (X[2,i] - Y[2,j])^2
#             d2 += (X[3,i] - Y[3,j])^2
#             d  = sqrt(d2)
#             s, c = sincos(k * d)
#             zr = inv(4π*d) * c
#             zi = inv(4π*d) * s
#             C_r[i] += (!iszero(d))*(zr*σ_r[j] - zi*σ_i[j])
#             C_i[i] += (!iszero(d))*(zi*σ_r[j] + zr*σ_i[j])
#         end
#     end
#     return C
# end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
import Base.show
import Base.zero
import Base.isless
type Route
    cities::Array{City}
    score::Float32
end

Route(cities::Array{City}) = Route(copy(cities), score(cities))
Route(route::Route) = Route(copy(route.cities), route.score)

Base.length(r::Route) = length(r.cities)
Base.getindex(r::Route, idx) = getindex(r.cities, idx)
Base.setindex!(r::Route, val, idx) = setindex!(r.cities, val, idx)
Base.endof(r::Route) = endof(r.cities)

zero(::Type{Route}) = Route(zeros(City, 0), -1)

function Base.show(io::IO, route::Route)
    print(io, "Route(")
    print(io, map(city->city.id, route.cities))
    print(io, " Total length: $(route.score))")
end

function score(cities::Array{City})
    total = 0.0f0
    for i in 0:length(cities) - 1
        total += distance(cities[i % length(cities) + 1], cities[(i + 1) % length(cities) + 1])
    end
    return total
end

function score(cities::Array{City}, grid::Grid)
    total = 0.0f0
    for i in 0:length(cities) - 1
        total += distance(cities[i % length(cities) + 1], cities[(i + 1) % length(cities) + 1], grid)
    end
    return total
end

function score!(route::Route)
    route.score = score(route.cities)
    route.score
end

function score!(route::Route, grid::Grid)
    route.score = score(route.cities, grid)
    route.score
end

isless(a::Route, b::Route) = a.score < b.score
function total_distance(route::Route)
    total = 0.0f0
    for i in 0:length(route) - 1
        total += distance(route[i % length(route) + 1], route[(i + 1) % length(route) + 1])
    end
    return total
end

function report(route::Route)
    print("<")
    for city in route[1:end-1]
        print("$(city.id),")
    end
    print("$(route[end].id)>\n")
    println("Total length: $(total_distance(route))")
end


text:
@testset "Inner product" begin
    @test v1 ⋅ v1 == 1v
    @test v1 ⋅ v2 == 0v
    @test v12 ⋅ v == 0v
    @test A ⋅ (B + C) ≈ A ⋅ B + A ⋅ C
    @test (B + C) ⋅ A ≈ B ⋅ A + C ⋅ A

    for r ∈ 0:4
        Aᵣ = As[r+1]
        for s ∈ 0:4
            Bₛ = Bs[s+1]
            for t ∈ 0:4
                Cₜ = Cs[t+1]
                if r + s ≤ t && r > 0 && s > 0
                    @test Aᵣ ⋅ (Bₛ ⋅ Cₜ) ≈ (Aᵣ ∧ Bₛ) ⋅ Cₜ
                elseif r + t ≤ s
                    @test Aᵣ ⋅ (Bₛ ⋅ Cₜ) ≈ (Aᵣ ⋅ Bₛ) ⋅ Cₜ
                end
            end
        end
        @test a ⋅ Aᵣ ≈ 0.5 * (a * Aᵣ + (-1) ^ (r + 1) * Aᵣ * a)
    end
end


text:
"""
    INT{N, AM} <: AbstractContinuousPost

Implementation of a reachability method for linear one-dimensional systems
interval arithmetic.

## Fields

- `δ`            -- step-size of the discretization
- `approx_model` -- (optional, default: `Forward`) approximation model;
                    see `Notes` below for possible options

## Notes

The type fields are:

- `N`  -- number type of the step-size
- `AM` -- approximation model

The default approximation model used in this algorithm is:

```julia
Forward(sih=:concrete, exp=:base, setops=:Interval)
```

In particular, the `setops=:Interval` flag specifies that intermediate computations
in the discretization are done using interval arithmetic. This allows for some
optimizations.

## References

This algorithm is essentially a non-decomposed version of the method in [[BFFPSV18]](@ref),
using intervals as set representation. For a general introduction we refer
to the dissertation [[LG09]](@ref).

Regarding the approximation model, by default we use an adaptation of the method
presented in [[FRE11]](@ref).

Interval arithmetic operations are performed using the `IntervalArithmetic.jl`
package. Hence, the results are guaranteed to comply to the IEE754 standard with
respect to the floating-point operations using intervals.
"""
struct INT{N, AM} <: AbstractContinuousPost
    δ::N
    approx_model::AM
end

# convenience constructor using symbols
function INT(; δ::N,
               approx_model::AM=Forward(sih=:concrete, exp=:base)) where {N, AM}
    return INT(δ, approx_model)
end

step_size(alg::INT) = alg.δ
numtype(::INT{N}) where {N} = N
setrep(::INT{N}) where {N} = Interval{N, IA.Interval{N}}
rsetrep(::INT{N}) where {N} = ReachSet{N, Interval{N, IA.Interval{N}}}

include("post.jl")
include("reach_homog.jl")
include("reach_inhomog.jl")


text:
using ExproniconLite
using Documenter
using Test
using Aqua
Aqua.test_all(ExproniconLite)

@testset "@test_expr" begin
    @test_expr quote
        x + 1
        $nothing
    end == quote
        x + 1
        $nothing
    end
end

@testset "printings" begin
    include("printing.jl")
end

@testset "analysis" begin
    include("analysis.jl")
end

@testset "transform" begin
    include("transform.jl")
end

@testset "match" begin
    nothing
end

@testset "codegen" begin
    include("codegen.jl")
end

# this feature is only available for 1.6+
@static if VERSION > v"1.6-" && Sys.isunix()
    @testset "expand" begin
        nothing
    end
end

doctest(ExproniconLite)


text:
"ω ∉ [0, 1]"
notunit(ω) = ω > 1.0 || ω < 0.0

function clip(x, eps = 1e-5)
  if x == 0.0
    eps
  elseif x == 1.0
    1 - eps
  else
    x
  end
end
@spec 0.0 < res < 1.0

"Bijective transformation from [0, 1] to the real line, T(x)=log(1/(1-x)-1)"
transform(x) = log(clip(x) / (1 - clip(x)))

"The inverse of the transformation above, T^(-1)(y)=1-1/(1+e^y)"
inv_transform(y) = 1 / (1 + exp(-y))

"Jacobian of the transformation above, J(x) = 1/x(1-x)"
jacobian(x) = inv_transform(x) * (1. -inv_transform(x))

"Bijective transformation from [0, 1] to the real line, T(x)=log(1/(1-x)-1)"
unbound(x) = log(clip(x)/(1-clip(x)))

"The inverse of the transformation above, T^(-1)(y)=1-1/(1+e^y)"
bound(y) = 1 / (1 + exp(-y))

"Jacobian of the transformation above, J(x) = 1/x(1-x)"
jac(x) = bound(x) * (1 - bound(x))


text:
#Tests des packages
include("../ALAS/include.jl")

using MPCCmod
using OutputRelaxationmod

using OptimizationProblems
using MathProgBase
using NLPModels
#8
sprob=names(OptimizationProblems)[91]
#sprob=:hs1
@time n=4
nlp=AbstractNLPModel
try
nlp = MathProgNLPModel(eval(sprob)(n))
catch
nlp = MathProgNLPModel(eval(sprob)())
end
nlp.meta.ncon==0 && @printf("Unconstrained problem? yes!\n")

#déclare le mpcc :
println("Create an MPCC")
@time exemple_nlp=MPCCmod.MPCC(nlp)
@printf("%i %i %i %i %i %i %i %i \n",nlp.counters.neval_obj,
        nlp.counters.neval_cons,nlp.counters.neval_grad,
        nlp.counters.neval_hess,exemple_nlp.G.counters.neval_cons,
        exemple_nlp.G.counters.neval_jac,exemple_nlp.H.counters.neval_cons
        ,exemple_nlp.H.counters.neval_jac)

println("Butterfly method:")
#résolution avec ALAS Butterfly
@time xb,fb,orb = MPCCsolve.solve(exemple_nlp)
@printf("%i %i %i %i %i %i %i %i \n",nlp.counters.neval_obj,
        nlp.counters.neval_cons,nlp.counters.neval_grad,
        nlp.counters.neval_hess,exemple_nlp.G.counters.neval_cons,
        exemple_nlp.G.counters.neval_jac,exemple_nlp.H.counters.neval_cons
        ,exemple_nlp.H.counters.neval_jac)


text:
import MathOptInterface
const MOI = MathOptInterface

"""
    strengthen(a, u0, v0, ::MOI.Nonnegatives)

Compute strenghtened coefficient for non-negative variable.

The strengthened coefficient is given by
```
α = min(-u0 * ⌊π⌋, a + v0 * ⌊π⌋)
```
where ``π = -a / (u0 + v0)``
"""
function strengthen(a::Float64, u0::Float64, v0::Float64, k::MOI.Nonnegatives)

    # Sanity checks
    MOI.dimension(k) == 1 || error("Trying to strengthen a $(typeof(k)) of dimension $(MOI.dimension(k)) (should be 1)")
    @assert u0 + v0 > 1e-4  "u0 + v0 = $(u0 + v0) (should be > 1e-4)"  # Otherwise cut is essentially a K* cut

    π = -a / (u0 + v0)
    return min(
        -u0 * floor(π),
        a + v0 * ceil(π)
    )
end

"""
    strengthen(a, u0, v0, ::MOI.Nonpositive)

Compute strenghtened coefficient for non-positive variable.

The strengthened coefficient is given by
```
α = min(-u0 * ⌊π⌋, a + v0 * ⌊π⌋)
```
where ``π = -a / (u0 + v0)``
"""
function strengthen(a::Float64, u0::Float64, v0::Float64, k::MOI.Nonpositives)

    # Sanity checks
    MOI.dimension(k) == 1 || error("Trying to strengthen a $(typeof(k)) of dimension $(MOI.dimension(k)) (should be 1)")
    @assert u0 + v0 > 1e-4  # Otherwise cut is essentially a K* cut

    π = -a / (u0 + v0)
    return max(
        -u0 * ceil(π),
        a + v0 * floor(π)
    )
end

# TODO: strengthening for conic components

max_stars_repo_name: mtanneau/CLaP

max_issues_repo_name: mtanneau/CLaP

max_forks_repo_name: mtanneau/CLaP


text:
using ..Pickle: Memo, PickleStack, HierarchicalTable, load, isdefer

using DataStructures
using Strided
using ZipFile

const DEFAULT_PROTO = 2

const MAGIC = BigInt(0x1950a86a20f9469cfc6c)
const TORCH_PROTOCOL = 1001

struct TorchPickler{PROTO} <: AbstractPickle
  memo::Memo
  stack::PickleStack
  mt::HierarchicalTable
  storage::StorageManager
end

function TorchPickler(proto=DEFAULT_PROTO, memo=Dict())
  st = StorageManager()
  mt = HierarchicalTable()

  # some corresponding methods
  mt["collections.OrderedDict"] = OrderedDict
  mt["torch._utils._rebuild_tensor_v2"] = (arg...) -> build_tensor(st, arg...)

  mt["__julia__.__torch__.rebuild_tensor"] = "torch._utils._rebuild_tensor_v2"
  mt["__julia__.OrderedDict"] = "collections.OrderedDict"
  mt["__julia__.__torch__.StorageType.Float32"] = "torch.FloatStorage"
  mt["__julia__.__torch__.StorageType.Float64"] = "torch.DoubleStorage"
  mt["__julia__.__torch__.StorageType.Float16"] = "torch.HalfStorage"
  mt["__julia__.__torch__.StorageType.UInt8"] = "torch.ByteStorage"
  mt["__julia__.__torch__.StorageType.Int8"] = "torch.CharStorage"
  mt["__julia__.__torch__.StorageType.Int16"] = "torch.ShortStorage"
  mt["__julia__.__torch__.StorageType.Int32"] = "torch.IntStorage"
  mt["__julia__.__torch__.StorageType.Int64"] = "torch.LongStorage"
  mt["__julia__.__torch__.StorageType.Bool"] = "torch.BoolStorage"

  # ingore state_dict version number
  mt["__build__.OrderedCollections.OrderedDict"] = (od, _meta) -> od
  mt["__build__.OrderedDict"] = (od, _meta...) -> od

  TorchPickler{proto}(Memo(memo), PickleStack(), mt, st)
end

protocol(::TorchPickler{P}) where {P} = P
isbinary(pklr::TorchPickler) = protocol(pklr) >= 1

"""
  THload(file::AbstractString)

load data that saved by `torch.save`. `torch.tensor`
will be load as `Array` or `Strided.StridedView`
dependent on the memory layout of that tensor.
"""
THload(file::AbstractString) = open(file) do io
  THload(TorchPickler(), io)
end

function THload(tp::TorchPickler, io)
  if peek(io) == 0x80
    return legacy_load(tp, io)
  elseif read(io, 4) == b"PK\x03\x04"
    z = ZipFile.Reader(io)
    if any(x->x.name=="constants.pkl", z.files)
      error("TorchScript archive not support.")
    end
    return zip_load(tp, z)
  else
    error("Unkown file format. Is this really a file from `torch.save`?")
  end
end

function get_record(zipfile, name)
  zipfile.files[findfirst(x->endswith(x.name, name), zipfile.files)]
end

function zip_load(tp::TorchPickler, zipfile)
  typeinfo = load(tp, get_record(zipfile, "data.pkl"))
  load_tensor_zip!(zipfile, tp.storage)
  return typeinfo
end

function unchecked_legacy_load(tp::TorchPickler, io)
  magic = load(tp, io)
  magic != MAGIC && error("Invalid magic number; corrupt file?")
  torch_protocol = load(tp, io)
  torch_protocol != TORCH_PROTOCOL && error("Invalid protocol version: $torch_protocol")

  _sys_info = load(tp, io)

  typeinfo = load(tp, io)
  tensor_key = load(tp, io)
  load_tensor!(io, tp.storage, tensor_key)
  return typeinfo
end

function legacy_load(tp::TorchPickler, io)
    typeinfo = unchecked_legacy_load(tp, io)
    @assert !isdefer(typeinfo)
    return typeinfo
end

function build_tensor(sm::StorageManager, fake_storage, offset, tsize, tstride, grad, _)
  @assert length(tsize) == length(tstride)
  @assert fake_storage.head == :persistent_load
  header, thtype, key, device, numel, = fake_storage.args[1]
  @assert header == "storage"
  dtype = thtype2dtype(thtype)
  jltype = dtype2jltype(dtype)
  tlength = prod(tsize)

  if haskey(sm, key)
    storage = sm[key][end]
  else
    storage = Array{jltype}(undef, numel)
    setindex!(sm, (dtype, numel, device, storage), key)
  end

  if (tlength == numel) && (isone(length(tsize)) || isempty(tsize) || isone(first(tstride)))
    tensor = reshape(storage, tsize) # f-contiguous
  else # otherwise use strided
    tensor = StridedView(storage, tsize, tstride, offset)
  end

  return tensor
end

function load_tensor!(io::IO, sm::StorageManager, tensor_key)
  for key in tensor_key
    type, numel, device, storage = sm[key]
    tsize = read(io, Int64)
    @assert tsize == numel
    tdata = read(io, tsize * bytewidth(type))
    storage .= reinterpret(dtype2jltype(type), tdata)
  end
end

function load_tensor_zip!(zipfile, sm::StorageManager)
  for (key, values) in pairs(sm)
    type, numel, device, storage = values
    zf = get_record(zipfile, "/$key")
    tdata = read(zf)
    storage .= reinterpret(dtype2jltype(type), tdata)
  end
end


text:
module TuringModels

using CategoricalArrays
using DataFrames

const project_root = string(pkgdir(TuringModels))::String

"""
    output_dir(name)

Return output directory for some model or script such as "africa".
This works around the OUTPUT macro from Franklin, because that macro doesn't work in Literate scripts.
"""
function output_dir(name)
    joinpath(project_root, "__site", "assets", "models", name, "code")
end

end # module


text:
abstract type AbstractStrategy end
struct BottomUp <: AbstractStrategy end
struct TopDown <: AbstractStrategy end

function initial_chart(tokens, grammar::AbstractGrammar, ::BottomUp)
    chart = Chart{rule_type(grammar), chart_key(grammar)}(length(tokens))
end

function initial_chart(tokens, grammar::AbstractGrammar{R}, ::TopDown) where {R}
    chart = Chart{R, chart_key(grammar)}(length(tokens))
    for arc in terminal_productions(grammar, tokens)
        push!(chart, PassiveArc(arc))
    end
    chart
end

struct Agenda{R, S <: SortedSet{ActiveArc{R}}}
    data::S
end

function Agenda{R}() where {R}
    s = SortedSet{ActiveArc{R}}(Base.Order.ReverseOrdering(Base.Order.By(arc -> (score(arc), objectid(arc)))))
    Agenda{R, typeof(s)}(s)
end
Base.push!(agenda::Agenda, arc::ActiveArc) = push!(agenda.data, arc)
Base.pop!(agenda::Agenda) = pop!(agenda.data)
Base.isempty(agenda::Agenda) = isempty(agenda.data)

function initial_agenda(tokens, grammar::AbstractGrammar{R}, ::BottomUp) where {R}
    agenda = Agenda{R}()
    for arc in terminal_productions(grammar, tokens)
        push!(agenda, ActiveArc(arc))
    end
    agenda
end

function initial_agenda(tokens, grammar::AbstractGrammar{R}, ::TopDown) where {R}
    agenda = Agenda{R}()
    for rule in productions(grammar)
        if lhs(rule) == start_symbol(grammar)
            push!(agenda, ActiveArc(Arc(0, 0, rule, [], score(rule))))
        end
    end
    agenda
end

struct PredictionCache{T}
    predictions::Set{Tuple{T, Int}}
end

PredictionCache{T}() where {T} = PredictionCache(Set{Tuple{T, Int}}())

"""
Returns `true` if the key was added, `false` otherwise.
"""
function maybe_push!(p::PredictionCache{T}, key::Tuple{T, Int}) where {T}
    # TODO: this can probably be done without two Set lookup operations
    if key in p.predictions
        return false
    else
        push!(p.predictions, key)
        return true
    end
end

struct ChartParser{R, G <: AbstractGrammar{R}, S <: AbstractStrategy, F}
    tokens::Vector{String}
    grammar::G
    strategy::S
    scoring_function::F
end

function ChartParser(tokens::AbstractVector{<:AbstractString},
            grammar::G,
            strategy::S=BottomUp(),
            scoring_function::F = arc -> 1.0) where {R, G <: AbstractGrammar{R}, S <: AbstractStrategy, F}
    ChartParser{R, G, S, F}(tokens, grammar, strategy, scoring_function)
end

struct ChartParserState{R, T}
    chart::Chart{R, T}
    agenda::Agenda{R}
    prediction_cache::PredictionCache{T}
end

function initial_state(parser::ChartParser{R}) where R
    chart = initial_chart(parser.tokens, parser.grammar, parser.strategy)
    agenda = initial_agenda(parser.tokens, parser.grammar, parser.strategy)
    prediction_cache = PredictionCache{chart_key(R)}()
    ChartParserState(chart, agenda, prediction_cache)
end

function Base.iterate(parser::ChartParser{R, T}, state=initial_state(parser)) where {R, T}
    while !isempty(state.agenda)
        candidate = pop!(state.agenda)
        if is_finished(candidate)
            arc = passive(candidate, parser.scoring_function)
            if score(arc) > 0
                update!(state, parser, arc)
                return (inner(arc), state)
            end
        else
            update!(state, parser, candidate)
        end
    end
    return nothing
end

Base.IteratorSize(::Type{<:ChartParser}) = Base.SizeUnknown()
Base.eltype(::Type{<:ChartParser{R}}) where {R} = Arc{R}

function is_complete(arc::Arc, parser::ChartParser)
    start(arc) == 0 && stop(arc) == length(parser.tokens) && head(arc) == start_symbol(parser.grammar)
end

is_complete(parser::ChartParser) = arc -> is_complete(arc, parser)

function update!(state::ChartParserState, parser::ChartParser, candidate::AbstractArc)
    push!(state.chart, candidate)
    for mate in mates(state.chart, candidate)
        push!(state.agenda, combine(candidate, mate))
    end
    predict!(state.agenda, state.chart, candidate, parser.grammar, state.prediction_cache, parser.strategy)
end

function predict!(agenda::Agenda, chart::Chart, candidate::ActiveArc,
                  grammar::AbstractGrammar{R}, prediction_cache::PredictionCache,
                  ::TopDown) where {R}
    is_new = maybe_push!(prediction_cache, (next_needed(candidate), stop(candidate)))
    if is_new
        for rule in productions(grammar)
            if lhs(rule) === next_needed(candidate)
                push!(agenda, ActiveArc(Arc(stop(candidate), stop(candidate), rule, Arc{R}[], score(rule))))
            end
        end
    end
end

function predict!(agenda::Agenda, chart::Chart, candidate::PassiveArc,
                  grammar::AbstractGrammar{R}, prediction_cache::PredictionCache,
                  ::TopDown) where {R}
    # Nothing to do here
end

function predict!(agenda::Agenda, chart::Chart, candidate::PassiveArc,
                  grammar::AbstractGrammar{R}, prediction_cache::PredictionCache,
                  ::BottomUp) where {R}
    is_new = maybe_push!(prediction_cache, (head(candidate), start(candidate)))
    if is_new
        for rule in productions(grammar)
            if first(rhs(rule)) === head(candidate)
                push!(agenda, ActiveArc(Arc(start(candidate), start(candidate), rule, Arc{R}[], score(rule))))
            end
        end
    end
end

function predict!(agenda::Agenda, chart::Chart, candidate::ActiveArc,
                  grammar::AbstractGrammar, prediction_cache::PredictionCache,
                  ::BottomUp)
    # Nothing to do here
end


text:
module Nucleotide_Essentials

import  CSV,
        DataFrames,
        Glob,
        StatsBase,
        DataFramesMeta,
        Plots,
        StatsPlots,
        Logging

export  readFastq,
        readFasta,
        FastqRecord,
        FastaRecord,
        writeFasta,
        FastqtoFasta,
        potential_mismatches,
        demultiplex_se,
        demultiplex_pe,
        reverse_complement,
        PlotQuality,
        FilterQuality_se,
        FilterQuality_pe

include("Reading/readFastq.jl")
include("Reading/readFasta.jl")
include("Writing/writeFasta.jl")
include("Writing/FastqtoFasta.jl")
include("Demultiplex/Demultiplex.jl")
include("QualityTesting/PlotQuality.jl")
include("QualityTesting/FilterQuality.jl")

end  # end module


text:
abstract type AbstractQueryOpt <: AbstractDict{Symbol,Any} end

Base.length(::T) where {T<:AbstractQueryOpt} = fieldcount(T)
Base.length(::Type{T}) where {T<:AbstractQueryOpt} = fieldcount(T)

function Base.iterate(aqo::T, state = 1) where {T<:AbstractQueryOpt}
  (state > length(T)) && return nothing
  (fieldname(T, state) => getfield(aqo, state), state + 1)
end

"""
    struct YahooOpt <: AbstractQueryOpt
      period1  # the start time
      period2  # the end time
      interval # "1d", "1wk" or "1mo"
      events   # currently only `:history` supported
    end

The Yahoo Finance HTTP API query object.

# Examples

```jl-repl
julia> t = Dates.now()
2020-08-09T01:38:04.735

julia> YahooOpt(period1 = t - Year(2), period2 = t)
YahooOpt{DateTime} with 4 entries:
  :period1  => 1533778685
  :period2  => 1596937085
  :interval => "1d"
  :events   => :history
```
"""
struct YahooOpt <: AbstractQueryOpt
  period1::DateTime
  period2::DateTime
  interval::String
  events::Symbol

  YahooOpt(; period1::DateTime = DateTime(1971, 2, 8),
          period2::DateTime = Dates.now(),
          interval::String  = "1d",
          events::Symbol    = :history) =
    new(period1, period2, interval, events)
end

function Base.iterate(opt::YahooOpt, state = 1)
  (state > length(YahooOpt)) && return nothing
  k = fieldname(YahooOpt, state)
  v = getfield(opt, state)
  v′ = (k ∈ (:period1, :period2)) ? round(Int, datetime2unix(v)) : v
  (k => v′, state + 1)
end

"""
    yahoo(symbol::AbstractString, opt::YahooOpt = YahooOpt())::TimeArray
    yahoo(symbol::Symbol, opt::YahooOpt = YahooOpt())::TimeArray

This is a wrapper for downloading historical stock prices from Yahoo Finance.

The yahoo method takes a stock name in the form of a string and returns a
`TimeSeries.TimeArray` corresponding to the Yahoo Finance ticker.
With no argument, the default historical time series is the S&P 500.

# Examples

```julia
AAPL = yahoo(:AAPL)
SPX = yahoo("^GSPC")
NQ = yahoo("^IXIC")
```

```jl-repl
julia> start = DateTime(2018, 1, 1)
2018-01-01T00:00:00

julia> yahoo(:AAPL, YahooOpt(period1 = start))
655×6 TimeArray{Float64,2,Date,Array{Float64,2}} 2018-01-02 to 2020-08-07
...
```

# References

  https://finance.yahoo.com

# See Also

- fred()  which accesses the St. Louis Federal Reserve financial and economic data sets.
- ons()   which is a wrapper to download financial and economic time series data from the Office for National Statistics (ONS).
"""
function yahoo(sym::AbstractString = "^GSPC", opt::YahooOpt = YahooOpt())
    host = rand(["query1", "query2"])
    url  = "https://$host.finance.yahoo.com/v7/finance/download/$sym"
    res  = HTTP.get(url, query = opt)
    @assert res.status == 200
    csv = CSV.File(res.body, missingstring = "null")
    sch = TimeSeries.Tables.schema(csv)
    TimeArray(csv, timestamp = first(sch.names)) |> cleanup_colname!
end

yahoo(s::Symbol, opt::YahooOpt = YahooOpt()) = yahoo(string(s), opt)

"""
  fred(data::String="CPIAUCNS")::TimeArray

The fred() method is a wrapper to download financial and economic time series data from the St. Louis Federal Reserve (FRED).

The fred() method takes a string argument that corresponds to a series code from the St. Louis Federal
Reserve (FRED) database. It returns the data in the TimeSeries.TimeArray data structure.  When no argument
is provided, the default data set is the Consumer Price Index for All Urban Consumers: All Items (CPIAUCNS).

# Examples

```julia
DGS = fred("DGS10")
CPI = fred()
```

# References

https://research.stlouisfed.org/fred2

# See Also

- yahoo() which is a wrapper to download financial time series for stocks from Yahoo Finance.
- ons()   which is a wrapper to download financial and economic time series data from the Office for National Statistics (ONS).
"""
function fred(data::AbstractString="CPIAUCNS")
    url = "http://research.stlouisfed.org/fred2/series/$data/downloaddata/$data.csv"
    res = HTTP.get(url)
    @assert res.status == 200
    csv = CSV.File(res.body)
    sch = TimeSeries.Tables.schema(csv)
    TimeArray(csv, timestamp = first(sch.names)) |> cleanup_colname!
end

"""
    ons(timeseries::String="L522", dataset::String="MM23")::TimeArray

The ons() method is a wrapper to download financial and economic time series data from the Office for National Statistics (ONS).

The ons() method takes two string arguments corresponding to a dataset and timeseries code from
the ONS database, to explore the database, you can use https://www.ons.gov.uk/timeseriestool.
It returns the data in the TimeSeries.TimeArray data structure with additional information in the meta field.
The data might include monthly values, quarterly averages and yearly averages, with column names
monthly, quarterly and yearly.
The timestamps are the first day of the period.  When no argument is provided, the default dataset
is the Consumer Price Index including housing costs (CPIH) which is the ONS’s headline measure of inflation.

# Examples

```julia
UK_RPI = ("CHAW","MM23")
UK_CPI = ("D7BT","MM23")
UK_CPIH = ("L522","MM23")
```

# References

https://www.ons.gov.uk/timeseriestool

# See Also

- fred() which accesses the St. Louis Federal Reserve financial and economic data sets.
- yahoo() which is a wrapper from downloading financial time series for stocks from Yahoo Finance.
"""
function ons(timeseries::AbstractString = "L522", dataset::AbstractString = "MM23")
    url = "https://api.ons.gov.uk/dataset/$dataset/timeseries/$timeseries/data"
    res = HTTP.get(url)
    @assert res.status == 200
    json = JSON3.read(HTTP.payload(res))
    ta = nothing
    if "months" in keys(json)
        data = json["months"]
        dates = [Date(x["year"]*" "*x["month"], dateformat"Y U") for x in data]
        values = [parse(Float64, x["value"]) for x in data]
        ta = TimeArray((timestamp = dates, monthly = values),timestamp=:timestamp)
    end
    if "quarters" in keys(json)
        data = json["quarters"]
        dates = [Date(x["year"]*" "*string(parse(Int, x["quarter"][2])*3-2),dateformat"Y m") for x in data]
        values = [parse(Float64,x["value"]) for x in data]
        tb = TimeArray((timestamp = dates, quarterly = values), timestamp = :timestamp)
        ta === nothing ? ta = tb : ta = merge(ta, tb; method = :outer)
    end
    if "years" in keys(json)
        data = json["years"]
        dates = [Date(x["year"]) for x in data]
        values = [parse(Float64, x["value"]) for x in data]
        tb = TimeArray((timestamp = dates, yearly = values), timestamp = :timestamp)
        ta === nothing ? ta = tb : ta = merge(ta, tb; method = :outer)
    end
    TimeArray(ta, meta = json["description"])
end


text:
"""
Module Email

Provides functions for operations in an Email client.

# Functions

###### Open inbox
- [`inbox`](@ref)

To see a description of a function type `?<functionname>`.
"""
module Email

import DefaultApplication
import ..JustSayIt: @voiceargs

## CONSTANTS

const EMAIL_ENGINE = "https://mail.google.com"

## Functions

@doc """
    email `inbox`

Manage e-mails, performing one of the following actions:
    - `inbox`
"""
email
@enum Action inbox
@voiceargs action=>(valid_input_auto=true, use_max_accuracy=true) function email(action::Action)
    if (action == inbox) open_inbox()
    else @info "unknown action" #NOTE: this should never happen.
    end
end

@doc "Search in internet: start search engine and automatically change to type command waiting for search keywords to be spoken."
function open_inbox()
    DefaultApplication.open(EMAIL_ENGINE)
end

end # module Email


text:
#Write function that reverses a list, preferably in place.

list = [1,23,3,4, "A","B"]

println(reverse!(list))

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


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

JLLWrappers.@generate_main_file_header("MagicVLSI")
JLLWrappers.@generate_main_file("MagicVLSI", UUID("8553b089-32f9-510e-b5d7-7da551ec2160"))
end  # module MagicVLSI_jll


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using BenchmarkTools

function get_input()
    parse.(Int, eachline("day09input.txt"))
end

function valid_next(window::AbstractSet{Int}, next::Int)::Bool
    for num ∈ window
        if (next - num ∈ window) && num*2 != next
            return true
        end
    end
    return false
end

function part1(nums::Array{Int}, window_size::Int)
    n = length(nums)
    if window_size >= n
        throw(ErrorException("Window size must be smaller than the input"))
    end
    window_set = Set(@view nums[1:window_size])
    for i ∈ (window_size+1):n
        val = nums[i]
        if !valid_next(window_set, val)
            return val
        end
        delete!(window_set, nums[i-window_size])
        push!(window_set, val)
    end
end

function part2(nums::Array{Int}, target::Int)::Int
    n = length(nums)
    left, right = n-1, n
    runningsum = nums[left] + nums[right]
    while runningsum != target
        if runningsum < target
            left -= 1
            runningsum += nums[left]
        else
            runningsum -= nums[right]
            right -= 1
        end
    end
    sum(extrema(@view nums[left:right]))
end

println("Parsing time:")
nums = @btime get_input()
println("Part 1:")
answer1 = @btime part1(nums, 25)
println(answer1)
println("Part 2:")
answer2 = @btime part2(nums, answer1)
println(answer2)

function solve_all()
    nums = get_input()
    answer1 = part1(nums, 25)
    answer2 = part2(nums, answer1)
end


text:
function equals(obj::GradientField, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_field(::Type{GradientField}, arg0::jint)
    return jcall(GradientField, "getField", GradientField, (jint,), arg0)
end

function get_one(obj::GradientField)
    return jcall(obj, "getOne", Gradient, ())
end

function get_runtime_class(obj::GradientField)
    return jcall(obj, "getRuntimeClass", Class, ())
end

function get_zero(obj::GradientField)
    return jcall(obj, "getZero", Gradient, ())
end

function hash_code(obj::GradientField)
    return jcall(obj, "hashCode", jint, ())
end


text:
using Test

for submodule in [
                  "Grids",
                  "Fields",
                  "FFTs",
                 ]

  println("Testing $submodule")
  t = @elapsed include(joinpath(submodule,"runtests.jl"))
  println("Completed tests for $submodule, $(round(Int, t)) seconds elapsed")
end


text:
using DataFrames, CSV

const postComBat = CSV.read("../transcriptomic_data/TPM_RLE_postComBat.csv", DataFrame);

function get_tpm(gene::String, cell_line::String)::Float64
    return postComBat[postComBat[!,Symbol("Description")] .== gene,Symbol(cell_line)][1]
end

function incorporate_ErbB1(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_EGFR] * get_tpm("EGFR", cell_line)
end

function incorporate_ErbB2(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_ERBB2] * get_tpm("ERBB2", cell_line)
end

function incorporate_ErbB3(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_ERBB3] * get_tpm("ERBB3", cell_line)
end

function incorporate_Erb4(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_ERBB4] * get_tpm("ERBB4", cell_line)
end

function incorporate_Grb2(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_GRB2] * get_tpm("GRB2", cell_line)
end

function incorporate_Shc(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_SHC1] * get_tpm("SHC1", cell_line) +
            p[C.w_SHC2] * get_tpm("SHC2", cell_line) +
            p[C.w_SHC3] * get_tpm("SHC3", cell_line) +
            p[C.w_SHC4] * get_tpm("SHC4", cell_line))
end

function incorporate_RasGAP(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_RASA1] * get_tpm("RASA1", cell_line) +
            p[C.w_RASA2] * get_tpm("RASA2", cell_line) +
            p[C.w_RASA3] * get_tpm("RASA3", cell_line))
end

function incorporate_PI3K(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_PIK3CA] * get_tpm("PIK3CA", cell_line) +
            p[C.w_PIK3CB] * get_tpm("PIK3CB", cell_line) +
            p[C.w_PIK3CD] * get_tpm("PIK3CD", cell_line) +
            p[C.w_PIK3CG] * get_tpm("PIK3CG", cell_line))
end

function incorporate_PTEN(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_PTEN] * get_tpm("PTEN", cell_line)
end

function incorporate_SOS(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_SOS1] * get_tpm("SOS1", cell_line) +
            p[C.w_SOS2] * get_tpm("SOS2", cell_line))
end

function incorporate_Gab1(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_GAB1] * get_tpm("GAB1", cell_line)
end

function incorporate_RasGDP(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_HRAS] * get_tpm("HRAS", cell_line) +
            p[C.w_KRAS] * get_tpm("KRAS", cell_line) +
            p[C.w_NRAS] * get_tpm("NRAS", cell_line))
end

function incorporate_Raf(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_ARAF] * get_tpm("ARAF", cell_line) +
            p[C.w_BRAF] * get_tpm("BRAF", cell_line) +
            p[C.w_RAF1] * get_tpm("RAF1", cell_line))
end

function incorporate_MEK(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_MAP2K1] * get_tpm("MAP2K1", cell_line) +
            p[C.w_MAP2K2] * get_tpm("MAP2K2", cell_line))
end

function incorporate_ERK(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_MAPK1] * get_tpm("MAPK1", cell_line) +
            p[C.w_MAPK3] * get_tpm("MAPK3", cell_line))
end

function incorporate_Akt(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_AKT1] * get_tpm("AKT1", cell_line) +
            p[C.w_AKT2] * get_tpm("AKT2", cell_line))
end

function incorporate_PTP1B(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_PTPN1] * get_tpm("PTPN1", cell_line)
end

function incorporate_GSK3b(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_GSK3B] * get_tpm("GSK3B", cell_line)
end

function incorporate_DUSP(p::Vector{Float64}, cell_line::String)::Float64
    return (p[C.w_DUSP5] * get_tpm("DUSP5", cell_line) +
            p[C.w_DUSP6] * get_tpm("DUSP6", cell_line) +
            p[C.w_DUSP7] * get_tpm("DUSP7", cell_line))
end

function incorporate_cMyc(p::Vector{Float64}, cell_line::String)::Float64
    return p[C.w_MYC] * get_tpm("MYC", cell_line)
end

function individualize!(
        p::Vector{Float64},
        u0::Vector{Float64},
        cell_line::String)::Tuple{Vector{Float64},Vector{Float64}}
    # as initial concentration of model species
    u0[V.ErbB1] *= incorporate_ErbB1(p, cell_line)
    u0[V.ErbB3] *= incorporate_ErbB3(p, cell_line)
    u0[V.ErbB4] *= incorporate_Erb4(p, cell_line)
    u0[V.ErbB2] *= incorporate_ErbB2(p, cell_line)
    u0[V.Grb2] *= incorporate_Grb2(p, cell_line)
    u0[V.Shc] *= incorporate_Shc(p, cell_line)
    u0[V.RasGAP] *= incorporate_RasGAP(p, cell_line)
    u0[V.PI3K] *= incorporate_PI3K(p, cell_line)
    u0[V.PTP1B] *= incorporate_PTP1B(p, cell_line)
    u0[V.SOS] *= incorporate_SOS(p, cell_line)
    u0[V.Gab1] *= incorporate_Gab1(p, cell_line)
    u0[V.RasGDP] *= incorporate_RasGDP(p, cell_line)
    u0[V.Raf] *= incorporate_Raf(p, cell_line)
    u0[V.MEK] *= incorporate_MEK(p, cell_line)
    u0[V.ERK] *= incorporate_ERK(p, cell_line)
    u0[V.PTEN] *= incorporate_PTEN(p, cell_line)
    u0[V.Akt] *= incorporate_Akt(p, cell_line)
    u0[V.GSK3b] *= incorporate_GSK3b(p, cell_line)

    # as transcription rates
    p[C.V291] *= incorporate_DUSP(p, cell_line)
    p[C.V310] *= incorporate_cMyc(p, cell_line)

    return p, u0
end

function scale_using_mcf7!(
        p::Vector{Float64},
        u0::Vector{Float64},
        cell_line_2::String)::Tuple{Vector{Float64},Vector{Float64}}
    # as initial concentration of model species
    u0[V.ErbB1] *= incorporate_ErbB1(p, cell_line_2) / incorporate_ErbB1(p, "MCF7_BREAST")
    u0[V.ErbB3] *= incorporate_ErbB3(p, cell_line_2) / incorporate_ErbB3(p, "MCF7_BREAST")
    u0[V.ErbB4] *= incorporate_Erb4(p, cell_line_2) / incorporate_Erb4(p, "MCF7_BREAST")
    u0[V.ErbB2] *= incorporate_ErbB2(p, cell_line_2) / incorporate_ErbB2(p, "MCF7_BREAST")
    u0[V.Grb2] *= incorporate_Grb2(p, cell_line_2) / incorporate_Grb2(p, "MCF7_BREAST")
    u0[V.Shc] *= incorporate_Shc(p, cell_line_2) / incorporate_Shc(p, "MCF7_BREAST")
    u0[V.RasGAP] *= incorporate_RasGAP(p, cell_line_2) / incorporate_RasGAP(p, "MCF7_BREAST")
    u0[V.PI3K] *= incorporate_PI3K(p, cell_line_2) / incorporate_PI3K(p, "MCF7_BREAST")
    u0[V.PTP1B] *= incorporate_PTP1B(p, cell_line_2) / incorporate_PTP1B(p, "MCF7_BREAST")
    u0[V.SOS] *= incorporate_SOS(p, cell_line_2) / incorporate_SOS(p, "MCF7_BREAST")
    u0[V.Gab1] *= incorporate_Gab1(p, cell_line_2) / incorporate_Gab1(p, "MCF7_BREAST")
    u0[V.RasGDP] *= incorporate_RasGDP(p, cell_line_2) / incorporate_RasGDP(p, "MCF7_BREAST")
    u0[V.Raf] *= incorporate_Raf(p, cell_line_2) / incorporate_Raf(p, "MCF7_BREAST")
    u0[V.MEK] *= incorporate_MEK(p, cell_line_2) / incorporate_MEK(p, "MCF7_BREAST")
    u0[V.ERK] *= incorporate_ERK(p, cell_line_2) / incorporate_ERK(p, "MCF7_BREAST")
    u0[V.PTEN] *= incorporate_PTEN(p, cell_line_2) / incorporate_PTEN(p, "MCF7_BREAST")
    u0[V.Akt] *= incorporate_Akt(p, cell_line_2) / incorporate_Akt(p, "MCF7_BREAST")
    u0[V.GSK3b] *= incorporate_GSK3b(p, cell_line_2) / incorporate_GSK3b(p, "MCF7_BREAST")

    # as transcription rates
    p[C.V291] *= incorporate_DUSP(p, cell_line_2) / incorporate_DUSP(p, "MCF7_BREAST")
    p[C.V310] *= incorporate_cMyc(p, cell_line_2) / incorporate_cMyc(p, "MCF7_BREAST")

    return p, u0
end