text:
#! format: off
requires_initialization(::AbstractHybridFormulation) = false

get_variable_multiplier(_, ::Type{<:PSY.HybridSystem}, ::AbstractHybridFormulation) = 1.0
get_expression_type_for_reserve(::ActivePowerReserveVariable, ::Type{<:PSY.HybridSystem}, ::Type{<:PSY.Reserve{PSY.ReserveUp}}) = ComponentReserveUpBalanceExpression
get_expression_type_for_reserve(::ActivePowerReserveVariable, ::Type{<:PSY.HybridSystem}, ::Type{<:PSY.Reserve{PSY.ReserveDown}}) = ComponentReserveDownBalanceExpression

########################### ActivePowerOutVariable, HybridSystem #################################
get_variable_binary(::ActivePowerVariable, ::Type{PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_warm_start_value(::ActivePowerVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_active_power(d)
get_variable_lower_bound(::ActivePowerVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = -1.0 * PSY.get_input_active_power_limits(d).max
get_variable_lower_bound(::ActivePowerVariable, d::PSY.HybridSystem, ::AbstractStandardHybridFormulation) = PSY.get_output_active_power_limits(d).min
get_variable_upper_bound(::ActivePowerVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_output_active_power_limits(d).max

############## ComponentActivePowerVariable, HybridSystem ####################
get_variable_binary(::ComponentActivePowerVariable, ::Type{PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_lower_bound(::ComponentActivePowerVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = 0.0

############## ActivePowerInVariable, HybridSystem ####################
get_variable_binary(::ActivePowerInVariable, ::Type{PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_lower_bound(::ActivePowerInVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_input_active_power_limits(d).min
get_variable_upper_bound(::ActivePowerInVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_input_active_power_limits(d).max

############## ActivePowerOutVariable, HybridSystem ####################
get_variable_binary(::ActivePowerOutVariable, ::Type{PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_upper_bound(::ActivePowerOutVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_output_active_power_limits(d).max
get_variable_lower_bound(::ActivePowerOutVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_output_active_power_limits(d).min

############## EnergyVariable, HybridSystem ####################
get_variable_binary(::EnergyVariable, ::Type{PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_upper_bound(::EnergyVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_state_of_charge_limits(PSY.get_storage(d)).max
get_variable_lower_bound(::EnergyVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_state_of_charge_limits(PSY.get_storage(d)).min
get_variable_warm_start_value(::EnergyVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_initial_energy(PSY.get_storage(d))

############## ReactivePowerVariable, HybridSystem ####################
get_variable_binary(::ReactivePowerVariable, ::Type{PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_upper_bound(::ReactivePowerVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_reactive_power_limits(PSY.get_storage(d)).max
get_variable_lower_bound(::ReactivePowerVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_reactive_power_limits(PSY.get_storage(d)).min
get_variable_warm_start_value(::ReactivePowerVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_reactive_power(PSY.get_storage(d))

############## ComponentReactivePowerVariable, ThermalGen ####################
get_variable_binary(::ComponentReactivePowerVariable, ::Type{PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_lower_bound(::ComponentReactivePowerVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = 0.0

############## ComponentActivePowerReserveUpVariable, HybridSystem ####################
get_variable_binary(::ComponentActivePowerReserveUpVariable, ::Type{<:PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_lower_bound(::ComponentActivePowerReserveUpVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = 0.0

############## ComponentActivePowerReserveDownVariable, HybridSystem ####################
get_variable_binary(::ComponentActivePowerReserveDownVariable, ::Type{<:PSY.HybridSystem}, ::AbstractHybridFormulation) = false
get_variable_lower_bound(::ComponentActivePowerReserveDownVariable, d::PSY.HybridSystem, ::AbstractHybridFormulation) = 0.0

############## ReservationVariable, HybridSystem ####################
get_variable_binary(::ReservationVariable, ::Type{<:PSY.HybridSystem}, ::AbstractHybridFormulation) = true

#################### Initial Conditions for models ###############

initial_condition_default(::InitialEnergyLevel, d::PSY.HybridSystem, ::AbstractHybridFormulation) = PSY.get_initial_energy(PSY.get_storage(d))
initial_condition_variable(::InitialEnergyLevel, d::PSY.HybridSystem, ::AbstractHybridFormulation) = EnergyVariable()

########################Objective Function##################################################
objective_function_multiplier(::VariableType, ::AbstractHybridFormulation)=OBJECTIVE_FUNCTION_POSITIVE

proportional_cost(cost::PSY.OperationalCost, ::OnVariable, ::PSY.HybridSystem, ::AbstractHybridFormulation)=PSY.get_fixed(cost)
proportional_cost(cost::PSY.MarketBidCost, ::OnVariable, ::PSY.HybridSystem, ::AbstractHybridFormulation)=PSY.get_no_load(cost)

sos_status(::PSY.HybridSystem, ::AbstractHybridFormulation)=SOSStatusVariable.NO_VARIABLE

uses_compact_power(::PSY.HybridSystem, ::AbstractHybridFormulation)=false
variable_cost(cost::PSY.OperationalCost, ::PSY.HybridSystem, ::AbstractHybridFormulation)=PSY.get_variable(cost)

get_multiplier_value(::ActivePowerTimeSeriesParameter, d::PSY.HybridSystem, ::Type{<:PSY.RenewableGen}, ::AbstractHybridFormulation) = PSY.get_max_active_power(PSY.get_renewable_unit(d))
get_multiplier_value(::ActivePowerTimeSeriesParameter, d::PSY.HybridSystem, ::Type{<:PSY.ElectricLoad}, ::AbstractHybridFormulation) = PSY.get_max_active_power(PSY.get_electric_load(d))
#! format: on
get_initial_conditions_device_model(
    ::OperationModel,
    ::DeviceModel{T, <:AbstractHybridFormulation},
) where {T <: PSY.HybridSystem} = DeviceModel(T, BasicHybridDisaptch)

does_subcomponent_exist(v::PSY.HybridSystem, ::Type{PSY.ThermalGen}) =
    !isnothing(PSY.get_thermal_unit(v))
does_subcomponent_exist(v::PSY.HybridSystem, ::Type{PSY.RenewableGen}) =
    !isnothing(PSY.get_renewable_unit(v))
does_subcomponent_exist(v::PSY.HybridSystem, ::Type{PSY.ElectricLoad}) =
    !isnothing(PSY.get_electric_load(v))
does_subcomponent_exist(v::PSY.HybridSystem, ::Type{PSY.Storage}) =
    !isnothing(PSY.get_storage(v))

get_subcomponent(v::PSY.HybridSystem, ::Type{PSY.ThermalGen}) = PSY.get_thermal_unit(v)
get_subcomponent(v::PSY.HybridSystem, ::Type{PSY.RenewableGen}) = PSY.get_renewable_unit(v)
get_subcomponent(v::PSY.HybridSystem, ::Type{PSY.ElectricLoad}) = PSY.get_electric_load(v)
get_subcomponent(v::PSY.HybridSystem, ::Type{PSY.Storage}) = PSY.get_storage(v)

function make_subsystem_time_series_name(subcomponent::PSY.Component, label::String)
    return make_subsystem_time_series_name(typeof(subcomponent), label)
end

function make_subsystem_time_series_name(subcomponent::Type{<:PSY.Component}, label::String)
    return IS.strip_module_name(subcomponent) * "__" * label
end

function get_default_time_series_names(
    ::Type{<:PSY.HybridSystem},
    ::Type{<:Union{FixedOutput, AbstractHybridFormulation}},
)
    return Dict{Type{<:TimeSeriesParameter}, String}(
        ActivePowerTimeSeriesParameter => "max_active_power",
    )
end

function get_default_attributes(
    ::Type{<:PSY.HybridSystem},
    ::Type{<:AbstractHybridFormulation},
)
    return Dict{String, Any}("reservation" => true)
end

################################ output power constraints ###########################

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{<:ReactivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_reactive_power_limits(device)

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{InputActivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_input_active_power_limits(PSY.get_storage(device))

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{OutputActivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_output_active_power_limits(PSY.get_storage(device))

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{PSY.RenewableGen},
    ::Type{ComponentActivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_active_power_limits(PSY.get_renewable_unit(device))

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{ActivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = (
    min=-1 * PSY.get_input_active_power_limits(device).max,
    max=PSY.get_output_active_power_limits(device).max,
)

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{PSY.ThermalGen},
    ::Type{ComponentActivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_active_power_limits(PSY.get_thermal_unit(device))

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{PSY.ThermalGen},
    ::Type{ComponentReactivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_reactive_power_limits(PSY.get_thermal_unit(device))

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{PSY.RenewableGen},
    ::Type{ComponentReactivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_reactive_power_limits(PSY.get_renewable_unit(device))

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{PSY.Storage},
    ::Type{ComponentReactivePowerVariableLimitsConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_reactive_power_limits(PSY.get_storage(device))

get_min_max_limits(
    device::PSY.HybridSystem,
    ::Type{EnergyCapacityConstraint},
    ::Type{<:AbstractHybridFormulation},
) = PSY.get_state_of_charge_limits(PSY.get_storage(device))

function add_lower_bound_range_constraints_impl!(
    container::OptimizationContainer,
    T::Type{ComponentActivePowerVariableLimitsConstraint},
    array,
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {V <: PSY.HybridSystem, W <: AbstractDeviceFormulation}
    constraint = T()
    component_type = V
    time_steps = get_time_steps(container)
    device_names =
        [PSY.get_name(d) for d in devices if does_subcomponent_exist(d, PSY.ThermalGen)]

    con_lb = add_constraints_container!(
        container,
        constraint,
        component_type,
        device_names,
        time_steps,
        meta="lb",
    )

    for device in devices, t in time_steps
        !does_subcomponent_exist(device, PSY.ThermalGen) && continue
        ci_name = PSY.get_name(device)
        subcomp_key = string(PSY.ThermalGen)
        limits = get_min_max_limits(device, PSY.ThermalGen, T, W) # depends on constraint type and formulation type
        con_lb[ci_name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[ci_name, subcomp_key, t] >= limits.min
        )
    end
end

function add_upper_bound_range_constraints_impl!(
    container::OptimizationContainer,
    T::Type{ComponentActivePowerVariableLimitsConstraint},
    array,
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {V <: PSY.HybridSystem, W <: AbstractDeviceFormulation}
    constraint = T()
    component_type = V
    time_steps = get_time_steps(container)
    device_names =
        [PSY.get_name(d) for d in devices if does_subcomponent_exist(d, PSY.ThermalGen)]

    con_ub = add_constraints_container!(
        container,
        constraint,
        component_type,
        device_names,
        time_steps,
        meta="ub",
    )

    for device in devices, t in time_steps
        !does_subcomponent_exist(device, PSY.ThermalGen) && continue
        ci_name = PSY.get_name(device)
        subcomp_key = string(PSY.ThermalGen)
        limits = get_min_max_limits(device, PSY.ThermalGen, T, W) # depends on constraint type and formulation type
        con_ub[ci_name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[ci_name, subcomp_key, t] <= limits.max
        )
    end
end

function add_parameterized_upper_bound_range_constraints_impl!(
    container::OptimizationContainer,
    T::Type{ComponentActivePowerVariableLimitsConstraint},
    array,
    P::Type{<:ParameterType},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {V <: PSY.HybridSystem, W <: AbstractDeviceFormulation}
    time_steps = get_time_steps(container)
    constraint = T()
    component_type = V
    names =
        [PSY.get_name(d) for d in devices if does_subcomponent_exist(d, PSY.RenewableGen)]

    constraint = add_constraints_container!(
        container,
        constraint,
        component_type,
        names,
        time_steps,
        meta="re ub",
    )

    parameter = get_parameter_array(container, P(), V)
    multiplier = get_parameter_multiplier_array(container, P(), V)
    for device in devices, t in time_steps
        !does_subcomponent_exist(device, PSY.RenewableGen) && continue
        subcomp_key = string(PSY.RenewableGen)
        name = PSY.get_name(device)
        constraint[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[name, subcomp_key, t] <=
            multiplier[name, subcomp_key, t] * parameter[name, subcomp_key, t]
        )
    end
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{<:PowerVariableLimitsConstraint},
    U::Type{<:Union{VariableType, ExpressionType}},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation}
    add_range_constraints!(container, T, U, devices, model, X)
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{ComponentActivePowerVariableLimitsConstraint},
    U::Type{<:VariableType},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation}
    array = get_variable(container, U(), V)
    add_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    add_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    add_parameterized_upper_bound_range_constraints_impl!(
        container,
        T,
        array,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
    )
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{ComponentActivePowerVariableLimitsConstraint},
    U::Type{<:RangeConstraintLBExpressions},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation}
    array = get_expression(container, U(), V)
    add_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{ComponentActivePowerVariableLimitsConstraint},
    U::Type{<:RangeConstraintUBExpressions},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation}
    array = get_expression(container, U(), V)
    add_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    add_parameterized_upper_bound_range_constraints_impl!(
        container,
        T,
        array,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
    )
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{InputActivePowerVariableLimitsConstraint},
    U::Type{<:Union{VariableType, ExpressionType}},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation}
    if get_attribute(model, "reservation")
        add_reserve_range_constraints!(container, T, U, devices, model, X)
    else
        add_range_constraints!(container, T, U, devices, model, X)
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{OutputActivePowerVariableLimitsConstraint},
    U::Type{<:Union{VariableType, ExpressionType}},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation}
    if get_attribute(model, "reservation")
        add_reserve_range_constraints!(container, T, U, devices, model, X)
    else
        add_range_constraints!(container, T, U, devices, model, X)
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{ComponentReactivePowerVariableLimitsConstraint},
    U::Type{<:VariableType},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation}
    time_steps = get_time_steps(container)
    var = get_variable(container, ComponentReactivePowerVariable(), V)
    device_names = [PSY.get_name(d) for d in devices]
    subcomp_types = get_subcomponent_types(U)

    constraint_ub = add_constraints_container!(
        container,
        ReactiveRangeConstraint(),
        V,
        device_names,
        subcomp_types,
        time_steps;
        meta="ub",
        sparse=true,
    )
    constraint_lb = add_constraints_container!(
        container,
        ReactiveRangeConstraint(),
        V,
        device_names,
        subcomp_types,
        time_steps;
        meta="lb",
        sparse=true,
    )

    for t in time_steps, d in devices, subcomp in subcomp_types
        !does_subcomponent_exist(d, subcomp) && continue
        name = PSY.get_name(d)
        limits = get_min_max_limits(d, subcomp, T, W)
        subcomp_key = string(subcomp)
        constraint_ub[name, subcomp, t] =
            JuMP.@constraint(container.JuMPmodel, var[name, subcomp_key, t] <= limits.max)
        constraint_lb[name, subcomp, t] =
            JuMP.@constraint(container.JuMPmodel, var[name, subcomp_key, t] >= limits.min)
    end
    return
end
######################## Energy balance constraints ############################

function add_constraints!(
    container::OptimizationContainer,
    ::Type{EnergyBalanceConstraint},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {V <: PSY.HybridSystem, W <: AbstractHybridFormulation, X <: PM.AbstractPowerModel}
    time_steps = get_time_steps(container)
    resolution = get_resolution(container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / MINUTES_IN_HOUR
    names = [PSY.get_name(x) for x in devices if !isnothing(PSY.get_storage(x))]
    initial_conditions = get_initial_condition(container, InitialEnergyLevel(), V)
    energy_var = get_variable(container, EnergyVariable(), V)
    powerin_var = get_variable(container, ActivePowerInVariable(), V)
    powerout_var = get_variable(container, ActivePowerOutVariable(), V)

    constraint = add_constraints_container!(
        container,
        EnergyBalanceConstraint(),
        V,
        names,
        time_steps,
    )

    for ic in initial_conditions
        device = get_component(ic)
        isnothing(PSY.get_storage(device)) && continue
        storage_device = PSY.get_storage(device)
        efficiency = PSY.get_efficiency(storage_device)
        name = PSY.get_name(device)
        constraint[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            energy_var[name, 1] ==
            get_value(ic) +
            (
                powerin_var[name, 1] * efficiency.in -
                (powerout_var[name, 1] / efficiency.out)
            ) * fraction_of_hour
        )

        for t in time_steps[2:end]
            constraint[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                energy_var[name, t] ==
                energy_var[name, t - 1] +
                (
                    powerin_var[name, t] * efficiency.in -
                    (powerout_var[name, t] / efficiency.out)
                ) * fraction_of_hour
            )
        end
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{PowerOutputRangeConstraint},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, D},
    ::Type{X},
) where {V <: PSY.HybridSystem, D <: AbstractHybridFormulation, X <: PM.AbstractPowerModel}
    time_steps = get_time_steps(container)
    name_index = [PSY.get_name(d) for d in devices]

    var_p = get_variable(container, ActivePowerVariable(), V)
    var_sub_p = get_variable(container, ComponentActivePowerVariable(), V)
    var_out = get_variable(container, ActivePowerOutVariable(), V)
    var_in = get_variable(container, ActivePowerInVariable(), V)

    constraint = add_constraints_container!(
        container,
        PowerOutputRangeConstraint(),
        V,
        name_index,
        time_steps,
    )

    for d in devices, t in time_steps
        name = PSY.get_name(d)

        constraint[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            var_p[name, t] ==
            var_sub_p[name, string(PSY.RenewableGen), t] +
            var_sub_p[name, string(PSY.ThermalGen), t] +
            var_out[name, t] - var_in[name, t]
        )
    end

    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{ReactivePowerConstraint},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, D},
    ::Type{X},
) where {V <: PSY.HybridSystem, D <: AbstractHybridFormulation, X <: PM.AbstractPowerModel}
    time_steps = get_time_steps(container)
    name_index = [PSY.get_name(d) for d in devices]

    var_q = get_variable(container, ReactivePowerVariable(), V)
    var_sub_q = get_variable(container, ComponentReactivePowerVariable(), V)

    constraint = add_constraints_container!(
        container,
        ReactivePowerConstraint(),
        V,
        name_index,
        time_steps,
    )

    for d in devices, t in time_steps
        name = PSY.get_name(d)

        constraint[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            var_q[name, t] ==
            var_sub_q[name, string(PSY.RenewableGen), t] +
            var_sub_q[name, string(PSY.ThermalGen), t] +
            var_sub_q[name, string(PSY.Storage), t]
        )
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{InterConnectionLimitConstraint},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, D},
    ::Type{X},
) where {V <: PSY.HybridSystem, D <: AbstractHybridFormulation, X <: PM.AbstractPowerModel}
    time_steps = get_time_steps(container)
    name_index = [PSY.get_name(d) for d in devices]

    var_q = get_variable(container, ReactivePowerVariable(), V)
    var_p = get_variable(container, ActivePowerVariable(), V)

    constraint = add_constraints_container!(
        container,
        InterConnectionLimitConstraint(),
        V,
        name_index,
        time_steps,
    )

    for d in devices, t in time_steps
        name = PSY.get_name(d)
        rating = PSY.get_interconnection_rating(d)
        constraint[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            rating^2 == var_q[name, t]^2 + var_p[name, t]^2
        )
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{ReserveEnergyConstraint},
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.HybridSystem, D <: AbstractHybridFormulation}
    time_steps = get_time_steps(container)
    var_e = get_variable(container, EnergyVariable(), T)
    r_up = get_variable(container, ComponentActivePowerReserveUpVariable(), T)
    r_dn = get_variable(container, ComponentActivePowerReserveDownVariable(), T)
    names = [PSY.get_name(x) for x in devices if does_subcomponent_exist(d, PSY.Storage)]
    con_up = add_constraints_container!(
        container,
        ReserveEnergyConstraint(),
        T,
        names,
        time_steps,
        meta="up",
    )
    con_dn = add_constraints_container!(
        container,
        ReserveEnergyConstraint(),
        T,
        names,
        time_steps,
        meta="dn",
    )

    for d in devices, t in time_steps
        !does_subcomponent_exist(d, PSY.Storage) && continue
        name = PSY.get_name(d)
        limits = PSY.get_state_of_charge_limits(PSY.get_storage(d))
        efficiency = PSY.get_efficiency(d)
        con_up[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            r_up[name, t] <= (var_e[name, t] - limits.min) * efficiency.out
        )
        con_dn[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            r_dn[name, t] <= (limits.max - var_e[name, t]) / efficiency.in
        )
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{RangeLimitConstraint},
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.HybridSystem, D <: AbstractHybridFormulation}
    time_steps = get_time_steps(container)
    var_in = get_variable(container, ActivePowerInVariable(), T)
    var_out = get_variable(container, ActivePowerOutVariable(), T)
    r_up = get_variable(container, ComponentActivePowerReserveUpVariable(), T)
    r_dn = get_variable(container, ComponentActivePowerReserveDownVariable(), T)
    names = [PSY.get_name(x) for x in devices]
    con_up = add_constraints_container!(
        container,
        RangeLimitConstraint(),
        T,
        names,
        time_steps,
        meta="up",
    )
    con_dn = add_constraints_container!(
        container,
        RangeLimitConstraint(),
        T,
        names,
        time_steps,
        meta="dn",
    )

    for d in devices, t in time_steps
        name = PSY.get_name(d)
        out_limits = PSY.get_output_active_power_limits(d)
        in_limits = PSY.get_input_active_power_limits(d)
        con_up[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            r_up[name, t] <= var_in[name, t] + (out_limits.max - var_out[name, t])
        )
        con_dn[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            r_dn[name, t] <= var_out[name, t] + (in_limits.max - var_in[name, t])
        )
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{ComponentReserveUpBalance},
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.HybridSystem, D <: AbstractHybridFormulation}
    time_steps = get_time_steps(container)
    sub_r_up = get_variable(container, ComponentActivePowerReserveUpVariable(), T)
    sub_expr_up = get_expression(container, ComponentReserveUpBalanceExpression(), T)
    names = [PSY.get_name(x) for x in devices]
    con_up = add_constraints_container!(
        container,
        ComponentReserveUpBalance(),
        T,
        names,
        time_steps,
    )

    for d in devices, t in time_steps
        name = PSY.get_name(d)
        con_up[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            sub_expr_up[name, t] == sum(
                sub_r_up[name, string(sub_comp_type), t] for
                sub_comp_type in [PSY.ThermalGen, PSY.RenewableGen, PSY.Storage]
            )
        )
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{ComponentReserveDownBalance},
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.HybridSystem, D <: AbstractHybridFormulation}
    time_steps = get_time_steps(container)
    sub_r_dn = get_variable(container, ComponentActivePowerReserveDownVariable(), T)
    sub_expr_dn = get_expression(container, ComponentReserveDownBalanceExpression(), T)
    names = [PSY.get_name(x) for x in devices]
    con_dn = add_constraints_container!(
        container,
        ComponentReserveDownBalance(),
        T,
        names,
        time_steps,
    )

    for d in devices, t in time_steps
        name = PSY.get_name(d)
        con_dn[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            sub_expr_dn[name, t] == sum(
                sub_r_dn[name, string(sub_comp_type), t] for
                sub_comp_type in [PSY.ThermalGen, PSY.RenewableGen, PSY.Storage]
            )
        )
    end
    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{D},
    formulation::AbstractHybridFormulation,
) where {D <: PSY.HybridSystem}
    add_initial_condition!(container, devices, formulation, InitialEnergyLevel())
    return
end

########################### Cost Function Calls#############################################
function objective_function!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.HybridSystem, U <: AbstractHybridFormulation}
    add_variable_cost!(container, ActivePowerVariable(), devices, U())
    add_proportional_cost!(container, OnVariable(), devices, U())
    return
end


text:
"""
    galerkin_tensor(
        B::AbstractBSplineBasis,
        (D₁::Derivative, D₂::Derivative, D₃::Derivative),
        [T = Float64],
    )

Compute 3D banded tensor appearing from quadratic terms in Galerkin method.

As with [`galerkin_matrix`](@ref), it is also possible to combine different
functional bases by passing, instead of `B`, a tuple `(B₁, B₂, B₃)` of three
`AbstractBSplineBasis`.
For now, the first two bases, `B₁` and `B₂`, must have the same length.

The tensor is efficiently stored in a [`BandedTensor3D`](@ref) object.
"""
function galerkin_tensor end

function galerkin_tensor(
        Bs::NTuple{3,AbstractBSplineBasis},
        deriv::DerivativeCombination{3},
        ::Type{T} = Float64,
    ) where {T}
    _check_bases(Bs)
    dims = length.(Bs)
    b = order(first(Bs)) - 1   # band width
    if length(Bs[1]) != length(Bs[2])
        throw(DimensionMismatch("the first two bases must have the same lengths"))
    end
    δl = first(num_constraints(Bs[3])) - first(num_constraints(Bs[1]))
    A = BandedTensor3D{T}(undef, dims, Val(b), bandshift=(0, 0, δl))
    galerkin_tensor!(A, Bs, deriv)
end

galerkin_tensor(B::AbstractBSplineBasis, args...) =
    galerkin_tensor((B, B, B), args...)

"""
    galerkin_tensor!(
        A::BandedTensor3D,
        B::AbstractBSplineBasis,
        (D₁::Derivative, D₂::Derivative, D₃::Derivative),
    )

Compute 3D Galerkin tensor in-place.

See [`galerkin_tensor`](@ref) for details.
"""
function galerkin_tensor! end

function galerkin_tensor!(
        A::BandedTensor3D, Bs::NTuple{3,AbstractBSplineBasis},
        deriv::DerivativeCombination{3},
    )
    _check_bases(Bs)

    Ns = size(A)
    if Ns != length.(Bs)
        throw(DimensionMismatch("wrong dimensions of Galerkin tensor"))
    end

    Bi, Bj, Bl = Bs
    Ni, Nj, Nl = Ns
    @assert Ni == Nj  # verified at construction of BandedTensor3D

    ioff = first(num_constraints(Bi))
    joff = first(num_constraints(Bj))
    loff = first(num_constraints(Bl))

    # Orders and knots are assumed to be the same (see _check_bases).
    k = order(Bi)
    ts = knots(Bi)

    if bandwidth(A) != k - 1
        throw(ArgumentError("BandedTensor3D must have bandwidth = $(k - 1)"))
    end

    δl, _ = num_constraints(Bl) .- num_constraints(Bi)

    if bandshift(A) != (0, 0, δl)
        throw(ArgumentError("BandedTensor3D must have bandshift = (0, 0, $δl)"))
    end

    same_12 = Bi == Bj && deriv[1] == deriv[2]
    same_13 = Bi == Bl && deriv[1] == deriv[3]
    same_23 = Bj == Bl && deriv[2] == deriv[3]

    # Quadrature information (weights, nodes).
    quadx, quadw = _quadrature_prod(Val(3k - 3))

    fill!(A, 0)
    T = eltype(A)
    Al = zero(MMatrix{2k - 1, 2k - 1, T})
    nlast = last(eachindex(ts))

    @inbounds for n in eachindex(ts)
        n == nlast && break
        tn, tn1 = ts[n], ts[n + 1]
        tn1 == tn && continue  # interval of length = 0

        metric = QuadratureMetric(tn, tn1)
        xs = metric .* quadx

        ilast = n - ioff
        jlast = n - joff
        llast = n - loff
        l₀ = llast - order(Bs[3])

        for (x, w) ∈ zip(xs, quadw)
            _, bis = evaluate_all(Bi, x, deriv[1], T; ileft = ilast)
            _, bjs = same_12 ?
                (ilast, bis) : evaluate_all(Bj, x, deriv[2], T; ileft = jlast)
            _, bls = same_13 ?
                (ilast, bis) : same_23 ?
                (jlast, bjs) : evaluate_all(Bl, x, deriv[3], T; ileft = llast)

            # Iterate in increasing order (not sure if it really helps performance)
            bls = reverse(bls)

            # We compute the submatrix A[:, :, l] for each `l`.
            for (δl, bl) in pairs(bls)
                iszero(bl) && continue  # can prevent problems at the borders
                l = l₀ + δl
                y₀ = metric.α * w * bl
                @inbounds fill!(Al, 0)
                inds = BandedTensors.band_indices(A, l)
                # @assert inds == (l - k + 1:l + k - 1) .+ bandshift(A)[3]
                # @assert (ilast - k + 1:ilast) ⊆ inds
                # @assert (jlast - k + 1:jlast) ⊆ inds
                Δi = first(inds) - 1
                Δj = Δi
                for (δj, bj) ∈ pairs(bjs)
                    y₁ = y₀ * bj
                    j = jlast + 1 - δj
                    jj = j - Δj
                    for (δi, bi) ∈ pairs(bis)
                        i = ilast + 1 - δi  # actual index of basis function
                        ii = i - Δi         # index in submatrix Al
                        @inbounds Al[ii, jj] = y₁ * bi
                    end
                end
                @inbounds A[:, :, l] += Al
            end
        end
    end

    A
end

ext: jl

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
"""
    Link

A WSTP link object.

# External links
 - [WSLINK](https://reference.wolfram.com/language/ref/c/WSLINK.html)
"""
mutable struct Link
    ptr::Ptr{Cvoid}
end

const CLink = Ptr{Cvoid}
Base.cconvert(::Type{CLink}, link::Link) = link.ptr

# default link to use if none other is provided
# this will be initialized lazily
const defaultlink = Link(C_NULL)

# we define an in-place version to operate on defaultlink
function open!(link::Link, args::AbstractString)
    # MLOpenString
    # local link
    err = Ref{Cint}()
    ptr = ccall((:MLOpenString, mlib), CLink,
                 (Env, Cstring, Ptr{Cint}),
                 env, args, err)
    if err[] != 0
        error("Could not open MathLink link")
    end
    link.ptr = ptr

    refcount_inc()
    finalizer(close, link)
    return link
end

open(args::AbstractString) = open!(Link(C_NULL), args)

function Base.close(link::Link)
    if link.ptr != C_NULL
        ccall((:MLClose, mlib), Cvoid, (CLink,), link)
        link.ptr = C_NULL
        refcount_dec()
    end
end

function _defaultlink()
    if defaultlink.ptr == C_NULL
        args = "-linkname '\"$mker\" -mathlink' -linkmode launch"
        open!(defaultlink, args)

        # Ignore first input packet
        @assert nextpacket(defaultlink) == PKT_INPUTNAME
        newpacket(defaultlink)
    end
    return defaultlink
end

mutable struct Mark
    link::Link
    ptr::Ptr{Cvoid}
end
const CMark = Ptr{Cvoid}
Base.cconvert(::Type{CMark}, mark::Mark) = mark.ptr

# we don't overload Base.mark as it has different behaviour
function Mark(link::Link)
    # WSMARK WSCreateMark(WSLINK link)
    ptr = ccall((:MLCreateMark, mlib), Ptr{Cvoid}, (CLink,), link)
    ptr == C_NULL && throw(MathLinkError(link))
    mark = Mark(link, ptr)
    refcount_inc()
    finalizer(close, mark)
    mark
end

function Base.seek(link::Link, mark::Mark, offset::Integer=Cint(0))
    # WSMARK WSSeekToMark(WSLINK link,WSMARK mark,int n)
    ptr = ccall((:MLSeekToMark, mlib), Ptr{Cvoid},
                (CLink, CMark, Cint), link, mark, offset)
    ptr == C_NULL && throw(MathLinkError(link))
    return nothing
end

function Base.close(mark::Mark)
    if mark.ptr != C_NULL
        # void WSDestroyMark(WSLINK link,WSMARK mark)
        ccall((:MLDestroyMark, mlib), Cvoid,
              (CLink, CMark), mark.link, mark)
        mark.ptr = C_NULL
        refcount_dec()
    end
end


text:
using Test

using Presentation

if get(ENV, "CI", "") == ""

    include("test_presentation.jl")

end


text:
module AccuracyBenchmark

using DataFrames
using Distributions
import FITSIO
import StaticArrays
import WCS

import ..Config
import ..DeterministicVI
import ..Infer
import ..Model
import ..ParallelRun
import ..SDSSIO
using ..Coordinates: angular_separation, match_coordinates

const ARCSEC_PER_DEGREE = 3600
const SDSS_ARCSEC_PER_PIXEL = 0.396
const SDSS_DATA_DIR = joinpath(Pkg.dir("Celeste"), "test", "data")
const STRIPE82_RCF = SDSSIO.RunCamcolField(4263, 5, 119)

struct BenchmarkFitsFileNotFound <: Exception
    filename::String
end

struct MatchException <: Exception
    msg::String
end

struct MissingColumnsError <: Exception
    missing_columns::Vector{Symbol}
end

################################################################################
# Load/store a catalog DF from/to a CSV.
################################################################################

CATALOG_COLUMNS = Set([
    :right_ascension_deg,
    :declination_deg,
    :is_star,
    :reference_band_flux_nmgy,
    :color_log_ratio_ug,
    :color_log_ratio_gr,
    :color_log_ratio_ri,
    :color_log_ratio_iz,
    :de_vaucouleurs_mixture_weight,
    :minor_major_axis_ratio,
    :half_light_radius_px,
    :angle_deg,
    :is_saturated,
])

STDERR_COLUMNS = Set([
    :log_reference_band_flux_stderr,
    :color_log_ratio_ug_stderr,
    :color_log_ratio_gr_stderr,
    :color_log_ratio_ri_stderr,
    :color_log_ratio_iz_stderr,
])

function assert_columns_are_present(catalog_df::DataFrame, required_columns::Set{Symbol})
    missing_columns = setdiff(required_columns, Set(names(catalog_df)))
    if !isempty(missing_columns)
        throw(MissingColumnsError([missing_columns...]))
    end
end

function read_catalog(csv_file::String)
    @printf("Reading '%s'...\n", csv_file)
    catalog_df = readtable(csv_file)
    assert_columns_are_present(catalog_df, CATALOG_COLUMNS)
    catalog_df
end

function write_catalog(csv_file::String, catalog_df::DataFrame)
    assert_columns_are_present(catalog_df, CATALOG_COLUMNS)
    @printf("Writing '%s'...\n", csv_file)
    writetable(csv_file, catalog_df)
end

function append_hash_to_file(filename::String)
    contents_hash = open(filename) do stream
        hash(read(stream))
    end
    hash_string = hex(contents_hash)[1:10]
    base, extension = splitext(filename)
    new_filename = @sprintf("%s_%s%s", base, hash_string, extension)
    @printf("Renaming %s -> %s\n", filename, new_filename)
    mv(filename, new_filename, remove_destination=true)
end

################################################################################
# Read various catalogs to a common catalog DF format
################################################################################

"""
Convert between SDSS asinh mags and SDSS flux (nMgy).
See http://www.sdss.org/dr12/algorithms/magnitudes/#asinh.
"""
const ASINH_SOFTENING_PARAMETERS = [
    1.4e-10, # for band 1 = u
    0.9e-10,
    1.2e-10,
    1.8e-10,
    7.4e-10, # for band 5 = z
]

function mag_to_flux(mags::AbstractFloat, band_index::Int)
    b = ASINH_SOFTENING_PARAMETERS[band_index]
    1e9 * 2 * b * sinh(-log(10) / 2.5 * mags - log(b))
end

function flux_to_mag(flux_nmgy::AbstractFloat, band_index::Int)
    b = ASINH_SOFTENING_PARAMETERS[band_index]
    -2.5 / log(10) * (asinh(flux_nmgy * 1e-9 / (2*b)) + log(b))
end

function color_from_fluxes(flux1::AbstractFloat, flux2::AbstractFloat)
    if flux1 <= 0 || flux2 <= 0
        NA
    else
        log(flux2 / flux1)
    end
end

function color_from_mags(mags1::AbstractFloat, band1::Int, mags2::AbstractFloat, band2::Int)
    color_from_fluxes(mag_to_flux(mags1, band1), mag_to_flux(mags2, band2))
end

canonical_angle(angle_deg) = angle_deg - floor(angle_deg / 180) * 180

STRIPE_82_CATALOG_KEYS = [
    :objid, :rerun, :run, :camcol, :field, :flags,
    :ra, :dec, :probpsf,
    :psfmag_u, :psfmag_g, :psfmag_r, :psfmag_i, :psfmag_z,
    :devmag_u, :devmag_g, :devmag_r, :devmag_i, :devmag_z,
    :expmag_u, :expmag_g, :expmag_r, :expmag_i, :expmag_z,
    :fracdev_r,
    :devab_r, :expab_r,
    :devphi_r, :expphi_r,
    :devrad_r, :exprad_r,
    :flags, :is_saturated,
]

function load_fits_table_as_data_frame(filename, extension_index, field_names)
    result_df = DataFrame()
    fits = FITSIO.FITS(filename)
    try
        for key in field_names
            result_df[key] = read(fits[extension_index], string(key))
        end
        return result_df
    finally
        close(fits)
    end
end

function load_stripe82_fits_catalog_as_data_frame(filename, extension_index)
    load_fits_table_as_data_frame(filename, extension_index, STRIPE_82_CATALOG_KEYS)
end

"""
Load Stripe 82 objects into a DataFrame. `fits_filename` should be a FITS file created by running a
CasJobs (skyserver.sdss.org/casjobs/) query on the Stripe82 database. See
https://github.com/jeff-regier/Celeste.jl/wiki/About-SDSS-and-Stripe-82.
"""
function load_coadd_catalog(fits_filename)
    raw_df = load_stripe82_fits_catalog_as_data_frame(fits_filename, 2)

    usedev = raw_df[:fracdev_r] .> 0.5  # true=> use dev, false=> use exp
    dev_or_exp(dev_column, exp_column) = ifelse.(usedev, raw_df[dev_column], raw_df[exp_column])
    is_star = [x != 0 for x in raw_df[:probpsf]]
    function star_or_galaxy(star_column, galaxy_dev_column, galaxy_exp_column)
        ifelse.(is_star, raw_df[star_column], dev_or_exp(galaxy_dev_column, galaxy_exp_column))
    end

    mag_u = star_or_galaxy(:psfmag_u, :devmag_u, :expmag_u)
    mag_g = star_or_galaxy(:psfmag_g, :devmag_g, :expmag_g)
    mag_r = star_or_galaxy(:psfmag_r, :devmag_r, :expmag_r)
    mag_i = star_or_galaxy(:psfmag_i, :devmag_i, :expmag_i)
    mag_z = star_or_galaxy(:psfmag_z, :devmag_z, :expmag_z)

    result = DataFrame()
    result[:objid] = [string(objid) for objid in raw_df[:objid]]
    result[:right_ascension_deg] = raw_df[:ra]
    result[:declination_deg] = raw_df[:dec]
    result[:is_saturated] = (raw_df[:is_saturated] .!= 0)
    result[:is_star] = is_star

    result[:reference_band_flux_nmgy] = mag_to_flux.(mag_r, 3)

    result[:color_log_ratio_ug] = color_from_mags.(mag_u, 1, mag_g, 2)
    result[:color_log_ratio_gr] = color_from_mags.(mag_g, 2, mag_r, 3)
    result[:color_log_ratio_ri] = color_from_mags.(mag_r, 3, mag_i, 4)
    result[:color_log_ratio_iz] = color_from_mags.(mag_i, 4, mag_z, 5)

    # gal shape -- fracdev
    result[:de_vaucouleurs_mixture_weight] = raw_df[:fracdev_r]

    # Note that the SDSS photo pipeline doesn't constrain the de Vaucouleur
    # profile parameters and exponential disk parameters (A/B, angle, scale)
    # to be the same, whereas Celeste does. Here, we pick one or the other
    # from SDSS, based on fracdev - we'll get the parameters corresponding
    # to the dominant component. Later, we limit comparison to objects with
    # fracdev close to 0 or 1 to ensure that we're comparing apples to apples.

    result[:minor_major_axis_ratio] = dev_or_exp(:devab_r, :expab_r)

    # gal effective radius (re)
    re_arcsec = dev_or_exp(:devrad_r, :exprad_r)
    re_pixel = re_arcsec ./ SDSS_ARCSEC_PER_PIXEL
    result[:half_light_radius_px] = convert(Vector{Float64}, re_pixel)

    # gal angle (degrees)
    raw_phi = dev_or_exp(:devphi_r, :expphi_r)
    result[:angle_deg] = canonical_angle.(raw_phi)

    return result
end

function object_dict_to_data_frame(objects::Dict)
    result_df = DataFrame()
    for (key, values) in objects
        result_df[Symbol(key)] = values
    end
    result_df
end

"""
Load the SDSS photoObj catalog used to initialize celeste, and reformat column
names to match what the rest of the scoring code expects.
"""
function load_primary(rcf::SDSSIO.RunCamcolField, stagedir::String)
    strategy = SDSSIO.PlainFITSStrategy(stagedir)
    raw_df = object_dict_to_data_frame(SDSSIO.read_photoobj(strategy, rcf))

    usedev = raw_df[:frac_dev] .> 0.5  # true=> use dev, false=> use exp
    dev_or_exp(dev_column, exp_column) = ifelse.(usedev, raw_df[dev_column], raw_df[exp_column])
    function star_or_galaxy(star_column, galaxy_dev_column, galaxy_exp_column)
        ifelse.(raw_df[:is_star], raw_df[star_column], dev_or_exp(galaxy_dev_column, galaxy_exp_column))
    end

    flux_u = star_or_galaxy(:psfflux_u, :devflux_u, :expflux_u)
    flux_g = star_or_galaxy(:psfflux_g, :devflux_g, :expflux_g)
    flux_r = star_or_galaxy(:psfflux_r, :devflux_r, :expflux_r)
    flux_i = star_or_galaxy(:psfflux_i, :devflux_i, :expflux_i)
    flux_z = star_or_galaxy(:psfflux_z, :devflux_z, :expflux_z)

    result = DataFrame()
    result[:objid] = raw_df[:objid]
    result[:right_ascension_deg] = raw_df[:ra]
    result[:declination_deg] = raw_df[:dec]
    result[:is_star] = raw_df[:is_star]

    result[:reference_band_flux_nmgy] = flux_r

    result[:color_log_ratio_ug] = color_from_fluxes.(flux_u, flux_g)
    result[:color_log_ratio_gr] = color_from_fluxes.(flux_g, flux_r)
    result[:color_log_ratio_ri] = color_from_fluxes.(flux_r, flux_i)
    result[:color_log_ratio_iz] = color_from_fluxes.(flux_i, flux_z)

    result[:de_vaucouleurs_mixture_weight] = raw_df[:frac_dev]
    result[:minor_major_axis_ratio] = dev_or_exp(:ab_dev, :ab_exp)

    # gal effective radius (re)
    re_arcsec = dev_or_exp(:theta_dev, :theta_exp)
    re_pixel = re_arcsec ./ SDSS_ARCSEC_PER_PIXEL
    result[:half_light_radius_px] = convert(Vector{Float64}, re_pixel)

    # gal angle (degrees)
    raw_phi = dev_or_exp(:phi_dev, :phi_exp)
    result[:angle_deg] = canonical_angle.(raw_phi)

    # primary is better at flagging oversaturated sources than coadd
    result[:is_saturated] = flux_to_mag.(raw_df[:psfflux_r], 3) .< 16

    return result
end

function fluxes_from_colors(reference_band_flux_nmgy::Float64, color_log_ratios::DataVector{Float64})
    @assert length(color_log_ratios) == 4
    color_ratios = exp.(color_log_ratios)
    fluxes = DataArray(Float64, 5)
    fluxes[3] = reference_band_flux_nmgy
    fluxes[4] = fluxes[3] * color_ratios[3]
    fluxes[5] = fluxes[4] * color_ratios[4]
    fluxes[2] = fluxes[3] / color_ratios[2]
    fluxes[1] = fluxes[2] / color_ratios[1]
    fluxes
end

function fluxes_from_colors(reference_band_flux_nmgy::Float64, color_log_ratios::Vector{Float64})
    fluxes = fluxes_from_colors(reference_band_flux_nmgy, DataArray{Float64}(color_log_ratios))
    convert(Vector{Float64}, fluxes)
end

function get_median_fluxes(variational_params::Vector{Float64}, source_type::Int64)
    fluxes_from_colors(
        exp(variational_params[Model.ids.flux_loc[source_type]]),
        variational_params[Model.ids.color_mean[:, source_type]],
    )
end

function variational_parameters_to_data_frame_row(variational_params::Vector{Float64})
    ids = Model.ids
    result = DataFrame()
    result[:right_ascension_deg] = variational_params[ids.pos[1]]
    result[:declination_deg] = variational_params[ids.pos[2]]
    result[:is_saturated] = false
    result[:is_star] = variational_params[ids.is_star[1, 1]]
    result[:de_vaucouleurs_mixture_weight] = variational_params[ids.gal_fracdev]
    result[:minor_major_axis_ratio] = variational_params[ids.gal_ab]
    result[:half_light_radius_px] = (
        variational_params[ids.gal_scale] * sqrt(variational_params[ids.gal_ab])
    )
    result[:angle_deg] = canonical_angle(180 / pi * variational_params[ids.gal_angle])

    star_galaxy_index = (result[1, :is_star] > 0.5 ? 1 : 2)
    fluxes = get_median_fluxes(variational_params, star_galaxy_index)
    result[:reference_band_flux_nmgy] = fluxes[3]
    result[:color_log_ratio_ug] = color_from_fluxes(fluxes[1], fluxes[2])
    result[:color_log_ratio_gr] = color_from_fluxes(fluxes[2], fluxes[3])
    result[:color_log_ratio_ri] = color_from_fluxes(fluxes[3], fluxes[4])
    result[:color_log_ratio_iz] = color_from_fluxes(fluxes[4], fluxes[5])

    result[:log_reference_band_flux_stderr] = sqrt(variational_params[ids.flux_scale[star_galaxy_index]])
    result[:color_log_ratio_ug_stderr] = sqrt(variational_params[ids.color_var[1, star_galaxy_index]])
    result[:color_log_ratio_gr_stderr] = sqrt(variational_params[ids.color_var[2, star_galaxy_index]])
    result[:color_log_ratio_ri_stderr] = sqrt(variational_params[ids.color_var[3, star_galaxy_index]])
    result[:color_log_ratio_iz_stderr] = sqrt(variational_params[ids.color_var[4, star_galaxy_index]])

    result
end

"""
Convert Celeste results to a dataframe.
"""
function celeste_to_df(results::Vector{ParallelRun.OptimizedSource})
    rows = [variational_parameters_to_data_frame_row(result.vs) for result in results]
    vcat(rows...)
end

################################################################################
# Generate a random catalog from the Celeste prior
################################################################################

# This is the ratio of stars derived from the catalog used the generate the prior; the value 0.99 is
# currently used in inference to account for the extra flexibility of the galaxy model
const PRIOR_PROBABILITY_OF_STAR = 0.28

"""
Draw a random source from the Celeste prior, returning a 1-row data frame.
"""
function draw_source_params(prior)
    is_star = (rand(Bernoulli(PRIOR_PROBABILITY_OF_STAR)) == 1)
    source_type_index = is_star ? 1 : 2
    reference_band_flux_nmgy = exp(rand(
        Normal(prior.flux_mean[source_type_index], sqrt(prior.flux_var[source_type_index]))
    ))

    color_mixture_weights = prior.k[:, source_type_index]
    NUM_COLOR_COMPONENTS = length(color_mixture_weights)
    color_components = MvNormal[
        MvNormal(prior.color_mean[:, k, source_type_index], prior.color_cov[:, :, k, source_type_index])
        for k in 1:NUM_COLOR_COMPONENTS
    ]
    color_log_ratios = rand(MixtureModel(color_components, color_mixture_weights))

    if !is_star
        half_light_radius_px = exp(rand(
            Normal(prior.gal_scale_mean, sqrt(prior.gal_scale_var))
        ))
        angle_deg = rand(Uniform(0, 180))
        minor_major_axis_ratio = rand(Beta(2, 2))
        de_vaucouleurs_mixture_weight = rand(Beta(0.5, 0.5))
    else
        half_light_radius_px = -1
        angle_deg = -1
        minor_major_axis_ratio = -1
        de_vaucouleurs_mixture_weight = -1
    end

    # Use approximate size of SDSS field in degrees
    right_ascension_deg = rand(Uniform(0, 0.14))
    declination_deg = rand(Uniform(0, 0.22))

    DataFrame(
        right_ascension_deg=right_ascension_deg,
        declination_deg=declination_deg,
        is_star=is_star,
        reference_band_flux_nmgy=reference_band_flux_nmgy,
        color_log_ratio_ug=color_log_ratios[1],
        color_log_ratio_gr=color_log_ratios[2],
        color_log_ratio_ri=color_log_ratios[3],
        color_log_ratio_iz=color_log_ratios[4],
        de_vaucouleurs_mixture_weight=de_vaucouleurs_mixture_weight,
        minor_major_axis_ratio=minor_major_axis_ratio,
        half_light_radius_px=half_light_radius_px,
        angle_deg=angle_deg,
        is_saturated=false,
    )
end

"""
Draw sources at random from Celeste prior, returning a catalog DF.
"""
function generate_catalog_from_celeste_prior(num_sources::Int64, seed::Int64)
    srand(seed)
    prior = Model.load_prior()
    vcat(
        [draw_source_params(prior) for index in 1:num_sources]...
    )
end

################################################################################
# Support for running Celeste on test imagery
################################################################################

## Load a multi-extension FITS imagery file

struct FitsImage
    pixels::Matrix{Float32}
    header::FITSIO.FITSHeader
    wcs::WCS.WCSTransform
end

function read_fits(filename::String)
    println("Reading '$filename'...")
    if !isfile(filename)
        throw(BenchmarkFitsFileNotFound(filename))
    end

    fits = FITSIO.FITS(filename)
    try
        println("Found $(length(fits)) extensions.")
        map(fits) do extension
            pixels = read(extension)
            header = FITSIO.read_header(extension)
            wcs = WCS.from_header(FITSIO.read_header(fits[1], String))[1]
            FitsImage(pixels, header, wcs)
        end
    finally
        close(fits)
    end
end

## Create Celeste images from FITS imagery

function make_simple_psf(psf_sigma_px)
    alphaBar = [1.; 0.]
    xiBar = [0.; 0.]
    tauBar = [psf_sigma_px^2 0.; 0. psf_sigma_px^2]
    [
        Model.PsfComponent(
            alphaBar[k],
            StaticArrays.SVector{2, Float64}(xiBar),
            StaticArrays.SMatrix{2, 2, Float64, 4}(tauBar)
        )
        for k in 1:2
    ]
end

function make_psf_from_header(header::FITSIO.FITSHeader)
    if haskey(header, "CLSIGMA")
        return make_simple_psf(header["CLSIGMA"])
    end

    num_components = header["PSFNCOMP"]
    map(1:num_components) do index
        get(field) = header["PSF$(index)_$(field)"]
        alpha_bar = get("A")
        xi_bar = Float64[get("X1"), get("X2")]
        tau_bar = Float64[get("T11") get("T12"); get("T12") get("T22")]
        Model.PsfComponent(
            alpha_bar,
            StaticArrays.SVector{2, Float64}(xi_bar),
            StaticArrays.SMatrix{2, 2, Float64, 4}(tau_bar),
        )
    end
end

function serialize_psf_to_header(psf::Vector{Model.PsfComponent}, header::FITSIO.FITSHeader)
    header["PSFNCOMP"] = length(psf)
    for (index, component) in enumerate(psf)
        @assert component.tauBar[1, 2] == component.tauBar[2, 1]
        function set(field, value)
            header["PSF$(index)_$(field)"] = value
        end
        set("A", component.alphaBar)
        set("X1", component.xiBar[1])
        set("X2", component.xiBar[2])
        set("T11", component.tauBar[1, 1])
        set("T12", component.tauBar[1, 2])
        set("T22", component.tauBar[2, 2])
    end
end

function make_image(
    pixels::Matrix{Float32}, band_index::Int, wcs::WCS.WCSTransform, psf::Vector{Model.PsfComponent},
    sky_level_nmgy::Float64, nelec_per_nmgy::Float64
)
    height_px, width_px = size(pixels)
    sky_intensity = Model.SkyIntensity(
        fill(sky_level_nmgy, height_px, width_px),
        collect(1:height_px),
        collect(1:width_px),
        ones(height_px),
    )
    Model.Image(
        height_px,
        width_px,
        pixels,
        band_index,
        wcs,
        psf,
        0, 0, 0, # run, camcol, field
        sky_intensity,
        fill(nelec_per_nmgy, height_px),
        Model.RawPSF(Matrix{Float64}(0, 0), 0, 0, Array{Float64,3}(0, 0, 0)),
    )
end

function make_images(band_extensions::Vector{FitsImage})
    map(enumerate(band_extensions)) do pair
        band_index, extension = pair
        make_image(
            extension.pixels,
            band_index,
            extension.wcs,
            make_psf_from_header(extension.header),
            convert(Float64, extension.header["CLSKY"]),
            convert(Float64, extension.header["CLIOTA"]),
        )
    end
end

## Create an initialization catalog for Celeste

function typical_band_fluxes(is_star::Bool)
    source_type_index = is_star ? 1 : 2
    prior_parameters::Model.PriorParams = Model.load_prior()
    # this is the mode. brightness is log normal.
    reference_band_flux = exp(
        prior_parameters.flux_mean[source_type_index] - prior_parameters.flux_var[source_type_index]
    )
    # Band relative intensities are a mixture of lognormals. Which mixture component has the most
    # weight?
    dominant_component = indmax(prior_parameters.k[:, source_type_index])
    # What are the most typical log relative intensities for that component?
    color_log_ratios = (
        prior_parameters.color_mean[:, dominant_component, source_type_index]
        - diag(prior_parameters.color_cov[:, :, dominant_component, source_type_index])
    )
    fluxes_from_colors(reference_band_flux, color_log_ratios)
end

function make_catalog_entry(x_position_world_coords::Float64,
                            y_position_world_coords::Float64)
    Model.CatalogEntry(
        [x_position_world_coords, y_position_world_coords],
        false, # is_star
        typical_band_fluxes(true),
        typical_band_fluxes(false),
        0.1, # gal_frac_dev
        0.7, # gal_ab
        pi / 4, # gal_angle
        4., # gal_scale
    )
end

ensure_small_flux(value) = (isna(value) || value <= 0) ? 1e-6 : value

na_to_default(value, default) = isna(value) ? default : value

function make_catalog_entry(row::DataFrameRow)
    color_log_ratios = DataArray{Float64}(DataArray(Any[
        row[:color_log_ratio_ug],
        row[:color_log_ratio_gr],
        row[:color_log_ratio_ri],
        row[:color_log_ratio_iz],
    ]))
    fluxes = fluxes_from_colors(row[:reference_band_flux_nmgy], color_log_ratios)
    fluxes = convert(Vector{Float64}, ensure_small_flux.(fluxes))
    minor_major_axis_ratio = na_to_default(row[:minor_major_axis_ratio], 0.8)
    Model.CatalogEntry(
        [row[:right_ascension_deg], row[:declination_deg]],
        row[:is_star] > 0.5,
        fluxes,
        fluxes,
        na_to_default(row[:de_vaucouleurs_mixture_weight], 0.5),
        minor_major_axis_ratio,
        na_to_default(row[:angle_deg], 0.) / 180.0 * pi,
        na_to_default(row[:half_light_radius_px] / sqrt(minor_major_axis_ratio), 2.)
        )
end

function make_initialization_catalog(catalog::DataFrame, use_full_initialzation::Bool)
    position_offset_width = SDSS_ARCSEC_PER_PIXEL / ARCSEC_PER_DEGREE # 1 pixel, in degrees
    map(eachrow(catalog)) do row
        if use_full_initialzation
            make_catalog_entry(row)
        else
            make_catalog_entry(
                row[:right_ascension_deg] + position_offset_width,
                row[:declination_deg] - 0.5 * position_offset_width
            )
        end
    end
end

################################################################################
# Support for generating imagery using Synthetic.jl
################################################################################

struct ImageGeometry
    height_px::Int64
    width_px::Int64
    world_coordinate_origin::Tuple{Float64, Float64} # (ra, dec)
    ra_degrees_per_pixel::Float64
    dec_degrees_per_pixel::Float64
end

function get_image_geometry(catalog_data::DataFrame; field_expand_arcsec=20.0)
    min_ra_deg = minimum(catalog_data[:right_ascension_deg])
    max_ra_deg = maximum(catalog_data[:right_ascension_deg])
    min_dec_deg = minimum(catalog_data[:declination_deg])
    max_dec_deg = maximum(catalog_data[:declination_deg])

    width_arcsec = (max_ra_deg - min_ra_deg) * ARCSEC_PER_DEGREE + 2 * field_expand_arcsec
    height_arcsec = (max_dec_deg - min_dec_deg) * ARCSEC_PER_DEGREE + 2 * field_expand_arcsec
    width_px = convert(Int64, round(width_arcsec / SDSS_ARCSEC_PER_PIXEL))
    height_px = convert(Int64, round(height_arcsec / SDSS_ARCSEC_PER_PIXEL))

    dec_degrees_per_pixel = SDSS_ARCSEC_PER_PIXEL / ARCSEC_PER_DEGREE
    ra_degrees_per_pixel = dec_degrees_per_pixel / cosd(min_dec_deg)

    ImageGeometry(
        height_px,
        width_px,
        (
            min_ra_deg - field_expand_arcsec / ARCSEC_PER_DEGREE,
            min_dec_deg - field_expand_arcsec / ARCSEC_PER_DEGREE,
        ),
        ra_degrees_per_pixel,
        dec_degrees_per_pixel,
    )
end

function make_template_images(
    catalog_data::DataFrame, psf_sigma_px::Float64, band_sky_level_nmgy::Vector{Float64},
    band_nelec_per_nmgy::Vector{Float64}
)
    @assert length(band_sky_level_nmgy) == 5
    @assert length(band_nelec_per_nmgy) == 5
    geometry = get_image_geometry(catalog_data)
    println("  Image dimensions $(geometry.height_px) H x $(geometry.width_px) W px")
    wcs = WCS.WCSTransform(
        2, # dimensions
        # reference pixel coordinates...
        crpix=[1., 1.],
        # ...and corresponding reference world coordinates
        crval=[geometry.world_coordinate_origin[1], geometry.world_coordinate_origin[2]],
        # this WCS is a simple linear transformation
        ctype=["RA---TAN", "DEC--TAN"],
        cunit=["deg", "deg"],
        # these are [du/dx du/dy; dv/dx dv/dy]. (u, v) = world coords, (x, y) = pixel coords.
        pc=[0. geometry.ra_degrees_per_pixel; geometry.dec_degrees_per_pixel 0.],
    )
    map(1:5) do band
        make_image(
            zeros(Float32, (geometry.height_px, geometry.width_px)),
            band,
            wcs,
            make_simple_psf(psf_sigma_px),
            band_sky_level_nmgy[band],
            band_nelec_per_nmgy[band],
        )
    end
end

## Parse a FITS header from a raw string representation

# It's silly that we have to do this, but the Julia FITSIO library doesn't support this (it can only
# read a header from a file), while the Julia WCS library will only generate a header in raw string
# form. (The Julia FITSIO library also won't write a raw string header to a file, so scrap that idea
# :).) Fortunately, the FITS header format is a simple fixed-width textual format with only a few
# data types. See https://fits.gsfc.nasa.gov/fits_primer.html.

function parse_fits_header_value(raw_value::AbstractString)
    if raw_value == "T"
        true
    elseif raw_value == "F"
        false
    elseif raw_value[1] == '\''
        convert(String, raw_value[2:(end - 1)])
    elseif ismatch(r"^[0-9]+$", raw_value)
        parse(Int64, raw_value)
    else
        parse(Float64, raw_value)
    end
end

function parse_fits_header_from_string(header_string::AbstractString)
    header = FITSIO.FITSHeader(String[], [], String[])
    start_position = 1
    while start_position < length(header_string)
        field_string = header_string[start_position:(start_position + 79)]
        field_name = strip(field_string[1:8])
        if length(field_name) != 0 && field_string[9:10] == "= "
            field_parts = split(field_string[11:80], "/", limit=2)
            raw_value = field_parts[1]
            header[field_name] = parse_fits_header_value(strip(raw_value))
            if length(field_parts) == 2
                comment = strip(field_parts[2])
                if length(comment) > 0
                    FITSIO.set_comment!(header, field_name, convert(String, comment))
                end
            end
        end
        start_position += 80
    end
    header
end

function save_images_to_fits(filename::String, images::Vector{Model.Image})
    println("Writing images to $filename...")
    fits_file = FITSIO.FITS(filename, "w")
    for band_image in images
        header = parse_fits_header_from_string(WCS.to_header(band_image.wcs))
        serialize_psf_to_header(band_image.psf, header)
        header["CLSKY"] = mean(band_image.sky.sky_small) * mean(band_image.sky.calibration)
        FITSIO.set_comment!(header, "CLSKY", "Mean sky background per pixel, nMgy")
        header["CLIOTA"] = mean(band_image.nelec_per_nmgy)
        FITSIO.set_comment!(header, "CLIOTA", "Gain, nelec per nMgy")
        write(
            fits_file,
            band_image.pixels,
            header=header,
            name="CELESTE_FIELD_$(band_image.b)",
        )
    end
    close(fits_file)
end

################################################################################
# Score a set of predictions against a ground truth catalog
################################################################################

COLOR_COLUMNS = [:color_log_ratio_ug, :color_log_ratio_gr, :color_log_ratio_ri, :color_log_ratio_iz]

ABSOLUTE_ERROR_COLUMNS = vcat(
    [:de_vaucouleurs_mixture_weight, :minor_major_axis_ratio, :half_light_radius_px],
    COLOR_COLUMNS,
)

function degrees_to_diff(a, b)
    angle_between = abs(a - b) % 180
    min.(angle_between, 180 - angle_between)
end

"""
Given two results data frame, one containing ground truth (i.e Coadd)
and one containing predictions (i.e., either Primary of Celeste),
compute an a data frame containing each prediction's error.
(It's not an average of the errors, it's each error.)
Let's call the return type of this function an \"error data frame\".
"""
function get_error_df(truth::DataFrame, predicted::DataFrame)
    errors = DataFrame(
        is_saturated=predicted[:is_saturated],
    )

    predicted_galaxy = predicted[:is_star] .< .5
    true_galaxy = truth[:is_star] .< .5
    errors[:missed_stars] = ifelse.(!true_galaxy, predicted_galaxy, NA)
    errors[:missed_galaxies] = ifelse.(true_galaxy, !predicted_galaxy, NA)

    errors[:position] =
        (ARCSEC_PER_DEGREE / SDSS_ARCSEC_PER_PIXEL) .*
        angular_separation.(truth[:right_ascension_deg],
                            truth[:declination_deg],
                            predicted[:right_ascension_deg],
                            predicted[:declination_deg])

    # compare flux in both mags and nMgy for now
    errors[:reference_band_flux_mag] = abs(
        flux_to_mag.(truth[:reference_band_flux_nmgy], 3)
        .- flux_to_mag.(predicted[:reference_band_flux_nmgy], 3)
    )
    errors[:reference_band_flux_nmgy] = abs(
        truth[:reference_band_flux_nmgy] .- predicted[:reference_band_flux_nmgy]
    )
    errors[:angle_deg] = degrees_to_diff(truth[:angle_deg], predicted[:angle_deg])

    for column_symbol in ABSOLUTE_ERROR_COLUMNS
        errors[column_symbol] = abs(truth[column_symbol] - predicted[column_symbol])
    end
    for color_column in COLOR_COLUMNS
        # to match up with Stripe82Score, which used differences of mags
        errors[color_column] *= 2.5 / log(10)
    end

    errors
end

function is_good_row(truth_row::DataFrameRow, error_row::DataFrameRow, column_name::Symbol)
    if isna(error_row[column_name]) || isnan(error_row[column_name])
        return false
    elseif !isna(truth_row[:half_light_radius_px]) && truth_row[:half_light_radius_px] > 20
        return false
    elseif truth_row[:is_saturated] || error_row[:is_saturated]
        return false
    end

    if column_name in [:minor_major_axis_ratio, :half_light_radius_px, :angle_deg,
                       :de_vaucouleurs_mixture_weight]
        has_mixture_weight = !isna(truth_row[:de_vaucouleurs_mixture_weight])
        if has_mixture_weight && (0.05 < truth_row[:de_vaucouleurs_mixture_weight] < 0.95)
            return false
        end
    end
    if column_name == :angle_deg
        if !isna(truth_row[:minor_major_axis_ratio]) && truth_row[:minor_major_axis_ratio] > .6
            return false
        end
    end
    return true
end

function filter_rows(truth::DataFrame, errors::DataFrame, column_name::Symbol)
    map(zip(eachrow(truth), eachrow(errors))) do rows
        (truth_row, error_row) = rows
        is_good_row(truth_row, error_row, column_name)
    end
end

function score_column(errors::DataArray)
    DataFrame(
        N=length(errors),
        first=mean(errors),
    )
end

function score_column(first_errors::DataArray, second_errors::DataArray)
    @assert length(first_errors) == length(second_errors)
    scores = score_column(first_errors)
    scores[:second] = mean(second_errors)
    diffs = first_errors .- second_errors
    scores[:diff] = mean(diffs)
    scores[:diff_sd] = std(abs.(diffs)) / sqrt(length(diffs))
    scores
end

function get_scores_df(
    truth::DataFrame, first_errors::DataFrame, second_errors::Nullable{DataFrame}
)
    score_rows = DataFrame[]
    for column_name in names(first_errors)
        good_row = filter_rows(truth, first_errors, column_name)
        if !isnull(second_errors)
            good_row .&= filter_rows(truth, get(second_errors), column_name)
        end
        if sum(good_row) <= 1
            continue
        else
            if isnull(second_errors)
                row = score_column(first_errors[good_row, column_name])
            else
                row = score_column(
                    first_errors[good_row, column_name],
                    get(second_errors)[good_row, column_name],
                )
            end

            row[:field] = column_name
            push!(score_rows, row)
        end
    end
    vcat(score_rows...)
end

"""
    match_catalogs(truth, predictions)

Return subset of rows in `truth` and `predictions`, such that each entry in
`truth` has a corresponding match in all `predictions`. A "match" is the
closest object within `tol` degrees. The output data frames all have the same
number of rows, and rows correspond to matching objects.
"""
function match_catalogs(truth::DataFrame, prediction_dfs::Vector{DataFrame};
                        tol::Float64=(0.396 / 3600.))

    matched = trues(nrow(truth))
    idxs = Vector{Int}[]
    for prediction in prediction_dfs
        idx, dists = match_coordinates(truth[:right_ascension_deg],
                                       truth[:declination_deg],
                                       prediction[:right_ascension_deg],
                                       prediction[:declination_deg])
        matched .&= dists .< tol
        push!(idxs, idx)
    end

    # For each row in `truth`, `matched` is now true if and only if all
    # predictions have a close match.
    #
    # `idxs[i]` is same length as `truth` and contains the index in
    # `prediction_dfs[i]` that is the closest match to each row in
    # `truth`.
    matched_truth = truth[matched, :]
    matched_predictions = [prediction[idx[matched], :]
                           for (prediction, idx) in zip(prediction_dfs, idxs)]

    return matched_truth, matched_predictions
end

function score_predictions(truth::DataFrame, prediction_dfs::Vector{DataFrame})
    matched_truth, matched_prediction_dfs = match_catalogs(truth, prediction_dfs)
    error_dfs = [get_error_df(matched_truth, predictions) for predictions in matched_prediction_dfs]
    @assert length(prediction_dfs) <= 2
    get_scores_df(
        matched_truth,
        error_dfs[1],
        length(error_dfs) > 1 ? Nullable(error_dfs[2]) : Nullable{DataFrame}(),
    )
end

function get_uncertainty_df(truth::DataFrame, predictions::DataFrame)
    assert_columns_are_present(predictions, STDERR_COLUMNS)
    matched_truth, matched_prediction_dfs = match_catalogs(truth, [predictions])
    matched_predictions = matched_prediction_dfs[1]

    valid_rows = (matched_truth[:reference_band_flux_nmgy] .> 0)
    matched_truth = matched_truth[valid_rows, :]
    matched_predictions = matched_predictions[valid_rows, :]

    get_errors(column, map_fn) =
        map_fn.(matched_predictions[column]) .- map_fn.(matched_truth[column])
    get_errors(column) = get_errors(column, x -> x)
    errors = [
        get_errors(:reference_band_flux_nmgy, log),
        get_errors(:color_log_ratio_ug),
        get_errors(:color_log_ratio_gr),
        get_errors(:color_log_ratio_ri),
        get_errors(:color_log_ratio_iz),
    ]
    std_errs = [
        matched_predictions[:log_reference_band_flux_stderr],
        matched_predictions[:color_log_ratio_ug_stderr],
        matched_predictions[:color_log_ratio_gr_stderr],
        matched_predictions[:color_log_ratio_ri_stderr],
        matched_predictions[:color_log_ratio_iz_stderr],
    ]
    names = [:log_reference_band_flux_nmgy, :color_log_ratio_ug, :color_log_ratio_gr,
             :color_log_ratio_ri, :color_log_ratio_iz]

    mapreduce(vcat, zip(names, errors, std_errs)) do values
        name, error, std_err = values
        DataFrame(name=fill(name, length(error)),
                  error=error,
                  posterior_std_err=std_err)
    end
end

function score_uncertainty(uncertainty_df::DataFrame)
    mapreduce(vcat, groupby(uncertainty_df, :name)) do group_df
        abs_error_sds = abs.(group_df[:error] ./ group_df[:posterior_std_err])
        abs_error_sds = abs_error_sds[!isna(abs_error_sds)]
        DataFrame(
            field=group_df[1, :name],
            within_half_sd=mean(abs_error_sds .<= 1/2),
            within_1_sd=mean(abs_error_sds .<= 1),
            within_2_sd=mean(abs_error_sds .<= 2),
            within_3_sd=mean(abs_error_sds .<= 3),
        )
    end
end

################################################################################
# Utilities
################################################################################

# Run Celeste with any combination of single/joint inference
function run_celeste(
    config::Config, catalog_entries, target_sources, images;
    use_joint_inference=false,
)
    neighbor_map = Infer.find_neighbors(target_sources, catalog_entries, images)
    if use_joint_inference
        ParallelRun.one_node_joint_infer(
            config,
            catalog_entries,
            target_sources,
            neighbor_map,
            images,
        )
    else
        ParallelRun.one_node_single_infer(
            config,
            catalog_entries,
            target_sources,
            neighbor_map,
            images,
        )
    end
end

end # module AccuracyBenchmark


text:
# exported exponential functions

"""
    ldexp(a, n)

Computes `a × 2^n`
"""
@inline ldexp(x::Union{Float32,Float64}, q::Int) = ldexpk(x, q)

const max_exp2(::Type{Float64}) = 1024
const max_exp2(::Type{Float32}) = 128f0

const min_exp2(::Type{Float64}) = -1075
const min_exp2(::Type{Float32}) = -150f0

@inline function exp2_kernel(x::Float64)
    c11 = 0.4434359082926529454e-9
    c10 = 0.7073164598085707425e-8
    c9  = 0.1017819260921760451e-6
    c8  = 0.1321543872511327615e-5
    c7  = 0.1525273353517584730e-4
    c6  = 0.1540353045101147808e-3
    c5  = 0.1333355814670499073e-2
    c4  = 0.9618129107597600536e-2
    c3  = 0.5550410866482046596e-1
    c2  = 0.2402265069591012214
    c1  = 0.6931471805599452862
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11
end

@inline function exp2_kernel(x::Float32)
    c6 = 0.1535920892f-3
    c5 = 0.1339262701f-2
    c4 = 0.9618384764f-2
    c3 = 0.5550347269f-1
    c2 = 0.2402264476f0
    c1 = 0.6931471825f0
    return @horner x c1 c2 c3 c4 c5 c6
end

"""
    exp2(x)

Compute the base-`2` exponential of `x`, that is `2ˣ`.
"""
@inline function exp2(d::T) where {T<:Union{Float32,Float64}}
    q = round(d)
    qi = unsafe_trunc(Int, q)

    s = d - q

    u = exp2_kernel(s)
    u = T(dnormalize(dadd(T(1.0), dmul(u,s))))

    u = ldexp2k(u, qi)

    d > max_exp2(T) && (u = T(Inf))
    d < min_exp2(T) && (u = T(0.0))
    return u
end

const max_exp10(::Type{Float64}) = 3.08254715559916743851e2 # log 2^1023*(2-2^-52)
const max_exp10(::Type{Float32}) = 38.531839419103626f0 # log 2^127 *(2-2^-23)

const min_exp10(::Type{Float64}) = -3.23607245338779784854769e2 # log10 2^-1075
const min_exp10(::Type{Float32}) = -45.15449934959718f0         # log10 2^-150

@inline function exp10_kernel(x::Float64)
    c11 = 0.2411463498334267652e-3
    c10 = 0.1157488415217187375e-2
    c9  = 0.5013975546789733659e-2
    c8  = 0.1959762320720533080e-1
    c7  = 0.6808936399446784138e-1
    c6  = 0.2069958494722676234e0
    c5  = 0.5393829292058536229e0
    c4  = 0.1171255148908541655e1
    c3  = 0.2034678592293432953e1
    c2  = 0.2650949055239205876e1
    c1  = 0.2302585092994045901e1
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11
end

@inline function exp10_kernel(x::Float32)
    c6 = 0.2064004987f0
    c5 = 0.5417877436f0
    c4 = 0.1171286821f1
    c3 = 0.2034656048f1
    c2 = 0.2650948763f1
    c1 = 0.2302585125f1
    return @horner x c1 c2 c3 c4 c5 c6
end

"""
    exp10(x)

Compute the base-`10` exponential of `x`, that is `10ˣ`.
"""
@inline function exp10(d::T) where {T<:Union{Float32,Float64}}
    q = round(T(MLOG10_2) * d)
    qi = unsafe_trunc(Int, q)

    s = muladd(q, -L10U(T), d)
    s = muladd(q, -L10L(T), s)

    u = exp10_kernel(s)
    u = T(dnormalize(dadd(T(1.0), dmul(u,s))))

    u = ldexp2k(u, qi)

    d > max_exp10(T) && (u = T(Inf))
    d < min_exp10(T) && (u = T(0.0))

    return u
end

const max_expm1(::Type{Float64}) = 7.09782712893383996732e2 # log 2^1023*(2-2^-52)
const max_expm1(::Type{Float32}) = 88.72283905206835f0 # log 2^127 *(2-2^-23)

const min_expm1(::Type{Float64}) = -37.42994775023704434602223
const min_expm1(::Type{Float32}) = -17.3286790847778338076068394f0

"""
    expm1(x)

Compute `eˣ- 1` accurately for small values of `x`.
"""
@inline function expm1(x::T) where {T<:Union{Float32,Float64}}
    u = T(dadd2(expk2(Double(x)), -T(1.0)))
    x > max_expm1(T) && (u = T(Inf))
    x < min_expm1(T) && (u = -T(1.0))
    isnegzero(x) && (u = T(-0.0))
    return u
end

const max_exp(::Type{Float64}) = 709.78271114955742909217217426  # log 2^1023*(2-2^-52)
const max_exp(::Type{Float32}) = 88.72283905206835f0             # log 2^127 *(2-2^-23)

const min_exp(::Type{Float64}) = -7.451332191019412076235e2 # log 2^-1075
const min_exp(::Type{Float32}) = -103.97208f0               # ≈ log 2^-150

@inline function exp_kernel(x::Float64)
    c11 = 2.08860621107283687536341e-09
    c10 = 2.51112930892876518610661e-08
    c9  = 2.75573911234900471893338e-07
    c8  = 2.75572362911928827629423e-06
    c7  = 2.4801587159235472998791e-05
    c6  = 0.000198412698960509205564975
    c5  = 0.00138888888889774492207962
    c4  = 0.00833333333331652721664984
    c3  = 0.0416666666666665047591422
    c2  = 0.166666666666666851703837
    c1  = 0.50
    return @horner x c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11
end

@inline function exp_kernel(x::Float32)
    c6 = 0.000198527617612853646278381f0
    c5 = 0.00139304355252534151077271f0
    c4 = 0.00833336077630519866943359f0
    c3 = 0.0416664853692054748535156f0
    c2 = 0.166666671633720397949219f0
    c1 = 0.5f0
    return @horner x c1 c2 c3 c4 c5 c6
end

"""
    exp(x)

Compute the base-`e` exponential of `x`, that is `eˣ`.
"""
@inline function exp(d::T) where {T<:Union{Float32,Float64}}
    q = round(T(MLN2E) * d)
    qi = unsafe_trunc(Int, q)

    s = muladd(q, -L2U(T), d)
    s = muladd(q, -L2L(T), s)

    u = exp_kernel(s)

    u = s * s * u + s + 1
    u = ldexp2k(u, qi)

    d > max_exp(T) && (u = T(Inf))
    d < min_exp(T) && (u = T(0))

    return u
end


text:
module DiffEqProblemLibrary

module ODEProblemLibrary
function importodeproblems()
  @isdefined(prob_ode_linear) ||
  include(joinpath(@__DIR__, "ode/ode_premade_problems.jl"));
  nothing
end
end # module

module DAEProblemLibrary
function importdaeproblems()
  @isdefined(prob_dae_resrob) ||
  include(joinpath(@__DIR__, "dae_premade_problems.jl"));
  nothing
end
end # module

module DDEProblemLibrary
function importddeproblems()
  @isdefined(prob_dde_1delay) ||
  include(joinpath(@__DIR__, "dde_premade_problems.jl"));
  nothing
end
end # module

module SDEProblemLibrary
function importsdeproblems()
  @isdefined(prob_sde_wave) ||
  include(joinpath(@__DIR__, "sde_premade_problems.jl"));
  nothing
end
end # module

module JumpProblemLibrary
function importjumpproblems()
  @isdefined(prob_jump_dnarepressor) ||
  include(joinpath(@__DIR__, "jump_premade_problems.jl"));
  nothing
end
end # module

end # module


text:
function getoperator(op, t)
    if op isa Tuple
        op = Semiring(op...)
    end
    #Default Semiring should be LOR_LAND for boolean
    if op == Semirings.PLUS_TIMES
        if t == Bool
            op = Semirings.LOR_LAND
        end
    end
    #Default BinaryOp should be LAND or LOR for boolean
    if op == BinaryOps.TIMES
        if t == Bool
            op = BinaryOps.LAND
        end
    elseif op == BinaryOps.PLUS
        if t == Bool
            op = BinaryOps.LOR
        end
    end
    #Default monoid should be LAND/LOR
    if op == Monoids.TIMES_MONOID
        if t == Bool
            op = Monoids.LAND_MONOID
        end
    elseif op == Monoids.PLUS_MONOID
        if t == Bool
            op = Monoids.LOR_MONOID
        end
    end

    if op isa AbstractOp
        return op[t]
    else
        return op
    end
end

getaccum(::Nothing, t) = C_NULL
getaccum(op::Function, t) = getoperator(BinaryOp(op), t)
getaccum(op, t) = getoperator(op, t)

_isloaded(op::AbstractOp) = !isempty(getfield(op, :typedops))
_isloaded(op::Union{AbstractSelectOp, AbstractDescriptor}) = getfield(op, :p) != C_NULL
"""
    validtypes(operator::AbstractOp)::Vector{DataType}
    validtypes(operator::SelectOp)::Nothing

Determine the types available as a domain for a particular operator.
Each operator is defined on a specific set of types, for instance the [`Semirings.LAND_LOR`](@ref)
    semiring is only defined for `Boolean` arguments.

When applied to an `AbstractSelectOp` this will return `nothing`.
When applied to certain operators like positional semirings it will return `[Any]`.
"""
function validtypes(o::AbstractOp)
    if !_isloaded(o)
        _load(o)
    end
    return collect(keys(o.typedops))
end

function Base.getindex(o::AbstractOp, t::DataType)
    _isloaded(o) || _load(o)
    if Any ∈ keys(o.typedops)
        getindex(o.typedops, Any)
    else
        if !haskey(o.typedops, t)
            addtoop(o, t)
        end
        getindex(o.typedops, t)
    end
end

function Base.getindex(o::AbstractBinaryOp, t1::DataType, t2::DataType)
    _isloaded(o) || _load(o)
    if (Any, Any) ∈ keys(o.typedops)
        getindex(o.typedops, (Any, Any))
    else
        if !haskey(o.typedops, (t1, t2))
            addtoop(o, (t1, t2))
        end
        getindex(o.typedops, (t1, t2))
    end
end
Base.getindex(o::AbstractBinaryOp, tup::Tuple{DataType, DataType}) = o[tup...]
Base.getindex(o::AbstractBinaryOp, t::DataType) = o[t, t]

Base.setindex!(o::AbstractBinaryOp, x, tup::Tuple{DataType, DataType}) = setindex!(o.typedops, x, tup)
Base.setindex!(o::AbstractBinaryOp, x, t1::DataType, t2::DataType) = setindex!(o, x, (t1, t2))
Base.setindex!(o::AbstractBinaryOp, x, t::DataType) = setindex!(o, x, (t, t))

function addtoop(op::AbstractUnaryOp, type)
    f = Base.invokelatest(juliaop, op)
    resulttypes = Base.return_types(f, (type,))
    if length(resulttypes) != 1
        throw(ArgumentError("Inferred more than one result type for function $(string(f)) on type $type."))
    end
    UnaryOps._addunaryop(op, f, resulttypes[1], type)
end

function addtoop(op::AbstractBinaryOp, type1, type2)
    f = Base.invokelatest(juliaop, op)
    resulttypes = Base.return_types(f, (type1, type2))
    if length(resulttypes) != 1
        throw(ArgumentError("Inferred more than one result type for function $(string(f)) on type $type."))
    end
    BinaryOps._addbinaryop(op, f, resulttypes[1], type1, type2)
end
addtoop(op::AbstractBinaryOp, type) = addtoop(op, type, type)

function Base.show(io::IO, ::MIME"text/plain", o::AbstractOp)
    print(io, o.name, ": ", validtypes(o))
end

juliaop(op) = nothing
juliaop(op::AbstractMonoid) = juliaop(op(op))


text:
function default_theme(scene::SceneLike, ::Type{<:Plot(Polytope)})
  Theme(
    color=:black,
    strokewidth=2
  )
end

import AbstractPlotting: Attributes

"""
  Function to plot general Polytope objects.
"""
function AbstractPlotting.plot!(plot::Plot(Polytope))
  pt = to_value(plot[1])
  v = coordinates.(vertices(pt))
  scatter!(plot, v, color=plot[:color])
  edges = facets(pt)
  for edge in edges
      v = edge.vertices
      vc = coordinates.(v)
      lines!(plot, vc)
  end
  plot
end


text:
"""
# module DataStructure

# Examples

```jldoctest
julia>
```
"""
module DataStructure

export Cell

struct Cell
    lattice::AbstractMatrix
    positions::AbstractMatrix
    numbers::AbstractVector
end

end


text:
"""
Perdew, Wang correlation from 1992 (10.1103/PhysRevB.45.13244)
"""
function energy_per_particle(::Val{:lda_c_pw}, ρ; improved_constants=false)
    α₁ = (0.21370,  0.20548,  0.11125)
    β₁ = (7.5957,  14.1189,  10.357)
    β₂ = (3.5876,   6.1977,   3.6231)
    β₃ = (1.6382,   3.3662,   0.88026)
    β₄ = (0.49294,  0.62517,  0.49671)
    # constant p = 1 is hard-coded in the expression for G below

    if !improved_constants  # Constants as given in original publication
        A   = (0.031091, 0.015545, 0.016887)
        f′′ = 1.709921  # f′′(0)
    else  # Modified constants, computed at improved accuracy
        A   = (0.0310907, 0.01554535, 0.0168869)
        f′′ = 8 / (9 * 2(cbrt(2) - 1))  # f′′(0)
    end

    energy_per_particle_c_pw(ρ; A, α₁, β₁, β₂, β₃, β₄, f′′)
end

function energy_per_particle_c_pw(ρ; A, α₁, β₁, β₂, β₃, β₄, f′′)
    function G(sqrt_rₛ, A, α₁, β₁, β₂, β₃, β₄)  # (10) with p = 1 hard-coded
        denom = β₁ * sqrt_rₛ + β₂ * sqrt_rₛ^2 + β₃ * sqrt_rₛ^3 + β₄ * sqrt_rₛ^4
        -2A * (1 + α₁*sqrt_rₛ^2) * log(1 + 1 / (2A * denom) )
    end

    # equation (9)
    f(ζ) = ((1+ζ)^(4/3) + (1-ζ)^(4/3) - 2)/(2^(4/3) - 2)  # == 0 for non-spin-polarised

    ε_0(rₛ) =  G(sqrt(rₛ), A[1], α₁[1], β₁[1], β₂[1], β₃[1], β₄[1])  # ε_c(rₛ, 0)
    ε_1(rₛ) =  G(sqrt(rₛ), A[2], α₁[2], β₁[2], β₂[2], β₃[2], β₄[2])  # ε_c(rₛ, 1)
    α(rₛ)   = -G(sqrt(rₛ), A[3], α₁[3], β₁[3], β₂[3], β₃[3], β₄[3])  # α_c(rₛ)

    # equation (8)
    ε(rₛ, ζ) = ε_0(rₛ) + α(rₛ) * f(ζ)/f′′ * (1 - ζ^4) + (ε_1(rₛ) - ε_0(rₛ)) * f(ζ) * ζ^4

    rₛ = cbrt(3 / (4π  * ρ))  # equation (1)
    ε(rₛ, #= ζ = =# 0)
end


text:
module Test

using Monads
using Base.Test
using Promisables

@testset "Tests for list monad (MList)" begin
    @test mreturn(MList, 1) == MList([1])
    @test mbind(x -> MList([x, 3]), MList([1,2])) == MList([1,3,2,3])
    @test (@mdo MList begin
                a <- MList([1,2,3])
                return (a + 1)
              end) == MList([2,3,4])
    @test (@mdo MList begin
                a <- MList([1,2,3])
                b <- MList([1,2,3])
                guard(a!=b)
                return (a,b)
              end) == MList([(1, 2), (1, 3), (2, 1), (2, 3), (3, 1), (3, 2)])

end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the `EnhancedTestSet` type.

This set of unit tests checks the output of `EnhancedTestSet` for failing tests.

Notes
-----
* For the unit tests in this files, failures and errors are expected.
"""

# --- Imports

# Standard library
using Test

# External packages
using Suppressor

# Local modules
using TestTools.jltest

# --- Tests

# ------ Failing tests with diffs

# Array equality test
output = strip(@capture_out begin
    @testset EnhancedTestSet "EnhancedTestSet: Array equality test" begin
        @test [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]
    end
end)

expected_prefix = strip("""
                  =====================================================
                  EnhancedTestSet: Array equality test: Test Failed
                    Expression: [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]

                    Diff:
                  [3, 5, 6, 1, (-)6, (+)9, 8]
                  """)
@test startswith(output, expected_prefix)

# Dict equality test
output = strip(
    @capture_out begin
        @testset EnhancedTestSet "EnhancedTestSet: Dict equality test" begin
            @test Dict(:foo => "bar", :baz => [1, 4, 5], :biz => nothing) ==
                Dict(:baz => [1, 7, 5], :biz => 42)
        end
    end
)

expected_prefix = strip(
    """
=====================================================
EnhancedTestSet: Dict equality test: Test Failed
  Expression: Dict(:foo => \"bar\", :baz => [1, 4, 5], :biz => nothing) == Dict(:baz => [1, 7, 5], :biz => 42)

  Diff:
[Dict{Symbol, Any}, (-):biz => nothing, (-):baz => [1, 4, 5], (-):foo => \"bar\", (+):biz => 42, (+):baz => [1, 7, 5]]

""",
)

@test startswith(output, expected_prefix)

# String equality test
output = strip(
    @capture_out begin
        @testset EnhancedTestSet "EnhancedTestSet: String equality test" begin
            @test """Lorem ipsum dolor sit amet,
                     consectetur adipiscing elit, sed do
                     eiusmod tempor incididunt ut
                     labore et dolore magna aliqua.
                     Ut enim ad minim veniam, quis nostrud
                     exercitation ullamco aboris.""" == """Lorem ipsum dolor sit amet,
                                                           consectetur adipiscing elit, sed do
                                                           eiusmod temper incididunt ut
                                                           labore et dolore magna aliqua.
                                                           Ut enim ad minim veniam, quis nostrud
                                                           exercitation ullamco aboris."""
        end
    end
)

expected_prefix = join(
    [
        "=====================================================",
        "EnhancedTestSet: String equality test: Test Failed",
        "  Expression: \"Lorem ipsum dolor sit amet,\\nconsectetur adipiscing elit, " *
        "sed do\\neiusmod tempor incididunt ut\\nlabore et dolore magna aliqua.\\nUt " *
        "enim ad minim veniam, quis nostrud\\nexercitation ullamco aboris.\" " *
        "== \"Lorem ipsum dolor sit amet,\\nconsectetur adipiscing elit, sed do\\n" *
        "eiusmod temper incididunt ut\\nlabore et dolore magna aliqua.\\nUt enim ad " *
        "minim veniam, quis nostrud\\nexercitation ullamco aboris.\"",
        "",
        "  Diff:",
        "\"\"\"",
        "  Lorem ipsum dolor sit amet,",
        "  consectetur adipiscing elit, sed do",
        "- eiusmod tempor incididunt ut",
        "+ eiusmod temper incididunt ut",
        "  labore et dolore magna aliqua.",
        "  Ut enim ad minim veniam, quis nostrud",
        "  exercitation ullamco aboris.\"\"\"",
        "",
    ],
    "\n",
)

@test startswith(output, expected_prefix)

# ------ Failing tests without diffs

# Boolean expression test
output = strip(@capture_out begin
    @testset EnhancedTestSet "EnhancedTestSet: Boolean expression test" begin
        @test iseven(7)
    end
end)

expected_prefix = strip("""
                  =====================================================
                  EnhancedTestSet: Boolean expression test: Test Failed at $(@__FILE__):132
                    Expression: iseven(7)

                  Stacktrace:

                  """)

@test startswith(output, expected_prefix)

# Exception test
output = strip(@capture_out begin
    @testset EnhancedTestSet "EnhancedTestSet: Exception test" begin
        throw(ErrorException("This test is supposed to throw an error"))
    end
end)

expected_prefix = strip("""
                  =====================================================
                  EnhancedTestSet: Exception test: Error During Test at $(@__FILE__):149
                    Got exception outside of a @test
                    This test is supposed to throw an error
                    Stacktrace:
                  """)

@test startswith(output, expected_prefix)

# Inequality test
output = strip(@capture_out begin
    @testset EnhancedTestSet "EnhancedTestSet: inequality test" begin
        @test 1 > 2
    end
end)

expected_prefix = strip("""
                  =====================================================
                  EnhancedTestSet: inequality test: Test Failed at $(@__FILE__):167
                    Expression: 1 > 2
                     Evaluated: 1 > 2

                  Stacktrace:
                  """)

@test startswith(output, expected_prefix)

# Matrix equality test
output = strip(@capture_out begin
    @testset EnhancedTestSet "EnhancedTestSet: Matrix equality test" begin
        @test [1 2; 3 4] == [1 4; 3 4]
    end
end)

expected_prefix = strip("""
                  =====================================================
                  EnhancedTestSet: Matrix equality test: Test Failed
                    Expression: [1 2; 3 4] == [1 4; 3 4]

                    Diff:
                  nothing
                  """)
@test startswith(output, expected_prefix)

# --- Emit message about expected failures and errors

println()
@info "For $(basename(@__FILE__)), 6 failures and 1 error are expected."


text:
module SSLTestGenerated
using SimpleSDMLayers
using Test

M = rand(Float64, (5, 10))
S = SimpleSDMPredictor(M, 0.0, 1.0, 0.0, 1.0)

@test maximum(S) ≤ 1.0
@test minimum(S) ≥ 0.0
@test extrema(S) == (minimum(S), maximum(S))

using Statistics

@test mean(S) ≈ mean(M)
@test std(S) ≈ std(M)
@test median(S) ≈ median(M)

M = ones(Float64, (5, 10))
S = SimpleSDMPredictor(M, 0.0, 1.0, 0.0, 1.0)
@test typeof(sqrt(S)) <: SimpleSDMLayer

S1 = SimpleSDMResponse([-1 0 1 nothing], 0.0, 1.0, 0.0, 1.0)
S2 = SimpleSDMPredictor([Float32.([-2.0 1.0 1.0])... nothing], 0.0, 1.0, 0.0, 1.0)
Sx = [S1, S2]
@test reduce(min, Sx).grid == [-2 0 1 nothing]
@test mean(Sx).grid == [Float32.([-1.5 0.5 1.0])... nothing]
@test reduce(+, Sx).grid == [-3 1 2 nothing]

end


text:
#
# Autogenerated by Thrift Compiler (1.0.0-dev)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING


text:
using SPIRV_LLVM_Translator_jll, SPIRV_Tools_jll

@testset "SPIR-V" begin

include("definitions/spirv.jl")

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

@testset "IR" begin

@testset "kernel functions" begin
@testset "calling convention" begin
    kernel() = return

    ir = sprint(io->spirv_code_llvm(io, kernel, Tuple{}; dump_module=true))
    @test !occursin("spir_kernel", ir)

    ir = sprint(io->spirv_code_llvm(io, kernel, Tuple{};
                                    dump_module=true, kernel=true))
    @test occursin("spir_kernel", ir)
end

@testset "byval workaround" begin
    kernel(x) = return

    ir = sprint(io->spirv_code_llvm(io, kernel, Tuple{Tuple{Int}}))
    @test occursin(r"@.*julia_.+_kernel.+\(({ i64 }|\[1 x i64\])\*", ir)

    ir = sprint(io->spirv_code_llvm(io, kernel, Tuple{Tuple{Int}}; kernel=true))
    @test occursin(r"@.*julia_.+_kernel.+\({ ({ i64 }|\[1 x i64\]) }\*.+byval", ir)
end

@testset "byval bug" begin
    # byval added alwaysinline, which could conflict with noinline and fail verification
    @noinline kernel() = return
    spirv_code_llvm(devnull, kernel, Tuple{}; kernel=true)
end
end

end

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

@testset "asm" begin

@testset "trap removal" begin
    function kernel(x)
        x && error()
        return
    end

    asm = sprint(io->spirv_code_native(io, kernel, Tuple{Bool}; kernel=true))
    @test occursin(r"OpFunctionCall %void %julia_error", asm)
end

end

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

end


text:
module Space

using ..Util
import ..Util: increment!, reset!
using Spec
import Flux     # FIXME: Can we excise this from this submodule?
import ForwardDiff
import Random
using Random: GLOBAL_RNG, AbstractRNG
using DocStringExtensions
using DataStructures: LinkedList, cons, nil, list, head, tail

export Ω, uid, @id, ID
export ΩBase, memrand, linearize, unlinearize
export append, base, combine, increment!, increment,
       Paired, pair
export ΩProj, parentω, memrand
export TaggedΩ, tag, Tags, hastags
export SimpleΩ, LinearΩ, update, nelem

include("idgen.jl")         # Id generation
include("index.jl")         # Pairing functions for omega ids
include("omega.jl")         # Sample Space
include("proj.jl")          # Sample Space Projection
include("tagged.jl")        # Space space Tagged with metadata
include("common.jl")
include("seed.jl")          # Setting the seed

# Omega Types
include("simple.jl")        # Simple Ω
include("linear.jl")        # Linear Ω
include("rng.jl")           # Linear Ω

end


text:
# Geometry which displays points at given (x, y) positions.
immutable PointGeometry <: Gadfly.GeometryElement
end

const point = PointGeometry

function element_aesthetics(::PointGeometry)
    [:x, :y, :size, :color, :opacity]
end

# Generate a form for a point geometry.
#
# Args:
#   geom: point geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::PointGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.point", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.point", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    default_aes.opacity = Float64[theme.default_point_opacity]
    default_aes.size = Measure[theme.default_point_size]
    aes = inherit(aes, default_aes)

    lw_hover_scale = 10
    lw_ratio = theme.line_width / aes.size[1]
    compose(circle(aes.x, aes.y, aes.size),
            fill(aes.color),
            opacity(aes.opacity),
            stroke([theme.highlight_color(c) for c in aes.color]),
            linewidth(theme.line_width),
            d3embed(@sprintf(".on(\"mouseover\", geom_point_mouseover(%0.2f, %0.2f), false)",
                             lw_hover_scale, lw_ratio)),
            d3embed(@sprintf(".on(\"mouseout\", geom_point_mouseout(%0.2f, %0.2f), false)",
                             lw_hover_scale, lw_ratio)),
            aes.color_key_continuous == true ?
                svgclass("geometry") :
                svgclass([@sprintf("geometry color_%s", escape_id(aes.color_label([c])[1]))
                          for c in aes.color]))
end