text:
########################### Thermal Generation Models ######################################
abstract type AbstractThermalFormulation <: AbstractDeviceFormulation end
abstract type AbstractThermalDispatchFormulation <: AbstractThermalFormulation end
abstract type AbstractThermalUnitCommitment <: AbstractThermalFormulation end
struct ThermalBasicUnitCommitment <: AbstractThermalUnitCommitment end
struct ThermalStandardUnitCommitment <: AbstractThermalUnitCommitment end
struct ThermalDispatch <: AbstractThermalDispatchFormulation end
struct ThermalRampLimited <: AbstractThermalDispatchFormulation end
struct ThermalDispatchNoMin <: AbstractThermalDispatchFormulation end
struct ThermalMultiStartUnitCommitment <: AbstractThermalUnitCommitment end

########################### Active Dispatch Variables ######################################
"""
This function add the variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: ActivePowerVariable, U <: PSY.ThermalGen}
    if get_warm_start(psi_container.settings)
        initial_value_func = d -> PSY.get_active_power(d)
    else
        initial_value_func = nothing
    end
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_active,
        initial_value_func = initial_value_func,
        lb_value_func = x -> PSY.get_active_power_limits(x).min,
        ub_value_func = x -> PSY.get_active_power_limits(x).max,
    )
end

function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: ActivePowerVariable, U <: PSY.ThermalMultiStart}
    if get_warm_start(psi_container.settings)
        initial_value_func = d -> PSY.get_active_power(d)
    else
        initial_value_func = nothing
    end
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_active,
        initial_value_func = initial_value_func,
        lb_value_func = x -> 0,
        ub_value_func = x -> PSY.get_active_power_limits(x).max,
    )
end

"""
This function add the variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: ReactivePowerVariable, U <: PSY.ThermalGen}
    if get_warm_start(psi_container.settings)
        initial_value_func = d -> PSY.get_reactive_power(d)
    else
        initial_value_func = nothing
    end
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_reactive,
        initial_value_func = initial_value_func,
        lb_value_func = x -> PSY.get_reactive_power_limits(x).min,
        ub_value_func = x -> PSY.get_reactive_power_limits(x).max,
    )
end

"""
This function add the variables for power generation commitment to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: OnVariable, U <: PSY.ThermalGen}
    if get_warm_start(psi_container.settings)
        initial_value_func = x -> (PSY.get_active_power(x) > 0 ? 1.0 : 0.0)
    else
        initial_value_func = nothing
    end
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        initial_value_func = initial_value_func,
        binary = true,
    )
end

function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: StopVariable, U <: PSY.ThermalGen}
    AddVariableSpec(; variable_name = make_variable_name(T, U), binary = true)
end

function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: StartVariable, U <: PSY.ThermalGen}
    AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = true, # TODO: This variable could be relaxed. Consider this as a different formulation
        lb_value_func = x -> 0.0,
        ub_value_func = x -> 1.0,
    )
end

function commitment_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.ThermalMultiStart},
)
    time_steps = model_time_steps(psi_container)
    if get_warm_start(psi_container.settings)
        initial_value = d -> (PSY.get_active_power(d) > 0 ? 1.0 : 0.0)
    else
        initial_value = nothing
    end

    add_variable!(
        psi_container,
        devices,
        make_variable_name(OnVariable, PSY.ThermalMultiStart),
        true,
    )
    varstatus = get_variable(psi_container, OnVariable, PSY.ThermalMultiStart)
    for t in time_steps, d in devices
        name = PSY.get_name(d)
        bus_number = PSY.get_number(PSY.get_bus(d))
        add_to_expression!(
            get_expression(psi_container, :nodal_balance_active),
            bus_number,
            t,
            varstatus[name, t],
            PSY.get_active_power_limits(d).min,
        )
    end

    var_names = (
        make_variable_name(StartVariable, PSY.ThermalMultiStart),
        make_variable_name(StopVariable, PSY.ThermalMultiStart),
    )
    for v in var_names
        add_variable!(psi_container, devices, v, true)
    end

    return
end

function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {
    T <: Union{ColdStartVariable, WarmStartVariable, HotStartVariable},
    U <: PSY.ThermalMultiStart,
}
    return AddVariableSpec(; variable_name = make_variable_name(T, U), binary = true)
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{<:VariableType},
    ::Type{T},
    ::Type{<:AbstractThermalFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::SemiContinuousFF,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ThermalGen}
    return DeviceRangeConstraintSpec()
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{<:VariableType},
    ::Type{T},
    ::Type{<:ThermalDispatchNoMin},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::SemiContinuousFF,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ThermalGen}
    return DeviceRangeConstraintSpec()
end

"""
This function adds the active power limits of generators when there are no CommitmentVariables
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractThermalDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ThermalGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function adds the active power limits of generators when there are CommitmentVariables
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractThermalFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ThermalGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            bin_variable_names = [make_variable_name(OnVariable, T)],
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function adds the active power limits of generators when there are
    no CommitmentVariables
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:ThermalDispatchNoMin},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ThermalGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            limits_func = x -> (min = 0.0, max = PSY.get_active_power_limits(x).max),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
        custom_psi_container_func = custom_active_power_constraints!,
    )
end

function custom_active_power_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{<:ThermalDispatchNoMin},
) where {T <: PSY.ThermalGen}
    var_key = make_variable_name(ActivePowerVariable, T)
    variable = get_variable(psi_container, var_key)
    # If the variable was a lower bound != 0, not removing the LB can cause infeasibilities
    for v in variable
        if JuMP.has_lower_bound(v)
            JuMP.set_lower_bound(v, 0.0)
        end
    end
end

"""
This function adds the active power limits of generators. Constraint (17) & (18) from PGLIB
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:ThermalMultiStartUnitCommitment},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ThermalMultiStart}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            limits_func = x -> (
                min = 0.0,
                max = PSY.get_active_power_limits(x).max -
                      PSY.get_active_power_limits(x).min,
            ),
            bin_variable_names = [
                make_variable_name(OnVariable, T),
                make_variable_name(StartVariable, T),
                make_variable_name(StopVariable, T),
            ],
            constraint_func = device_multistart_range!,
            constraint_struct = DeviceMultiStartRangeConstraintsInfo,
            lag_limits_func = PSY.get_power_trajectory,
        ),
    )
end

function _get_data_for_range_ic(
    initial_conditions_power::Vector{InitialCondition},
    initial_conditions_status::Vector{InitialCondition},
)
    lenght_devices_power = length(initial_conditions_power)
    lenght_devices_status = length(initial_conditions_status)
    @assert lenght_devices_power == lenght_devices_status
    ini_conds = Matrix{InitialCondition}(undef, lenght_devices_power, 2)
    idx = 0
    for (ix, ic) in enumerate(initial_conditions_power)
        g = ic.device
        @assert g == initial_conditions_status[ix].device
        idx += 1
        ini_conds[idx, 1] = ic
        ini_conds[idx, 2] = initial_conditions_status[ix]
    end
    return ini_conds
end

"""
This function adds range constraint for the first time period. Constraint (10) from PGLIB formulation
"""
function initial_range_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.ThermalMultiStart},
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {S <: PM.AbstractPowerModel}
    time_steps = model_time_steps(psi_container)
    resolution = model_resolution(psi_container)
    key_power = ICKey(DevicePower, PSY.ThermalMultiStart)
    key_status = ICKey(DeviceStatus, PSY.ThermalMultiStart)
    initial_conditions_power = get_initial_conditions(psi_container, key_power)
    initial_conditions_status = get_initial_conditions(psi_container, key_status)
    ini_conds = _get_data_for_range_ic(initial_conditions_power, initial_conditions_status)

    constraint_data = Vector{DeviceMultiStartRangeConstraintsInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        limits = PSY.get_active_power_limits(d)
        name = PSY.get_name(d)
        @assert name == PSY.get_name(ini_conds[ix, 1].device)
        lag_ramp_limits = PSY.get_power_trajectory(d)
        range_data = DeviceMultiStartRangeConstraintsInfo(name, limits, lag_ramp_limits)
        add_device_services!(range_data, d, model)
        constraint_data[ix] = range_data
    end

    if !isempty(ini_conds)
        device_multistart_range_ic!(
            psi_container,
            constraint_data,
            ini_conds,
            make_constraint_name(ACTIVE_RANGE_IC, PSY.ThermalMultiStart),
            make_variable_name(StopVariable, PSY.ThermalMultiStart),
        )
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end
    return
end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractThermalDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ThermalGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ReactivePowerVariable,
                T,
            ),
            variable_name = make_variable_name(ReactivePowerVariable, T),
            limits_func = x -> PSY.get_reactive_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function adds the reactive power limits of generators when there CommitmentVariables
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractThermalFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.ThermalGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ReactivePowerVariable,
                T,
            ),
            variable_name = make_variable_name(ReactivePowerVariable, T),
            bin_variable_names = [make_variable_name(OnVariable, T)],
            limits_func = x -> PSY.get_reactive_power_limits(x),
            constraint_func = device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

### Constraints for Thermal Generation without commitment variables ####
"""
This function adds the Commitment Status constraint when there are CommitmentVariables
"""
function commitment_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, D <: AbstractThermalFormulation, S <: PM.AbstractPowerModel}
    device_commitment!(
        psi_container,
        get_initial_conditions(psi_container, ICKey(DeviceStatus, T)),
        make_constraint_name(COMMITMENT, T),
        (
            make_variable_name(StartVariable, T),
            make_variable_name(StopVariable, T),
            make_variable_name(OnVariable, T),
        ),
    )
    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{D},
) where {T <: PSY.ThermalGen, D <: AbstractThermalUnitCommitment}
    status_init(psi_container, devices)
    output_init(psi_container, devices)
    duration_init(psi_container, devices)
    return
end

function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{ThermalBasicUnitCommitment},
) where {T <: PSY.ThermalGen}
    status_init(psi_container, devices)
    output_init(psi_container, devices)
    return
end

function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{D},
) where {T <: PSY.ThermalGen, D <: AbstractThermalDispatchFormulation}
    output_init(psi_container, devices)
    return
end
########################### Ramp/Rate of Change Constraints ################################
"""
This function gets the data for the generators for ramping constraints of thermal generators
"""
function _get_data_for_rocc(
    psi_container::PSIContainer,
    ::Type{T},
) where {T <: PSY.ThermalGen}
    resolution = model_resolution(psi_container)
    if resolution > Dates.Minute(1)
        minutes_per_period = Dates.value(Dates.Minute(resolution))
    else
        @warn("Not all formulations support under 1-minute resolutions. Exercise caution.")
        minutes_per_period = Dates.value(Dates.Second(resolution)) / 60
    end

    initial_conditions_power = get_initial_conditions(psi_container, DevicePower, T)
    lenght_devices_power = length(initial_conditions_power)
    data = Vector{DeviceRampConstraintInfo}(undef, lenght_devices_power)
    idx = 0
    for (ix, ic) in enumerate(initial_conditions_power)
        g = ic.device
        name = PSY.get_name(g)
        non_binding_up = false
        non_binding_down = false
        ramp_limits = PSY.get_ramp_limits(g)
        if !isnothing(ramp_limits)
            p_lims = PSY.get_active_power_limits(g)
            max_rate = abs(p_lims.min - p_lims.max) / minutes_per_period
            if (ramp_limits.up >= max_rate) & (ramp_limits.down >= max_rate)
                @debug "Generator $(name) has a nonbinding ramp limits. Constraints Skipped"
                continue
            else
                idx += 1
            end
            ramp = (
                up = ramp_limits.up * minutes_per_period,
                down = ramp_limits.down * minutes_per_period,
            )
            data[idx] = DeviceRampConstraintInfo(name, p_lims, ic, ramp)
        end
    end
    if idx < lenght_devices_power
        deleteat!(data, (idx + 1):lenght_devices_power)
    end
    return data
end

"""
This function adds the ramping limits of generators when there are CommitmentVariables
"""
function ramp_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, D <: AbstractThermalFormulation, S <: PM.AbstractPowerModel}
    time_steps = model_time_steps(psi_container)
    data = _get_data_for_rocc(psi_container, T)
    if !isempty(data)
        # Here goes the reactive power ramp limits when versions for AC and DC are added
        device_mixedinteger_rateofchange!(
            psi_container,
            data,
            make_constraint_name(RAMP, T),
            (
                make_variable_name(ActivePowerVariable, T),
                make_variable_name(StartVariable, T),
                make_variable_name(StopVariable, T),
            ),
        )
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end
    return
end

function ramp_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    time_steps = model_time_steps(psi_container)
    data = _get_data_for_rocc(psi_container, T)
    if !isempty(data)
        # Here goes the reactive power ramp limits when versions for AC and DC are added
        device_linear_rateofchange!(
            psi_container,
            data,
            make_constraint_name(RAMP, T),
            make_variable_name(ActivePowerVariable, T),
        )
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end
    return
end

function ramp_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.ThermalMultiStart},
    model::DeviceModel{T, ThermalMultiStartUnitCommitment},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalMultiStart, S <: PM.AbstractPowerModel}
    time_steps = model_time_steps(psi_container)
    data = _get_data_for_rocc(psi_container, T)

    # TODO: Refactor this to a cleaner format that doesn't require passing the device and rate_data this way
    for r in data
        add_device_services!(r, r.ic_status.device, model)
    end
    if !isempty(data)
        device_multistart_rateofchange!(
            psi_container,
            constaint_data,
            make_constraint_name(RAMP, PSY.ThermalMultiStart),
            make_variable_name(ActivePowerVariable, PSY.ThermalMultiStart),
        )
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end
    return
end

########################### start up trajectory constraints ######################################

@doc raw"""
    turbine_temperature(psi_container::PSIContainer,
                            startup_data::Vector{DeviceStartUpConstraintInfo},
                            cons_name::Symbol,
                            var_stop::Symbol,
                            var_starts::Tuple{Symbol, Symbol})

Constructs contraints for different types of starts based on generator down-time

# Equations
for t in time_limits[s+1]:T

``` var_starts[name, s, t] <= sum( var_stop[name, t-i] for i in time_limits[s]:(time_limits[s+1]-1)  ```

# LaTeX

``  δ^{s}(t)  \leq \sum_{i=TS^{s}_{g}}^{TS^{s+1}_{g}} x^{stop}(t-i) ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* rate_data::Vector{DeviceStartUpConstraintInfo},
* cons_name::Symbol : name of the constraint
* var_stop::Symbol : name of the stop variable
* var_starts::Tuple{Symbol, Symbol} : the names of the different start variables
"""
function turbine_temperature(
    psi_container::PSIContainer,
    startup_data::Vector{DeviceStartUpConstraintInfo},
    cons_name::Symbol,
    var_stop::Symbol,
    var_starts::Tuple{Symbol, Symbol},
)
    time_steps = model_time_steps(psi_container)
    start_vars = [
        get_variable(psi_container, var_starts[1]),
        get_variable(psi_container, var_starts[2]),
    ]
    varstop = get_variable(psi_container, var_stop)

    hot_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "hot")
    warm_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "warm")

    names = [get_component_name(st) for st in startup_data]

    con = [
        add_cons_container!(psi_container, hot_name, names, time_steps; sparse = true),
        add_cons_container!(psi_container, warm_name, names, time_steps; sparse = true),
    ]

    # constraint (15)
    for t in time_steps, st in startup_data
        for ix in 1:(st.startup_types - 1)
            if t >= st.time_limits[ix + 1]
                name = get_component_name(st)
                con[ix][name, t] = JuMP.@constraint(
                    psi_container.JuMPmodel,
                    start_vars[ix][name, t] <= sum(
                        varstop[name, t - i]
                        for i in st.time_limits[ix]:(st.time_limits[ix + 1] - 1)
                    )
                )
            end
        end
    end
    return
end

@doc raw"""
    device_start_type_constraint(psi_container::PSIContainer,
                            data::Vector{DeviceStartTypesConstraintInfo},
                            cons_name::Symbol,
                            var_start::Symbol,
                            var_names::Tuple{Symbol, Symbol, Symbol},)

Constructs contraints that restricts devices to one type of start at a time

# Equations

``` sum(var_starts[name, s, t] for s in starts) = var_start[name, t]  ```

# LaTeX

``  \sum^{S_g}_{s=1} δ^{s}(t)  \eq  x^{start}(t) ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* data::Vector{DeviceStartTypesConstraintInfo},
* cons_name::Symbol : name of the constraint
* var_start::Symbol : name of the startup variable
* var_starts::Tuple{Symbol, Symbol} : the names of the different start variables
"""
function device_start_type_constraint(
    psi_container::PSIContainer,
    data::Vector{DeviceStartTypesConstraintInfo},
    cons_name::Symbol,
    var_start::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
)
    time_steps = model_time_steps(psi_container)
    varstart = get_variable(psi_container, var_start)
    start_vars = [
        get_variable(psi_container, var_names[1]),
        get_variable(psi_container, var_names[2]),
        get_variable(psi_container, var_names[3]),
    ]

    set_name = [get_component_name(d) for d in data]
    con = add_cons_container!(psi_container, cons_name, set_name, time_steps)

    for t in time_steps, d in data
        # constraint (16)
        name = get_component_name(d)
        con[name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            varstart[name, t] == sum(start_vars[ix][name, t] for ix in 1:(d.startup_types))
        )
    end
    return
end

@doc raw"""
    device_startup_initial_condition(psi_container::PSIContainer,
                            data::Vector{DeviceStartUpConstraintInfo},
                            initial_conditions::Vector{InitialCondition},
                            cons_name::Symbol,
                            var_names::Tuple{Symbol, Symbol},
                            bin_name::Symbol,)

Constructs contraints that restricts devices to one type of start at a time

# Equations
ub:
``` (time_limits[st+1]-1)*δ^{s}(t) + (1 - δ^{s}(t)) * M_VALUE >= sum(1-varbin[name, i]) for i in 1:t) + initial_condition_offtime  ```
lb:
``` (time_limits[st]-1)*δ^{s}(t) =< sum(1-varbin[name, i]) for i in 1:t) + initial_condition_offtime  ```

# LaTeX

`` TS^{s+1}_{g} δ^{s}(t) + (1-δ^{s}(t)) M_VALUE   \geq  \sum^{t}_{i=1} x^{status}(i)  +  DT_{g}^{0}  \forall t in \{1, \ldots,  TS^{s+1}_{g}``

`` TS^{s}_{g} δ^{s}(t) \leq  \sum^{t}_{i=1} x^{status}(i)  +  DT_{g}^{0}  \forall t in \{1, \ldots,  TS^{s+1}_{g}``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* data::Vector{DeviceStartTypesConstraintInfo},
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol} : the names of the different start variables
* bin_name::Symbol : name of the status variable
"""
function device_startup_initial_condition(
    psi_container::PSIContainer,
    data::Vector{DeviceStartUpConstraintInfo},
    initial_conditions::Vector{InitialCondition},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol},
    bin_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    T = length(time_steps)

    set_name = [device_name(ic) for ic in initial_conditions]
    up_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    down_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")
    varbin = get_variable(psi_container, bin_name)
    varstarts = [
        get_variable(psi_container, var_names[1]),
        get_variable(psi_container, var_names[2]),
    ]

    con_ub = add_cons_container!(
        psi_container,
        up_name,
        set_name,
        time_steps,
        1:(MAX_START_STAGES - 1);
        sparse = true,
    )
    con_lb = add_cons_container!(
        psi_container,
        down_name,
        set_name,
        time_steps,
        1:(MAX_START_STAGES - 1);
        sparse = true,
    )

    for t in time_steps, (ix, d) in enumerate(data)
        name = get_component_name(d)
        ic = initial_conditions[ix]
        for st in 1:(d.startup_types - 1)
            var = varstarts[st]
            if t < (d.time_limits[st + 1] - 1)
                con_ub[name, t, st] = JuMP.@constraint(
                    psi_container.JuMPmodel,
                    (d.time_limits[st + 1] - 1) * var[name, t] +
                    (1 - var[name, t]) * M_VALUE >=
                    sum((1 - varbin[name, i]) for i in 1:t) + ic.value
                )
                con_lb[name, t, st] = JuMP.@constraint(
                    psi_container.JuMPmodel,
                    d.time_limits[st] * var[name, t] <=
                    sum((1 - varbin[name, i]) for i in 1:t) + ic.value
                )
            end
        end
    end
    return
end

"""
This function creates the contraints for different types of starts based on generator down-time
"""
function startup_time_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.ThermalMultiStart},
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {S <: PM.AbstractPowerModel}
    time_steps = model_time_steps(psi_container)
    constraint_data = Vector{DeviceStartUpConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        starttime = PSY.get_start_time_limits(d)
        name = PSY.get_name(d)
        start_types = PSY.get_start_types(d)
        range_data = DeviceStartUpConstraintInfo(name, starttime, start_types)
        constraint_data[ix] = range_data
    end

    turbine_temperature(
        psi_container,
        constraint_data,
        make_constraint_name(STARTUP_TIMELIMIT, PSY.ThermalMultiStart),
        make_variable_name(StopVariable, PSY.ThermalMultiStart),
        (
            make_variable_name(HotStartVariable, PSY.ThermalMultiStart),
            make_variable_name(WarmStartVariable, PSY.ThermalMultiStart),
        ),
    )
    return
end

"""
This function creates constraints to select a single type of startup based on off-time
"""
function startup_type_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.ThermalMultiStart},
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {S <: PM.AbstractPowerModel}
    time_steps = model_time_steps(psi_container)
    constraint_data = Vector{DeviceStartTypesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        start_types = PSY.get_start_types(d)
        range_data = DeviceStartTypesConstraintInfo(name, start_types)
        constraint_data[ix] = range_data
    end

    device_start_type_constraint(
        psi_container,
        constraint_data,
        make_constraint_name(START_TYPE, PSY.ThermalMultiStart),
        make_variable_name(START, PSY.ThermalMultiStart),
        (
            make_variable_name(HotStartVariable, PSY.ThermalMultiStart),
            make_variable_name(WarmStartVariable, PSY.ThermalMultiStart),
            make_variable_name(ColdStartVariable, PSY.ThermalMultiStart),
        ),
    )
    return
end

"""
This function gets the data for startup initial condition
"""
function _get_data_startup_ic(initial_conditions::Vector{InitialCondition})
    lenght_devices = length(initial_conditions)
    data = Vector{DeviceStartUpConstraintInfo}(undef, lenght_devices)
    idx = 0
    for ic in initial_conditions
        g = ic.device
        if PSY.get_start_types(g) > 1
            idx = +1
            name = PSY.get_name(g)
            data[idx] = DeviceStartUpConstraintInfo(
                name,
                PSY.get_start_time_limits(g),
                PSY.get_start_types(g),
            )
        end
    end
    if idx < lenght_devices
        deleteat!(data, (idx + 1):lenght_devices)
    end

    return data
end

"""
This function creates the initial conditions for multi-start devices
"""
function startup_initial_condition_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.ThermalMultiStart},
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {S <: PM.AbstractPowerModel}
    time_steps = model_time_steps(psi_container)
    resolution = model_resolution(psi_container)
    key_off = ICKey(TimeDurationOFF, PSY.ThermalMultiStart)
    initial_conditions_offtime = get_initial_conditions(psi_container, key_off)
    constraint_data = _get_data_startup_ic(initial_conditions_offtime)

    device_startup_initial_condition(
        psi_container,
        constraint_data,
        initial_conditions_offtime,
        make_constraint_name(STARTUP_INITIAL_CONDITION, PSY.ThermalMultiStart),
        (
            make_variable_name(HotStartVariable, PSY.ThermalMultiStart),
            make_variable_name(WarmStartVariable, PSY.ThermalMultiStart),
        ),
        make_variable_name(OnVariable, PSY.ThermalMultiStart),
    )
    return
end

"""
This function creates constraints that keep must run devices online
"""
function must_run_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.ThermalMultiStart},
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {S <: PM.AbstractPowerModel}
    time_steps = model_time_steps(psi_container)
    forecast_label = "must_run"
    constraint_infos = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = ones(time_steps[end])
        timeseries_data =
            DeviceTimeSeriesConstraintInfo(d, x -> PSY.get_must_run(x), ts_vector)
        constraint_infos[ix] = timeseries_data
    end
    ts_inputs = TimeSeriesConstraintSpecInternal(
        constraint_infos,
        make_constraint_name(MUST_RUN, PSY.ThermalMultiStart),
        make_variable_name(OnVariable, PSY.ThermalMultiStart),
        nothing,
        nothing,
    )

    device_timeseries_lb!(psi_container, ts_inputs)
    return
end

########################### time duration constraints ######################################
"""
If the fraction of hours that a generator has a duration constraint is less than
the fraction of hours that a single time_step represents then it is not binding.
"""
function _get_data_for_tdc(
    initial_conditions_on::Vector{InitialCondition},
    initial_conditions_off::Vector{InitialCondition},
    resolution::Dates.TimePeriod,
)
    steps_per_hour = 60 / Dates.value(Dates.Minute(resolution))
    fraction_of_hour = 1 / steps_per_hour
    lenght_devices_on = length(initial_conditions_on)
    lenght_devices_off = length(initial_conditions_off)
    @assert lenght_devices_off == lenght_devices_on
    time_params = Vector{UpDown}(undef, lenght_devices_on)
    ini_conds = Matrix{InitialCondition}(undef, lenght_devices_on, 2)
    idx = 0
    for (ix, ic) in enumerate(initial_conditions_on)
        g = ic.device
        @assert g == initial_conditions_off[ix].device
        non_binding_up = false
        non_binding_down = false
        time_limits = PSY.get_time_limits(g)
        name = PSY.get_name(g)
        if !isnothing(time_limits)
            if (time_limits.up <= fraction_of_hour) & (time_limits.down <= fraction_of_hour)
                @debug "Generator $(name) has a nonbinding time limits. Constraints Skipped"
                continue
            else
                idx += 1
            end
            ini_conds[idx, 1] = ic
            ini_conds[idx, 2] = initial_conditions_off[ix]
            up_val = round(time_limits.up * steps_per_hour, RoundUp)
            down_val = round(time_limits.down * steps_per_hour, RoundUp)
            time_params[idx] = time_params[idx] = (up = up_val, down = down_val)
        end
    end
    if idx < lenght_devices_on
        ini_conds = ini_conds[1:idx, :]
        deleteat!(time_params, (idx + 1):lenght_devices_on)
    end
    return ini_conds, time_params
end

function time_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, D <: AbstractThermalFormulation, S <: PM.AbstractPowerModel}
    parameters = model_has_parameters(psi_container)
    resolution = model_resolution(psi_container)
    initial_conditions_on = get_initial_conditions(psi_container, ICKey(TimeDurationON, T))
    initial_conditions_off =
        get_initial_conditions(psi_container, ICKey(TimeDurationOFF, T))
    ini_conds, time_params =
        _get_data_for_tdc(initial_conditions_on, initial_conditions_off, resolution)
    if !(isempty(ini_conds))
        if parameters
            device_duration_parameters!(
                psi_container,
                time_params,
                ini_conds,
                make_constraint_name(DURATION, T),
                (
                    make_variable_name(OnVariable, T),
                    make_variable_name(StartVariable, T),
                    make_variable_name(StopVariable, T),
                ),
            )
        else
            device_duration_retrospective!(
                psi_container,
                time_params,
                ini_conds,
                make_constraint_name(DURATION, T),
                (
                    make_variable_name(OnVariable, T),
                    make_variable_name(StartVariable, T),
                    make_variable_name(StopVariable, T),
                ),
            )
        end
    else
        @warn "Data doesn't contain generators with time-up/down limits, consider adjusting your formulation"
    end
    return
end

function time_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, ThermalMultiStartUnitCommitment},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    parameters = model_has_parameters(psi_container)
    resolution = model_resolution(psi_container)
    initial_conditions_on = get_initial_conditions(psi_container, ICKey(TimeDurationON, T))
    initial_conditions_off =
        get_initial_conditions(psi_container, ICKey(TimeDurationOFF, T))
    ini_conds, time_params =
        _get_data_for_tdc(initial_conditions_on, initial_conditions_off, resolution)
    if !(isempty(ini_conds))
        if parameters
            device_duration_parameters!(
                psi_container,
                time_params,
                ini_conds,
                make_constraint_name(DURATION, T),
                (
                    make_variable_name(OnVariable, T),
                    make_variable_name(StartVariable, T),
                    make_variable_name(StopVariable, T),
                ),
            )
        else
            device_duration_compact_retrospective!(
                psi_container,
                time_params,
                ini_conds,
                make_constraint_name(DURATION, T),
                (
                    make_variable_name(OnVariable, T),
                    make_variable_name(StartVariable, T),
                    make_variable_name(StopVariable, T),
                ),
            )
        end
    else
        @warn "Data doesn't contain generators with time-up/down limits, consider adjusting your formulation"
    end
    return
end

########################### Cost Function Calls#############################################
# These functions are custom implementations of the cost data. In the file cost_functions.jl there are default implementations. Define these only if needed.

function AddCostSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: PSY.ThermalGen, U <: AbstractThermalFormulation}
    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        has_status_variable = has_on_variable(psi_container, T),
        has_status_parameter = has_on_parameter(psi_container, T),
        variable_cost = PSY.get_variable,
        start_up_cost = PSY.get_start_up,
        shut_down_cost = PSY.get_shut_down,
        fixed_cost = PSY.get_fixed,
        sos_status = VARIABLE,
    )
end

function AddCostSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: PSY.ThermalGen, U <: AbstractThermalDispatchFormulation}
    if has_on_parameter(psi_container, T)
        sos_status = PARAMETER
    else
        sos_status = NO_VARIABLE
    end

    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        has_status_variable = has_on_variable(psi_container, T),
        has_status_parameter = has_on_parameter(psi_container, T),
        variable_cost = PSY.get_variable,
        fixed_cost = PSY.get_fixed,
        sos_status = sos_status,
    )
end

function AddCostSpec(
    ::Type{T},
    ::Type{ThermalMultiStartUnitCommitment},
    psi_container::PSIContainer,
) where {T <: PSY.ThermalGen}
    fixed_cost_func = x -> PSY.get_fixed(x) + PSY.get_no_load(x)
    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        has_status_variable = has_on_variable(psi_container, T),
        has_status_parameter = has_on_parameter(psi_container, T),
        #variable_cost = PSY.get_variable, uses SOS by default
        shut_down_cost = PSY.get_shut_down,
        fixed_cost = fixed_cost_func,
        sos_status = VARIABLE,
    )
end

"""
Cost function for generators formulated as No-Min
"""
function cost_function!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, ThermalDispatchNoMin},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen}
    no_min_spec = AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        has_status_variable = has_on_variable(psi_container, T),
        has_status_parameter = has_on_parameter(psi_container, T),
    )

    for g in devices
        component_name = PSY.get_name(g)
        op_cost = PSY.get_operation_cost(g)
        cost_component = PSY.get_variable(op_cost)
        if isa(cost_component, PSY.VariableCost{Array{Tuple{Float64, Float64}, 1}})
            @debug "PWL cost function detected for device $(component_name) using ThermalDispatchNoMin"
            slopes = PSY.get_slopes(cost_component)
            if any(slopes .< 0) || !pwlparamcheck(cost_component)
                throw(IS.InvalidValue("The PWL cost data provided for generator $(PSY.get_name(g)) is not compatible with a No Min Cost."))
            end
            if slopes[1] != 0.0
                @debug "PWL has no 0.0 intercept for generator $(PSY.get_name(g))"
                # adds a first intercept a x = 0.0 and Y below the intercept of the first tuple to make convex equivalent
                first_pair = PSY.get_cost(cost_component)[1]
                cost_function_data = deepcopy(cost_component.cost)
                intercept_point = (0.0, first_pair - COST_EPSILON)
                cost_function_data = vcat(intercept_point, cost_function_data)
                corrected_slopes = PSY.get_slopes(cost_function_data)
                @assert slope_convexity_check(slopes)
            else
                cost_function_data = cost_component.cost
            end
            time_steps = model_time_steps(psi_container)
            for t in time_steps
                pwl_gencost_linear!(
                    psi_container,
                    no_min_spec,
                    component_name,
                    cost_function_data,
                    t,
                )
            end
        else
            add_to_cost!(psi_container, no_min_spec, op_cost, g)
        end
    end
    return
end

"""
Adds to the models costs represented by PowerSystems Market-Bid costs
"""
function add_to_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    cost_data::PSY.MarketBidCost,
    component::PSY.ThermalMultiStart,
)
    component_name = PSY.get_name(component)
    @debug "Market Bid" component_name
    resolution = model_resolution(psi_container)
    dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
    time_steps = model_time_steps(psi_container)
    initial_time = model_initial_time(psi_container)
    # TODO: Use mthods from PowerSystems to eliminate this step
    variable_cost_forecast = get_time_series(psi_container, component, "variable_cost")
    variable_cost_forecast = map(PSY.VariableCost, variable_cost_forecast)
    for t in time_steps
        variable_cost!(psi_container, spec, component_name, variable_cost_forecast[t], t)
    end

    if !isnothing(spec.start_up_cost)
        # Start-up costs
        start_cost_data = PSY.get_start_up(cost_data)
        @debug "start up cost" component_name
        for (st, var_type) in
            enumerate((HotStartVariable, WarmStartVariable, ColdStartVariable))
            var_name = make_variable_name(var_type, spec.component_type)
            for t in time_steps
                linear_gen_cost!(
                    psi_container,
                    var_name,
                    component_name,
                    start_cost_data[st] * spec.multiplier,
                    t,
                )
            end
        end
    end

    if spec.has_status_variable
        @debug "no_load cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(OnVariable, spec.component_type),
                component_name,
                PSY.get_no_load(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    if !isnothing(spec.shut_down_cost)
        @debug "Shut down cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(StopVariable, spec.component_type),
                component_name,
                spec.shut_down_cost(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    #Service Cost Bid
    ancillary_services = PSY.get_ancillary_services(cost_data)
    for service in ancillary_services
        add_service_bid_cost!(psi_container, spec, component, service)
    end
    return
end

# TODO: Define for now just for Area Balance and reason about others later. This will
# be needed and useful for PowerFlow
function NodalExpressionSpec(
    ::Type{T},
    ::Type{AreaBalancePowerModel},
    use_forecasts::Bool,
) where {T <: PSY.ThermalGen}
    return NodalExpressionSpec(
        "max_active_power",
        ACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_active_power(x) : x -> PSY.get_active_power(x),
        1.0,
        T,
    )
end


text:
export TranslationOperation

export apply_operation
export domaintype
export isidentity, istranslation, ispoint

"""
    TranslationOperation{S<:Real}

Represents translation symmetry operation

# Fields
* `displacement`: displacement vector

# Examples
```jldoctest
julia> using TightBindingLattice

julia> TranslationOperation([1, 2])
TranslationOperation{Int64}([1, 2])

julia> TranslationOperation([0.5, 0.0, 0.5])
TranslationOperation{Float64}([0.5, 0.0, 0.5])
```
"""
struct TranslationOperation{S<:Real} <:AbstractSpaceSymmetryOperation{S}
    displacement::Vector{S}

    @doc """
        TranslationOperation{S}(displacement::AbstractVector{<:Real}) where {S}

    Construct a translation operation of `displacement` on coordinate space of type `S`.
    """
    function TranslationOperation{S}(displacement::AbstractVector{<:Real}) where {S}
        new{S}(displacement)
    end

    @doc """
        TranslationOperation(displacement::AbstractVector{S}) where {S<:Real}

    Construct a translation operation of `displacement` on coordinate space of type `S`.
    """
    function TranslationOperation(displacement::AbstractVector{S}) where {S<:Real}
        new{S}(displacement)
    end
end

function Base.convert(::Type{TranslationOperation{S}}, obj::IdentityOperation{S}) where S
    dim = dimension(obj)
    return TranslationOperation{S}(zeros(S, dim))
end

function Base.convert(::Type{TranslationOperation{S}}, displacement::AbstractVector{S}) where S
    return TranslationOperation{S}(displacement)
end

function Base.promote_rule(::Type{TranslationOperation{S}}, ::Type{IdentityOperation{S}}) where S
    return TranslationOperation{S}
end

## properties
"""
    isidentity(t::TranslationOperation)

Return true if the translation operation is an identity, i.e. `iszero(t.displacement)`
"""
isidentity(arg::TranslationOperation) = iszero(arg.displacement)

"""
    istranslation(arg::TranslationOperation)

Always return `true`, since `arg` is already a translation operation.
"""
istranslation(arg::TranslationOperation) = true

"""
    ispoint(arg::TranslationOperation)

Return `true` if `arg` is a point operation. The only way this can be true is when `arg` is
an identity operation.
"""
ispoint(arg::TranslationOperation) = iszero(arg.displacement)

"""
    dimension(arg::TranslationOperation)

Spatial dimension of the translation operation
"""
dimension(arg::TranslationOperation) = length(arg.displacement)

Base.hash(arg::TranslationOperation{S}) where S = hash(arg.displacement, hash(TranslationOperation{S}))

## operators
function Base.:(==)(lhs::TranslationOperation{S}, rhs::TranslationOperation{S}) where S
    return lhs.displacement == rhs.displacement
end
function Base.:(==)(top::TranslationOperation{S}, iden::IdentityOperation{S}) where S
    return iszero(top.displacement)
end
function Base.:(==)(iden::IdentityOperation{S}, top::TranslationOperation{S}) where S
    return iszero(top.displacement)
end

function Base.:(*)(lhs::TranslationOperation{S}, rhs::TranslationOperation{S}) where S
    return TranslationOperation{S}(lhs.displacement .+ rhs.displacement)
end

function Base.:(^)(lhs::TranslationOperation{S}, rhs::Real) where S
    return TranslationOperation{S}(lhs.displacement .* rhs)
end

function Base.inv(arg::TranslationOperation{S}) where S
    return TranslationOperation{S}(-arg.displacement)
end

## apply
"""
    apply_operation(symop::TranslationOperation{S}, coordinate::AbstractArray{S}) where S

Return the translated coordinate.
"""
function apply_operation(symop::TranslationOperation{S},
                         coord::AbstractArray{S}) where {S<:Real}
    return coord .+ symop.displacement
end

function (symop::TranslationOperation{S})(coord::AbstractVector{S}) where S
    return coord .+ symop.displacement
end

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Error

export GitError

@enum(Code, GIT_OK          = Cint(0),   # no error
            ERROR           = Cint(-01), # generic error
            ENOTFOUND       = Cint(-03), # requested object could not be found
            EEXISTS         = Cint(-04), # object exits preventing op
            EAMBIGUOUS      = Cint(-05), # more than one object matches
            EBUFS           = Cint(-06), # output buffer too small to hold data
            EUSER           = Cint(-07), # user callback generated error
            EBAREREPO       = Cint(-08), # operation not allowed on bare repo
            EUNBORNBRANCH   = Cint(-09), # HEAD refers to branch with 0 commits
            EUNMERGED       = Cint(-10), # merge in progress prevented op
            ENONFASTFORWARD = Cint(-11), # ref not fast-forwardable
            EINVALIDSPEC    = Cint(-12), # name / ref not in valid format
            EMERGECONFLICT  = Cint(-13), # merge conflict prevented op
            ELOCKED         = Cint(-14), # lock file prevented op
            EMODIFIED       = Cint(-15), # ref value does not match expected
            EAUTH           = Cint(-16), # authentication error
            ECERTIFICATE    = Cint(-17), # server certificate is invalid
            EAPPLIED        = Cint(-18), # patch/merge has already been applied
            EPEEL           = Cint(-19), # the requested peel operation is not possible
            EEOF            = Cint(-20), # Unexpted EOF
            PASSTHROUGH     = Cint(-30), # internal only
            ITEROVER        = Cint(-31)) # signals end of iteration

@enum(Class, None,
             NoMemory,
             OS,
             Invalid,
             Reference,
             Zlib,
             Repository,
             Config,
             Regex,
             Odb,
             Index,
             Object,
             Net,
             Tag,
             Tree,
             Indexer,
             SSL,
             Submodule,
             Thread,
             Stash,
             Checkout,
             FetchHead,
             Merge,
             SSH,
             Filter,
             Revert,
             Callback,
             CherryPick,
             Describe,
             Rebase)

struct ErrorStruct
    message::Ptr{UInt8}
    class::Cint
end

struct GitError <: Exception
    class::Class
    code::Code
    msg::AbstractString
end
Base.show(io::IO, err::GitError) = print(io, "GitError(Code:$(err.code), Class:$(err.class), $(err.msg))")

function last_error()
    err = ccall((:giterr_last, :libgit2), Ptr{ErrorStruct}, ())
    if err != C_NULL
        err_obj   = unsafe_load(err)
        err_class = Class(err_obj.class)
        err_msg   = unsafe_string(err_obj.message)
    else
        err_class = Class(0)
        err_msg = "No errors"
    end
    return (err_class, err_msg)
end

function GitError(code::Integer)
    err_code = Code(code)
    err_class, err_msg = last_error()
    return GitError(err_class, err_code, err_msg)
end

end # Error module

macro check(git_func)
    quote
        err = Cint($(esc(git_func::Expr)))
        if err < 0
            throw(Error.GitError(err))
        end
        err
    end
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function solve_005()
    D = Dict{Int64, Int64}()
    for i = 1:20
        D[i] = 0
    end
    for n = 4:20
        V = get_prime_factors(n)
        for p = 1:20
            multiplicity = sum(V.==p)
            if multiplicity > D[p] D[p] = multiplicity end
        end
    end
    prod = 1
    for p = 1:20
        prod *= p^D[p]
    end
    return prod
end


text:
include("electrical_components.jl")

R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R = R)
@named capacitor = Capacitor(C = C)
@named source = ConstantVoltage(V = V)
@named ground = Ground()

rc_eqs = [connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n)
          connect(capacitor.n, ground.g)]

@named rc_model = ODESystem(rc_eqs, t)
rc_model = compose(rc_model, [resistor, capacitor, source, ground])


text:
@testset "regression.jl" begin

    # number of input features
    N = 200
    # homoscedastic observation noise
    σ = 0.3

    # experiment of Murray et al.
    @testset "GP regression" begin
        # added to the diagonal of the covariance matrix due to numerical issues of
        # the Cholesky decomposition
        jitter = 1e-9

        # for different dimensions of input features
        for d in 1:10
            # sample input features
            inputs = rand(d, N)

            # define covariance matrix of latent variable
            # add noise to the diagonal due to numerical issues
            prior_Σ =
                Symmetric(exp.((-0.5) .* pairwise(SqEuclidean(), inputs; dims=2))) +
                jitter * I
            prior = MvNormal(prior_Σ)

            # sample noisy observations
            observations = rand(prior) .+ σ .* randn(N)

            # define log likelihood function
            ℓ(f) = logpdf(MvNormal(f, σ^2 * I), observations)

            # run elliptical slice sampler for 100 000 time steps
            samples = sample(ESSModel(prior, ℓ), ESS(), 100_000; progress=false)

            # compute analytical posterior of GP
            posterior_Σ = prior_Σ * (I - (prior_Σ + σ^2 * I) \ prior_Σ)
            posterior_μ = posterior_Σ * observations / σ^2

            # compare with empirical estimates
            @test mean(samples) ≈ posterior_μ rtol = 0.05
        end
    end

    # extreme case with independent observations
    @testset "Independent components" begin
        # define  distribution of latent variables
        prior = MvNormal(Zeros(N), I)

        # sample noisy observations
        observations = rand(prior) .+ σ .* randn(N)

        # define log likelihood function
        ℓ(f) = logpdf(MvNormal(f, σ^2 * I), observations)

        # run elliptical slice sampling for 100 000 time steps
        samples = sample(ESSModel(prior, ℓ), ESS(), 100_000; progress=false)

        # compute analytical posterior
        posterior_μ = observations / (1 + σ^2)

        # compare with empirical estimates
        @test mean(samples) ≈ posterior_μ rtol = 0.02
    end
end


text:
function deref(pointer::Int, program::Array{Int,1}, mode::Int, base::Int)::Int
    if mode == 0
        return program[pointer+1]
    elseif mode == 1
        return pointer
    elseif mode == 2
        return program[pointer+base+1]
    else
        throw(ErrorException("Invalid mode: " * mode))
    end
end

function write_pointer(pointer::Int, mode::Int, base::Int)::Int
    if mode == 0
        return pointer+1
    elseif mode == 2
        return pointer+base+1
    else
        throw(ErrorException("Invalid mode for writing: " * mode))
    end
end

function run_intcode(program::Array{Int,1}, input::Channel, output::Channel)::Int
    program = vcat(program, zeros(Int, 1000))
    pointer = 1
    ret = nothing
    base = 0
    while pointer <= length(program)
        # println(program)
        instruction = reverse(digits(program[pointer], pad=5))
        # println(instruction)
        opcode = instruction[end] + instruction[end-1]*10
        modes = reverse(instruction[1:end-2])
        # println(join([opcode,": ",modes]))
        if opcode == 1 # add
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            program[write_pointer(pt3, modes[3], base)] = val1 + val2
            pointer = pointer + 4
        elseif opcode == 2 # multiply
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            program[write_pointer(pt3, modes[3], base)] = val1 * val2
            pointer = pointer + 4
        elseif opcode == 3 # store input
            pt1 = program[pointer+1]
            program[write_pointer(pt1, modes[1], base)] = take!(input)
            pointer = pointer + 2
        elseif opcode == 4 # output
            pt1 = program[pointer+1]
            val1 = deref(pt1, program, modes[1], base)
            pointer = pointer + 2
            ret = val1
            put!(output, val1)
        elseif opcode == 5 # jump if true
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            if val1 != 0
                pointer = val2 + 1
            else
                pointer = pointer + 3
            end
        elseif opcode == 6 # jump if false
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            if val1 == 0
                pointer = val2 + 1
            else
                pointer = pointer + 3
            end
        elseif opcode == 7 # less than
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            program[write_pointer(pt3, modes[3], base)] = val1 < val2 ? 1 : 0
            pointer = pointer + 4
        elseif opcode == 8 # equal
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            program[write_pointer(pt3, modes[3], base)] = val1 == val2 ? 1 : 0
            pointer = pointer + 4
        elseif opcode == 9 # modify base
            pt1 = program[pointer+1]
            val1 = deref(pt1, program, modes[1], base)
            base = base + val1
            pointer = pointer + 2
        elseif opcode == 99
            return ret
        else
            throw(ErrorException("Invalid opcode: " * string(opcode)))
        end
    end
end

function print_output(ch::Channel)
    while true
        println(take!(ch))
    end
end

program = open("boost.csv") do f
    parse.(Int, split(readlines(f)[1], ","))
end

# program = [109,1,204,-1,1001,100,1,100,1008,100,16,101,1006,101,0,99]
# program = [1102,34915192,34915192,7,4,7,99,0]
# program = [104,1125899906842624,99]

in = Channel(1)
out = Channel(2)

@async print_output(out)
put!(in, 1)
run_intcode(deepcopy(program), in, out)
put!(in, 2)
run_intcode(deepcopy(program), in, out)


text:
export BasicDirectory

"""
    BasicDirectory(map::BlockMap)
Creates a BasicDirectory, which implements the methods of Directory with
basic implmentations
"""
mutable struct BasicDirectory{GID <: Integer, PID <:Integer, LID <: Integer} <: Directory{GID, PID, LID}
    map::BlockMap{GID, PID, LID}

    directoryMap::Union{BlockMap, Nothing}

    procList::Vector{PID}
    procListLists::Vector{Vector{PID}}

    entryOnMultipleProcs::Bool

    localIndexList::Vector{LID}
    allMinGIDs::Vector{GID}

    function BasicDirectory{GID, PID, LID}(map::BlockMap{GID, PID, LID}) where {GID, PID, LID}

        if !(distributedGlobal(map))
            new(map, nothing, [], [], numProc(getComm(map))!=1, [], [])
        elseif linearMap(map)
            commObj = getComm(map)

            allMinGIDs = gatherAll(commObj, minMyGID(map))
            allMinGIDs = vcat(allMinGIDs, [1+maxAllGID(map)])

            entryOnMultipleProcs = length(Set(allMinGIDs)) != length(allMinGIDs)

            new(map, nothing, [], [], entryOnMultipleProcs, [], allMinGIDs)
        else
            generateContent(
                new(map, nothing, [], [], false, [], []),
                map)
        end
    end
end

"""
internal method to assist constructor
"""
function generateContent(dir::BasicDirectory{GID, PID, LID}, map::BlockMap{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
    minAllGID = minAllGID(map)
    maxAllGID = maxAllGID(map)

    dirNumGlobalElements = maxAllGID - minAllGID + 1

    directoryMap = BlockMap(dirNumGlobalElements, minAllGID, commObj)
    dir.directoryMap = directoryMap

    dirNumMyElements = numMyElements(dir.directoryMap)

    if dirNumMyElements > 0
        dir.procList = Vector{PID}(undef, dirNumMyElements)
        dir.localIndexList = Vector{LID}(undef, dirNumMyElements)

        fill!(dir.procList, -1)
        fill!(dir.localIndexList, -1)
    else
        dir.procList = []
        dir.localIndexList = []
    end

    map_numMyElements = numMyElements(map)
    map_myGlobalElements = myGlobalElements(map)

    sendProcs = remoteIDList(directoryMap, map_numMyElements, myGlobalElements(map))

    distributor = createDistributor(commObj)
    numRecvs = createFromSends(distributor, map_numMyElements, sendProcs)

    exportElements = Array{Tuple{GID, PID, LID}}(numMyElements)

    myPIDVal = myPID(commObj)
    for i = 1:numMyElements
        exportElements[i] = (map_myGlobalElements[i], myPIDVal, i)
    end

    importElements = resolve(distributor, exportElements)

    for i = 1:numRecvs
        currLID = lid(directoryMap, importElements[i][1])
        @assert currLID > 0 //internal error

        proc = importElements[i][2]
        if dir.procList[currLID] >= 0
            if dir.procList[currLID] != proc
                if dir.procListLists == []
                    numProcLists = numMyElements(directoryMap)
                    procListLists = Vector{Vector{PID}}(undef, numProcLists)
                    fill!(procListLists, PID[])
                end

                l = procListLists[currLID]

                index = searchsortedfirst(l, procList[currLID])
                insert(l, index, procList[currLID])

                index = searchsortedfirst(l, proc)
                insert(l, index, proc)

                dir.procList[currLID] = dir.procListLists[curr_LID][1]
            end
        else
            dir.procList[currLID] = proc
        end

        dir.localIndexList[currLID] = importElements[i][3]
    end

    globalVal = maxAll(commObj, numProcLists)
    dir.entryOnMultipleProcs = globalval > 0 ? true : false;

    dir
end

function getDirectoryEntries(directory::BasicDirectory{GID, PID, LID}, map::BlockMap{GID, PID, LID}, globalEntries::AbstractVector{GID},
        high_rank_sharing_procs::Bool)::Tuple{Vector{PID}, Vector{LID}} where GID <: Integer where PID <: Integer where LID <: Integer
    numEntries = length(globalEntries)
    procs = Vector{PID}(undef, numEntries)
    localEntries = Vector{LID}(undef, numEntries)

    if !distributedGlobal(map)
        myPIDVal = myPid(getComm(map))

        for i = 1:numEntries
            lidVal = lid(map, globalEntries[i])

            if lidVal == 0
                procs[i] = 0
                warn("GID $(globalEntries[i]) is not part of this map")
            else
                procs[i] = myPIDVal
            end
            localEntries[i] = lidVal
        end
    elseif linearMap(map)
        minAllGIDVal = minAllGID(map)
        maxAllGIDVal = maxAllGID(map)

        numProcVal = numProc(getComm(map))

        n_over_p = numGlobalElements(map)/numProcVal

        allMinGIDs_list = copy(directory.allMinGIDs)
        order = sortperm(allMinGIDs_list)
        permute!(allMinGIDs_list, order)

        for i = 1:numEntries
            lid  = 0
            proc = 0

            gid = globalEntries[i]
            if gid < minAllGIDVal || gid > maxAllGIDVal
                throw(InvalidArgumentError("GID=$gid out of valid range [$minAllGIDVal, $maxAllGIDVal]"))
            end
            #guess uniform distribution and start a little above it
            proc1 = min(GID(fld(gid, max(n_over_p, 1)) + 2), numProcVal)
            proc1 = 1
            found = false

            while proc1 >= 1 && proc1 <= numProcVal
                if allMinGIDs_list[proc1] <= gid
                    if (gid < allMinGIDs_list[proc1+1])
                        found = true
                        break
                    else
                        proc1 += 1
                    end
                else
                    proc1 -= 1
                end
            end
            if found
                proc = order[proc1]
                lid = gid - allMinGIDs_list[proc1] + 1
            end

            procs[i] = proc
            localEntries[i] = lid
        end
    else #general case
        distributor = createDistributor(getComm(map))

        dirProcs = remoteIDList(map, numEntries, globalEntries)

        numMissing = 0
        for i = 1:numEntries
            if dirProcs[i] == 0
                procs[i] = 0
                localEntries[i] = 0
                numMissing += 1
            end
        end

        (sendGIDs, sendPIDs) = createFromRecvs(distrbutor, globalEntries, dirProcs)
        numSends = length(sendGIDs)

        if numSends > 0
            exports = Array{Tuple{GID, PID, LID}}(numSends)
            for i = 1:numSends
                currGID = sendGIDs[i]
                exports[i][1] = currGID

                currLID = lid(map, currGID)
                @assert currLID > 0 #internal error
                if !high_rank_sharing_procs
                    exports[i][2] = procList[currLID]
                else
                    numProcLists = numMyElements(directory.directoryMap)
                    if numProcLists > 0
                        num = length(directory.procListLists[currLID])
                        if num > 1
                            exports[i][2] = directory.procListLists[currLID][num]
                        else
                            exports[i][2] = directory.procList[currLID]
                        end
                    else
                        exports[i][2] = directory.procList[currLID]
                    end
                end
                exports[i][3] = directory.localIndexList[currLID]
            end
        end

        numRecv = numEntries - numMissing
        imports = resolve(distributor, exports)

        offsets = sortperm(globalEntries)
        sortedGE = globalEntries[offsets]

        for i = 1:numRecv
            currLID = imports[i][1]
            j = searchsortedfirst(sortedGE, currLID)
            if j > 0
                procs[offsets[j]] = imports[i][2]
                localEntries[offsets[j]] = imports[i][3]
            end
        end
    end
    (procs, localEntries)
end

function gidsAllUniquelyOwned(directory::BasicDirectory)
    !directory.entryOnMultipleProcs
end


text:
module NeuralNet

export Network, feedforward, backpropagate, cost

type Network
	# array of layer sizes, layers[1] is the size of the input layer
	layers::Array{Int64,1}

	# weights[i] is the weights connecting layer i to i+1
	# so has dimension layers[i+1] x layers[i]
	weights::Array{Array{Float64,2},1}

	# bias[i] is the output bias for neurons in layer i
	bias::Array{Array{Float64,1},1}

	function Network(layers::Array{Int64,1})
		bias = [randn(s) for s in layers[2:end] ]
		weights = [ randn(b,a) for (a,b) in zip(layers[1:end-1], layers[2:end]) ]
		new(layers, weights, bias)
	end
end

function sigmoid(z::Real)
    one(z)/(one(z) + exp(-z))
end
@vectorize_1arg Real sigmoid

function dsigmoid(z::Real)
    exp(-z)/((one(z)+exp(-z))*(one(z)+exp(-z)))
end
@vectorize_1arg Real dsigmoid

function feedforward(n::Network, input)
  # run the input through the neural network
  # i.e., if a[1] = input then iterate
  #   a[i+1] = sigmoid(w[i] * a[i] + b[i])
  a = input
  for idx in 1:size(n.bias,1)
	a = sigmoid( n.weights[idx] * a .+ n.bias[idx] )
  end
  a
end

function sum_weights_sq(n::Network)
    res = 0
    for wt in n.weights
        res += sum(wt.^2)
    end
    res
end

function cost(n::Network, input, output, lambda = 0)
    # standard neural network cost function with optional weight regularization lambda
    noutput = NeuralNet.feedforward(n,input)
    err = (noutput - output).^2
    1./size(input,2) * 0.5 * sum(err) + lambda/2.0 * sum_weights_sq(n)
end

function backpropagate(n::NeuralNet.Network, input, output)
  # Run back-propagation to calculate derivatives of square-error cost without regularization term
  #  n is the network
  #  input is num_features * num_examples
  #  output is num_outpust * num_examples

  num_examples = size(input,2)
  @assert size(input,2) == size(output,2)

  ## Run forward prop and store results
  # a[l] are the activation values
  a = cell(size(n.bias,1)+1)
  # z[l] are the sigmoid inputs
  z = cell(size(n.bias,1))

  a[1] = input
  for idx=1:length(z)
    z[idx] = n.weights[idx] * a[idx] .+ n.bias[idx]
    a[idx+1] = sigmoid( z[idx] )
  end

  ## Backpropogate
  dWeights = cell(size(n.weights,1))
  dBias = cell(size(n.bias,1))

  dBias[end] = (a[end] .- output) .* dsigmoid(z[end])
  dWeights[end] = dBias[end] * a[end-1]'
  for idx=size(dBias,1)-1:-1:1
    dBias[idx] = (n.weights[idx+1]' * dBias[idx+1]) .* dsigmoid(z[idx])
    dWeights[idx] = dBias[idx] * a[idx]'
  end

  # normalize and sum out examples dimension
  for idx in 1:length(dWeights)
    dBias[idx] = 1/num_examples * sum(dBias[idx],2)
    dWeights[idx] = 1/num_examples * dWeights[idx]
  end
  (dBias,dWeights)
end

function dummy_callback(epoch, n::Network)
# do nothing
end

function update_weights!(weight_mat, f1, f2, dweight_mat)
    for idx in 1:length(weight_mat)
        weight_mat[idx] = f1 * weight_mat[idx] - f2 * dweight_mat[idx]
    end
end

function update_bias!(bias, eta, dbias)
    for idx in 1:length(bias)
        bias[idx] = bias[idx] - eta * dbias[idx]
    end
end

function train(n::Network, data, outputs, epochs = 10, batch_size = 1, eta = 0.95, lambda = 0., callback = dummy_callback)
   # Train the network via stochastic gradient descent on training set defined
   # by data and labels
   #
   #  data is num_features * num_samples
   #  outputs is num_outputs  * num_samples
   #  epochs is the number of iteratations to run
   #  batch_size is size of mini-batches in each epoch
   #  eta is learning rate
   #  lambda is regularization coefficient
   #  callback is called at the end of each epoch iteration with the current epoch number
   #        and the current network n

   num_samples = size(data,2)
   @assert size(data,2) == size(outputs,2)

   # For stochastic mini-batch run, we need to pick random subsets of samples
   # However, we do not "own" the data, so we should not shuffle in place.
   # We also do not want to make a copy since data may be large, so instead
   # I make an indexing vector and shuffle that
   idx = [ x for x in 1:num_samples ]
   batch_offsets = vcat([x for x in 1:batch_size:num_samples],[num_samples+1])

   for epoch in 1:epochs
       shuffle!(idx)
       for bid in 1:length(batch_offsets)-1
           batch_ids = idx[batch_offsets[bid]:(batch_offsets[bid+1]-1)]
	       (dBias, dWeights) = backpropagate(n, data[:,batch_ids], outputs[:,batch_ids])
           for layer in 1:length(n.layers)-1
               #n.weights[layer] = ( 1 - eta .* lambda) .* n.weights[layer] - eta .* dWeights[layer]
               update_weights!(n.weights[layer], 1-eta * lambda, eta, dWeights[layer])
               #n.bias[layer] = n.bias[layer] - eta * squeeze(dBias[layer],2)
               update_bias!(n.bias[layer],-eta, dBias[layer])
           end
       end
       callback(epoch, n)
   end
   n
end

end

max_stars_repo_name: bachase/nnadl-julia

max_issues_repo_name: bachase/nnadl-julia

max_forks_repo_name: bachase/nnadl-julia


text:
using Printf

function set_examples()
    println("--------------------------------------")
    println("Set examples")

    println()

    println("Sets are like Arrays with the exception that all items must be unique")
    set1 = Set(["Apple", "Banana", "Carrot", "Donut"])
    set2 = Set(["Banana", "Cheese", "Bread", "Donut"])
    println(set1)
    println(set2)

    println("We can add items with push")
    push!(set1, "Orange")
    println(set1)

    println("We can check to see if an item is in the set")
    println(in("Orange", set1))

    println("We can add sets togethe by using union")
    set3 = union(set1, Set(["Grapefruit", "Beans"]))
    println(set3)

    println("We can find whats in both sets by using intersect")
    println("set1 = ", set1)
    println("set2 = ", set2)
    println(intersect(set1, set2))

    println("We can find whats only in first one but not second")
    println("set1 = ", set1)
    println("set2 = ", set2)
    println(setdiff(set1, set2))
end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
# This is multiplicative analog of diff
function reldiff{T}(v::Vector{T})
    n = length(v)
    res = Array(T, n - 1)
    for i in 2:n
        res[i - 1] = v[i] / v[i - 1]
    end
    return res
end

# Diff scaled by previous value
function percent_change{T}(v::Vector{T})
    n = length(v)
    res = Array(T, n - 1)
    for i in 2:n
        res[i - 1] = (v[i] - v[i - 1]) / v[i - 1]
    end
    return res
end

autocor{T}(dv::DataVector{T}, lag::Int) = cor(dv[1:(end - lag)], dv[(1 + lag):end])
autocor{T}(dv::DataVector{T}) = autocor(dv, 1)

# Generate levels - see the R documentation for gl
function gl(n::Integer, k::Integer, l::Integer)
    nk = n * k
    d, r = divrem(l, nk)
    r == 0 || throw(ArgumentError("length out must be a multiple of n * k"))
    aa = Array(Int, l)
    for j = 0:(d - 1), i = 1:n
        aa[j * nk + (i - 1) * k + (1:k)] = i
    end
    compact(PooledDataArray(aa))
end

gl(n::Integer, k::Integer) = gl(n, k, n*k)

# A cross-tabulation type. Currently just a one-way table
type xtab{T}
    vals::Array{T}
    counts::Vector{Int}
end

function xtab{T}(x::AbstractArray{T})
    d = Dict{T, Int}()
    for el in x
        d[el] = get(d, el, 0) + 1
    end
    kk = sort(keys(d))
    cc = Array(Int, length(kk))
    for i in 1:length(kk)
        cc[i] = d[kk[i]]
    end
    return xtab(kk, cc)
end

# Another cross-tabulation function, this one leaves the result as a Dict
# Again, this is currently just for one-way tables.
function xtabs{T}(x::AbstractArray{T})
    d = Dict{T, Int}()
    for el in x
        d[el] = get(d, el, 0) + 1
    end
    return d
end


text:
import Base.*, Base.sqrt
export langevinSampler, IdLangevin

struct IdLangevin end
function *(P::IdLangevin, x)
    return x
end
function sqrt(P::IdLangevin)
    return IdLangevin()
end

function langevinSampler(
    neglogprob::Function,
    nsamples::Int64, x0::Array{Float32}, ϵ::Union{Float32, NTuple{3, Float32}};
    flag_MALA::Bool = false,
    P::Union{IdLangevin, AbstractArray{Float32, 2}} = IdLangevin(),
    verbose::Bool = true)

    # Initialize output variables
    fval = zeros(Float32, nsamples)
    samples = Array{Float32, 2}(undef, prod(size(x0)), nsamples)

    # Computing sqrt of preconditioner
    sqrtP = sqrt(P)

    # Noise parameters
    ϵ_vec = compute_ϵvec(ϵ, nsamples)

    # Starting model
    x = x0
    samples[:, 1] = x0

    # Loop over iterations
    n = 1
    accmem = 10
    accrej_array = Array{Bool, 1}(undef, accmem); accrej_array .= false
    while n < nsamples

        # Evaluate objective and gradients
        fval[n], g = neglogprob(x)

        # Print current iteration
        if verbose
            if !flag_MALA
                println("Iter: [", n, "/", nsamples, "] --- fval = ", fval[n])
            else
                println("Iter: [", n, "/", nsamples, "] --- fval = ", fval[n], ", acc_rate = ", sum(accrej_array)/accmem)
            end
        end

        # Proposed model
        Pg = P*g
        x_ = x-ϵ_vec[n]*Pg+sqrt(2*ϵ_vec[n])*(sqrtP*randn(Float32, size(x)))

        # Acceptance-rejection step
        accrej_array[1:end-1] = accrej_array[2:end]
        if flag_MALA
            neglogq = Float32(0.25/ϵ_vec[n])*norm(x_-x+ϵ_vec[n]*Pg)^2
            fval_, g_ = neglogprob(x_); Pg_ = P*g_
            neglogq_ = Float32(0.25/ϵ_vec[n])*norm(x-x_+ϵ_vec[n]*Pg_)^2
            α = min(Float32(0), -(fval_+neglogq-fval[n]-neglogq_))
            if α >= log(rand(Float32))
                x = x_
                accrej_array[end] = true
            else
                n -= 1
                accrej_array[end] = false
            end
        else
            x = x_
        end

        # Collect sample
        n += 1
        samples[:, n] = x

    end

    return fval, samples

end

function compute_ϵvec(ϵ::Union{Float32, NTuple{3, Float32}}, n::Int64)
    if isa(ϵ, NTuple{3, Float32})
        b = (n-1)/((ϵ[1]/ϵ[2])^(1/ϵ[3])-1)
        a = b^gamma*ϵ0
        return a./(b.+collect(1:n).-1).^ϵ[3]
    else
        return ϵ*ones(Float32, n)
    end
end


text:
const COPY_SYMBOLS = (-1:9...)
const REVERSE_SYMBOLS = (-2:9...)
const WALK_SYMBOLS = (-4:9...)
const SYMBOLS = (0:9...)
const NO_SYMBOL = -1
const WRITE = "write"
const NOT_WRITE = "not-write"

const TAPE_ACTIONS = ("mr","ml")
const GRID_ACTIONS = ("mr","ml")
const WRITE_ACTIONS = (WRITE, NOT_WRITE)
const STOP_ACTION = ("<s>", NOT_WRITE)

function get_symbols(task)
    if in(task,("copy","reverse","walk"))
        return eval(parse(uppercase(string(task,"_symbols"))))
    end
    return SYMBOLS
end

function get_actions(task)
    actions = nothing
    if in(task,("copy","reverse"))
        actions = TAPE_ACTIONS
    else
        actions = GRID_ACTIONS
    end
    retval = []
    for ma in actions
        for wa in WRITE_ACTIONS
            push!(retval, (ma,wa))
        end
    end
    return [retval..., STOP_ACTION]
end

function initvocab(symbols)
    symbols = collect(symbols)
    s2i, i2s = Dict(), Dict()
    c = 1
    for sym in symbols
        s2i[sym] = c
        i2s[c] = sym
        c += 1
    end
    return s2i, i2s
end

max_stars_repo_name: ereday/algorithm-learning

max_issues_repo_name: ereday/algorithm-learning

max_forks_repo_name: ereday/algorithm-learning


text:
module TestScratchSpace

using FLoops: Cleared, ScratchSpace, allocate
using Serialization
using Test

function roundtrip(x)
    @nospecialize
    buf = IOBuffer()
    serialize(buf, x)
    seek(buf, 0)
    return deserialize(buf)
end

function check_invariance(x)
    @nospecialize
    @assert !(x.value isa Cleared)
    @testset "allocate" begin
        @test allocate(x) === x
    end
    @testset "serialize: y = roundtrip(x)" begin
        y = roundtrip(x)
        @test y isa typeof(x)
        @test y.f === x.f
        @test y.value === Cleared()
        @testset "re-allocate: z = allocate(y)" begin
            z = allocate(y)
            @test z isa typeof(x)
            @test z.f === x.f
            @test z.value isa typeof(x.value)
        end
    end
    @testset "s = Some(x); t = roundtrip(s)" begin
        s = Some(x)
        t = roundtrip(s)
        @test typeof(t) === typeof(s)
        @test something(t).f === x.f
        @test something(t).value === Cleared()
    end
    @testset "y = roundtrip((x,))[1]" begin
        y = roundtrip((x,))[1]
        @test y isa typeof(x)
        @test y.f === x.f
        @test y.value === Cleared()
    end
end

makevector() = [123]
makeint() = 123

function test_invariance_mutable()
    @testset "x = ScratchSpace(makevector, makevector())" begin
        x = ScratchSpace(makevector, makevector())
        check_invariance(x)
    end
end

function test_invariance_immutable()
    @testset "x = ScratchSpace(makeint, makeint())" begin
        x = ScratchSpace(makeint, makeint())
        check_invariance(x)
    end
end

end  # module


text:
include("jump_maximin.jl")
include("jump_robust_maximin.jl")


text:
module Benchmark
using LRUCache

# This benchmark is a simple recursive Fibonnaci calculation. While not
# comparable to most real world problems, the simplicity of the underlying
# calculation means that most of the time is spent in the cacheing operations,
# rather than on the function itself. Which is optimal for comparing
# improvements to the cache speed.

const FIBCACHE = LRU{Int, Int}(10)

function fib(a::Int)
    @get! FIBCACHE a begin
        if a < 2
            a
        else
            fib(a - 1) + fib(a - 2)
        end
    end
end

function fib_benchmark(cachesize)
    resize!(FIBCACHE, cachesize)
    empty!(FIBCACHE)
    println("Cache Size = $cachesize")
    @time fib(50)
end

println("========== BENCHMARKS ==========")
println()

println("Toss this timing, things are still compiling")
fib_benchmark(5)
println()
println("Fibonnaci Benchmarks")
println("--------------------")
for n in (10, 50, 100)
    fib_benchmark(n)
end
println()

# Now we benchmark individual operations
function setup_cache(cachesize, items)
    resize!(FIBCACHE, cachesize)
    empty!(FIBCACHE)
    for i in items
        FIBCACHE[i] = i
    end
end
function benchmark_get(n)
    @time FIBCACHE[n]
end
function benchmark_set(k, v)
    @time FIBCACHE[k] = v
end

println("Access Benchmarks")
println("-----------------")
# Cache is [10, 9, 8, ...., 1] in order of use, and is full
setup_cache(10, 1:10)
println("Access for elements at head of cache")
for i in 1:5
    benchmark_get(10)
end
println()

# Cache is [10, 9, 8, ...., 1] in order of use, and is full
println("Access for elements not at head of cache")
for i in 1:5
    benchmark_get(i)
end
println()

println("Insertion Benchmarks")
println("--------------------")
# Cache is empty, with size 5
resize!(FIBCACHE, 5)
empty!(FIBCACHE)
println("Insertion when cache is not full, element not in cache")
for i in 1:5
    benchmark_set(i, i)
end
println()

# Cache is [5, 4, 3, 2, 1]
println("Insertion when element already in cache, not at head")
for i in 1:5
    benchmark_set(i, i)
end
println()

# Cache is [5, 4, 3, 2, 1]
println("Insertion when element already in cache, at head")
for i in 1:5
    benchmark_set(5, 5)
end
println()

# Cache is [5, 4, 3, 2, 1]
println("Insertion when cache is full, element not in cache")
for i in 6:10
    benchmark_set(i, i)
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9d286b46-73ae-49f7-b834-791b9ef95420
begin
	using PlutoUI
end

# ╔═╡ 4abf8eda-2c36-11ec-3cdb-4104751a59ad
md"""

#### Solution: Flash Practice Problem

"""

# ╔═╡ c017ef4d-a828-4f46-922a-453d72636c39
md"""$(PlutoUI.LocalResource("./figs/Fig-VLE-Ideal-Pxy-P2-F19.png"))"""

# ╔═╡ d5b781d8-6e6a-4bc9-9b45-1272e1b86e61
md"""
a) What temperature $T$ (K) is the Flash drum operating at?

To find the temperature, we can solve the Antoine equation for $T$.
As a reminder, the Antoine equation is given by (note this is the $\log_{10}$ version):

$$\begin{equation}
\log_{10}\left(P_{i}^{sat}~[\mathrm{bar}]\right) = A - \frac{B}{C+T[K]}
\end{equation}$$

where the Antoine parameters are given in the problem. We can solve for $T$ to get the expression:

$$T = - \left(C + \frac{B}{\log_{10}P^{sat}_{i} - A}\right)$$

We can choose _either_ pure component endpoints, let's choose pure component 1 which gives $P_{1}^{sat}\simeq$ 15.9 bar.
"""

# ╔═╡ e86e41b3-d39b-4986-b1cd-4ed34eaabae5
begin

	# what are the A, B, and C for component 1?
	A₁ = 4.45 		# units: AU
	B₁ = 718.1 		# units: 1/K
	C₁ = -22.01 	# units: K

	A₂ = 3.98 		# units: AU
	B₂ = 677.1 		# units: 1/K
	C₂ = -24.51 	# units: K

	P1sat = 15.9 	# units: bar
	P2sat = 7.5 	# units: bar
	P = 10.5 		# units: bar

	# esimate T -
	T = -(C₁ + B₁/(log10(P1sat) - A₁))

	with_terminal() do
		println("T = $(round(T; digits=2)) K")
	end
end

# ╔═╡ 796fd591-bf7f-4b94-af8f-f5b7ab25e563
md"""
b) Estimate the outlet composition and the mol flow rates for the input, liqud and vapor streams

The equlibrium compositions can be estimated directly from the _Pxy_ diagram:

* The composition of the liquid stream exiting the drum $x^{eq}_{1}\simeq$ 0.34
* The composition of the vapor stream exiting the drum $y^{eq}_{1}\simeq$ 0.53

To estimate the mol flow rate of liquid and vapor that leaves the flash drum, we can use the lever rule and the species mol balance (as discussed in lecture and discussion). The Flash Lever Rule (FLR) is given as:

$$\frac{\dot{L}}{\dot{V}} = \frac{y^{eq}_{i} - z_{i}}{z_{i} - x_{i}^{eq}}$$

_Substitution_: We can re-arrange the FLR, to solve for $\dot{L}$ in terms of $\dot{V}$, and then subsutute that into the total mol balance to solve for $\dot{V}$:

$$\dot{V} = \dot{F}\left(\frac{b}{a} + 1\right)^{-1}$$
Once we have $\dot{V}$ we can then solve for $\dot{L}$.

_Matrix-vector_: Alternatively, we could set this problem up in matrix vector form:

$$\begin{pmatrix}
1 & 1 \\
1 & -\frac{b}{a}
\end{pmatrix}
\begin{pmatrix}
\dot{L} \\
\dot{V}
\end{pmatrix} =
\begin{pmatrix}
\dot{F} \\
0
\end{pmatrix}$$

and solve for $\dot{L}$ and $\dot{V}$ directly by inverting the 2$\times$2 matrix.

"""

# ╔═╡ 453aba4c-fc32-46a7-ad1c-0c360a5dc886
begin

	# from the diagram -
	y1_eq = 0.53
	x1_eq = 0.34
	z1 = 0.42
	F_dot = 10 # units: mol/t

	# compte b and a -
	b = (y1_eq - z1)
	a = (z1 - x1_eq)

	# setup Ax = b
	A = [1 1 ; 1 -(b/a)]
	b = [F_dot; 0];

	# solve for L_dot and V_dot -
	x_graph = inv(A)*b
	L_dot = x_graph[1]
	V_dot = x_graph[2]

	with_terminal() do
		println("L_dot = $(round(L_dot; digits=1)) mol/t and V_dot = $(round(V_dot; digits=1)) mol/t")
	end
end

# ╔═╡ c4f4a653-cd23-4c9a-9b3d-10da27ca9d79
md"""
c) Numerical check on graphically values. If the pressure summation expression:

$$\begin{equation}
        \sum_{i=1}^{\mathcal{M}}z_{i}\left[\frac{P_{i}^{sat}}{P}\left(\frac{\dot{V}}{\dot{F}}\right)+\frac{\dot{L}}{\dot{F}}\right]^{-1} = 1
\end{equation}$$

is _significantly_ different than 1 (greater than $\pm$ 10% difference), please re-estimate your graphical values (show your work).
"""

# ╔═╡ 2c1378b9-08e3-42ef-9239-fd86fd85cd38
begin

	# we only have two terms, since we are binary
	term_1 = z1*((P1sat/P)*(V_dot/F_dot)+(L_dot/F_dot))^-1
	term_2 = (1-z1)*((P2sat/P)*(V_dot/F_dot)+(L_dot/F_dot))^-1
	test_value = (term_1 + term_2)
	frac_error = (1.0 - test_value)*100

	with_terminal() do
		println("Percentage diff: $(round(frac_error; digits=2))%")
	end
end

# ╔═╡ 91f257b4-b480-4179-b68a-760824600949
md"""

d) _Analytically_ check your graphical composition estimates by computing $x^{eq}_{1}$ and $y^{eq}_{1}$ using the your VLE model and mol balances.

Strategy:
* Solve for $P_{i}^{sat}$ using the estimated $T$ from part a)
* Solve for $x^{eq}_{1}$ using the pressure expression (then $x^{eq}_{2}$ = $1 - x^{eq}_{1}$).
* Solve for $y^{eq}_{1}$ using the VLE model: $y^{eq}_{1}P$ = $x^{eq}_{1}P_{1}^{sat}$. (and then $y^{eq}_{2}$ = $1 - y^{eq}_{1}$).
* Solve 2$\times$2 system for $\dot{L}$ and $\dot{V}$ (either 2 species balances, or 1 species balance and the total mol balance)

"""

# ╔═╡ f167be1a-2801-4cc6-b036-6d6688dddb44
begin

	# 1) compute the saturation pressures -
	P1sat_aeq = 10^(A₁ - (B₁/(C₁ + T)))
	P2sat_aeq = 10^(A₂ - (B₂/(C₂ + T)))

	with_terminal() do
		println("P1sat_aeq = $(round(P1sat_aeq; digits=2)) bar and P2sat_aeq = $(round(P2sat_aeq; digits=2)) bar")
	end
end

# ╔═╡ b8002941-f892-46b8-9c27-1cb08cc3c71b
begin

	# 2) estimate x1_eq from the pressure expression -
	x1_eq_a = (P - P2sat_aeq)/(P1sat_aeq - P2sat_aeq)
	x2_eq_a = 1 - x1_eq_a

	with_terminal() do
		println("x1_eq = $(round(x1_eq_a; digits=2)) and x2_eq = $(round(x2_eq_a; digits=2))")
	end
end

# ╔═╡ 52fba159-cafc-4293-ab32-6b965ca957b8
begin

	# 3) estimate y1_eq and y2_eq -
	y1_eq_a = (x1_eq_a)*(P1sat_aeq/P)
	y2_eq_a = 1 - y1_eq_a

	with_terminal() do
		println("y1_eq = $(round(y1_eq_a; digits=2)) and y2_eq = $(round(y2_eq_a; digits=2))")
	end
end

# ╔═╡ ed9fe0d8-0d9b-4ad2-b327-44d9798cc666
begin

	# 4) solve for the flow rates. We could do this as a traditional substitution, but lets be quick and setup as matrix vector
	A_flow = [x1_eq_a y1_eq_a ; x2_eq_a y2_eq_a]
	b_flow = [F_dot*z1 ; F_dot*(1 - z1)]
	flow_vec = inv(A_flow)*b_flow

	# out -
	with_terminal() do
		println("L_dot = $(round(flow_vec[1]; digits=1)) mol/t and V_dot = $(round(flow_vec[2]; digits=1)) mol/t")
	end
end

# ╔═╡ dbfacc23-773b-49bb-a531-f460aacf5f98
html"""
<style>
main {
    max-width: 1200px;
    width: 90%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.16"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "f6532909bf3d40b308a0f360b6a0e626c0e263a8"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.1"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.5"

[[PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─4abf8eda-2c36-11ec-3cdb-4104751a59ad
# ╟─c017ef4d-a828-4f46-922a-453d72636c39
# ╠═d5b781d8-6e6a-4bc9-9b45-1272e1b86e61
# ╠═e86e41b3-d39b-4986-b1cd-4ed34eaabae5
# ╟─796fd591-bf7f-4b94-af8f-f5b7ab25e563
# ╠═453aba4c-fc32-46a7-ad1c-0c360a5dc886
# ╟─c4f4a653-cd23-4c9a-9b3d-10da27ca9d79
# ╠═2c1378b9-08e3-42ef-9239-fd86fd85cd38
# ╟─91f257b4-b480-4179-b68a-760824600949
# ╠═f167be1a-2801-4cc6-b036-6d6688dddb44
# ╠═b8002941-f892-46b8-9c27-1cb08cc3c71b
# ╠═52fba159-cafc-4293-ab32-6b965ca957b8
# ╠═ed9fe0d8-0d9b-4ad2-b327-44d9798cc666
# ╟─dbfacc23-773b-49bb-a531-f460aacf5f98
# ╟─9d286b46-73ae-49f7-b834-791b9ef95420
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using LearnJulia
using Test

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

max_stars_repo_name: BradLyman/learn-julia

max_issues_repo_name: BradLyman/learn-julia

max_forks_repo_name: BradLyman/learn-julia


text:
struct Labeled{T}
    obs::T
    label::Union{Missing,Int}
end

Base.broadcastable(x::Labeled) = Ref(x)

Labeled(x) = Labeled(x, missing)

putatives(p::InfiniteParticle, x::Labeled) =
    x.label === missing ? putatives(p, x.obs) : (fit(p, x.obs, x.label), )