text:
abstract type AbstractHydroFormulation <: AbstractDeviceFormulation end
abstract type AbstractHydroDispatchFormulation <: AbstractHydroFormulation end
abstract type AbstractHydroUnitCommitment <: AbstractHydroFormulation end
abstract type AbstractHydroReservoirFormulation <: AbstractHydroDispatchFormulation end
struct HydroDispatchRunOfRiver <: AbstractHydroDispatchFormulation end
struct HydroDispatchReservoirBudget <: AbstractHydroReservoirFormulation end
struct HydroDispatchReservoirStorage <: AbstractHydroReservoirFormulation end
struct HydroDispatchPumpedStorage <: AbstractHydroReservoirFormulation end
struct HydroDispatchPumpedStoragewReservation <: AbstractHydroReservoirFormulation end
struct HydroCommitmentRunOfRiver <: AbstractHydroUnitCommitment end
struct HydroCommitmentReservoirBudget <: AbstractHydroUnitCommitment end
struct HydroCommitmentReservoirStorage <: AbstractHydroUnitCommitment end

########################### Hydro generation variables #################################

"""
This function add the variables for active power to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ActivePowerVariable, U <: PSY.HydroGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_active,
        initial_value_func = x -> PSY.get_active_power(x),
        lb_value_func = x -> PSY.get_active_power_limits(x).min,
        ub_value_func = x -> PSY.get_active_power_limits(x).max,
    )
end

"""
This function add the variables for reactive power to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ReactivePowerVariable, U <: PSY.HydroGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_reactive,
        initial_value_func = x -> PSY.get_reactive_power(x),
        lb_value_func = x -> PSY.get_reactive_power_limits(x).min,
        ub_value_func = x -> PSY.get_reactive_power_limits(x).max,
    )
end

"""
This function add the variables for energy storage to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: EnergyVariable, U <: PSY.HydroGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        initial_value_func = x -> PSY.get_initial_storage(x),
        lb_value_func = x -> 0.0,
        ub_value_func = x -> PSY.get_storage_capacity(x),
    )
end

"""
This function add the variables for upper energy storage to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: EnergyVariableUp, U <: PSY.HydroGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        initial_value_func = x -> PSY.get_initial_storage(x).up,
        lb_value_func = x -> 0.0,
        ub_value_func = x -> PSY.get_storage_capacity(x).up,
    )
end

"""
This function add the variables for lower energy storage to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: EnergyVariableDown, U <: PSY.HydroGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        initial_value_func = x -> PSY.get_initial_storage(x).down,
        lb_value_func = x -> 0.0,
        ub_value_func = x -> PSY.get_storage_capacity(x).down,
    )
end

"""
This function add the variables for active power withdrawl to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ActivePowerInVariable, U <: PSY.HydroGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_active,
        sign = -1.0,
        lb_value_func = x -> 0.0,
    )
end

"""
This function add the variables for active power injection to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ActivePowerOutVariable, U <: PSY.HydroGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_active,
        lb_value_func = x -> 0.0,
    )
end

"""
This function add the variables for power generation commitment to the model
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    psi_container::PSIContainer,
) where {T <: OnVariable, U <: PSY.HydroGen}
    return AddVariableSpec(; variable_name = make_variable_name(T, U), binary = true)
end

"""
This function add the spillage variable for storage models
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: SpillageVariable, U <: PSY.HydroGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        lb_value_func = x -> 0.0,
    )
end

"""
This function adds the reservation variable for storage models
"""
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ReserveVariable, U <: PSY.HydroGen}
    return AddVariableSpec(; variable_name = make_variable_name(T, U), binary = true)
end

"""
This function define the range constraint specs for the
reactive power for dispatch formulations.
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ReactivePowerVariable,
                T,
            ),
            variable_name = make_variable_name(ReactivePowerVariable, T),
            limits_func = x -> PSY.get_reactive_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function define the range constraint specs for the
active power for dispatch Run of River formulations.
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    if !use_parameters && !use_forecasts
        return DeviceRangeConstraintSpec(;
            range_constraint_spec = RangeConstraintSpec(;
                constraint_name = make_constraint_name(
                    RangeConstraint,
                    ActivePowerVariable,
                    T,
                ),
                variable_name = make_variable_name(ActivePowerVariable, T),
                limits_func = x -> (min = 0.0, max = PSY.get_active_power(x)),
                constraint_func = device_range!,
                constraint_struct = DeviceRangeConstraintInfo,
            ),
        )
    end

    return DeviceRangeConstraintSpec(;
        timeseries_range_constraint_spec = TimeSeriesConstraintSpec(
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            parameter_name = use_parameters ? ACTIVE_POWER : nothing,
            forecast_label = "max_active_power",
            multiplier_func = x -> PSY.get_max_active_power(x),
            constraint_func = use_parameters ? device_timeseries_param_ub! :
                              device_timeseries_ub!,
        ),
    )
end

"""
This function define the range constraint specs for the
active power for dispatch Reservoir formulations.
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroReservoirFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function define the range constraint specs for the
active power for commitment formulations (semi continuous).
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroUnitCommitment},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            bin_variable_names = [make_variable_name(OnVariable, T)],
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function define the range constraint specs for the
reactive power for commitment formulations (semi continuous).
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroUnitCommitment},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ReactivePowerVariable,
                T,
            ),
            variable_name = make_variable_name(ReactivePowerVariable, T),
            bin_variable_names = [make_variable_name(OnVariable, T)],
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerOutVariable},
    ::Type{T},
    ::Type{<:HydroDispatchPumpedStorage},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerOutVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerOutVariable, T),
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerInVariable},
    ::Type{T},
    ::Type{<:HydroDispatchPumpedStorage},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerInVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerInVariable, T),
            limits_func = x -> PSY.get_active_power_limits_pump(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerOutVariable},
    ::Type{T},
    ::Type{<:HydroDispatchPumpedStoragewReservation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerOutVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerOutVariable, T),
            bin_variable_names = [make_variable_name(ReserveVariable, T)],
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = reserve_device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerInVariable},
    ::Type{T},
    ::Type{<:HydroDispatchPumpedStoragewReservation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerInVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerInVariable, T),
            bin_variable_names = [make_variable_name(ReserveVariable, T)],
            limits_func = x -> PSY.get_active_power_limits_pump(x),
            constraint_func = reserve_device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end
######################## RoR constraints ############################

"""
This function define the range constraint specs for the
reactive power for Commitment Run of River formulation.
    `` P <= multiplier * P_max ``
"""
function commit_hydro_active_power_ub!(
    psi_container::PSIContainer,
    devices,
    model::DeviceModel{V, W},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {V <: PSY.HydroGen, W <: AbstractHydroUnitCommitment}
    use_parameters = model_has_parameters(psi_container)
    use_forecasts = model_uses_forecasts(psi_container)
    if use_parameters || use_forecasts
        spec = DeviceRangeConstraintSpec(;
            timeseries_range_constraint_spec = TimeSeriesConstraintSpec(
                constraint_name = make_constraint_name(
                    RangeConstraint,
                    ActivePowerVariable,
                    V,
                ),
                variable_name = make_variable_name(ActivePowerVariable, V),
                parameter_name = use_parameters ? ACTIVE_POWER : nothing,
                forecast_label = "max_active_power",
                multiplier_func = x -> PSY.get_max_active_power(x),
                constraint_func = use_parameters ? device_timeseries_param_ub! :
                                  device_timeseries_ub!,
            ),
        )
        device_range_constraints!(psi_container, devices, model, feedforward, spec)
    end
end

######################## Energy balance constraints ############################

"""
This function defines the constraints for the water level (or state of charge)
for the Hydro Reservoir.
"""
function energy_balance_constraint!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, S},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {
    H <: PSY.HydroEnergyReservoir,
    S <: Union{HydroDispatchReservoirStorage, HydroCommitmentReservoirStorage},
}
    key = ICKey(EnergyLevel, H)
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)

    if !has_initial_conditions(psi_container.initial_conditions, key)
        throw(IS.DataFormatError("Initial Conditions for $(H) Energy Constraints not in the model"))
    end

    inflow_forecast_label = "inflow"
    target_forecast_label = "storage_target"
    constraint_infos_inflow = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    constraint_infos_target = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector_inflow = get_time_series(psi_container, d, inflow_forecast_label)
        constraint_info_inflow = DeviceTimeSeriesConstraintInfo(
            d,
            x -> PSY.get_inflow(x) * PSY.get_conversion_factor(x),
            ts_vector_inflow,
        )
        add_device_services!(constraint_info_inflow.range, d, model)
        constraint_infos_inflow[ix] = constraint_info_inflow

        ts_vector_target = get_time_series(psi_container, d, target_forecast_label)
        constraint_info_target = DeviceTimeSeriesConstraintInfo(
            d,
            x -> PSY.get_storage_target(x) * PSY.get_storage_capacity(x),
            ts_vector_target,
        )
        constraint_infos_target[ix] = constraint_info_target
    end

    if parameters
        energy_balance_hydro_param!(
            psi_container,
            get_initial_conditions(psi_container, key),
            (constraint_infos_inflow, constraint_infos_target),
            (
                make_constraint_name(ENERGY_CAPACITY, H),
                make_constraint_name(ENERGY_TARGET, H),
            ),
            (
                make_variable_name(SPILLAGE, H),
                make_variable_name(ACTIVE_POWER, H),
                make_variable_name(ENERGY, H),
            ),
            (
                UpdateRef{H}(INFLOW, inflow_forecast_label),
                UpdateRef{H}(TARGET, target_forecast_label),
            ),
        )
    else
        energy_balance_hydro!(
            psi_container,
            get_initial_conditions(psi_container, key),
            (constraint_infos_inflow, constraint_infos_target),
            (
                make_constraint_name(ENERGY_CAPACITY, H),
                make_constraint_name(ENERGY_TARGET, H),
            ),
            (
                make_variable_name(SPILLAGE, H),
                make_variable_name(ACTIVE_POWER, H),
                make_variable_name(ENERGY, H),
            ),
        )
    end
    return
end

"""
This function defines the constraints for the water level (or state of charge)
for the HydroPumpedStorage.
"""
function energy_balance_constraint!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, S},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {
    H <: PSY.HydroPumpedStorage,
    S <: Union{HydroDispatchPumpedStorage, HydroDispatchPumpedStoragewReservation},
}
    key = ICKey(EnergyLevelUP, H)
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)

    if !has_initial_conditions(psi_container.initial_conditions, key)
        throw(IS.DataFormatError("Initial Conditions for $(H) Energy Constraints not in the model"))
    end

    forecast_label_in = "inflow"
    constraint_infos = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label_in)
        constraint_info = DeviceTimeSeriesConstraintInfo(
            d,
            x -> PSY.get_inflow(x) * PSY.get_conversion_factor(x),
            ts_vector,
        )
        add_device_services!(constraint_info.range, d, model)
        constraint_infos[ix] = constraint_info
    end

    forecast_label_out = "outflow"
    constraint_infos_outflow =
        Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label_out)
        constraint_info = DeviceTimeSeriesConstraintInfo(
            d,
            x -> PSY.get_outflow(x) * PSY.get_conversion_factor(x),
            ts_vector,
        )
        add_device_services!(constraint_info.range, d, model)
        constraint_infos_outflow[ix] = constraint_info
    end

    if parameters
        energy_balance_hydro_param!(
            psi_container,
            get_initial_conditions(psi_container, key),
            (constraint_infos, constraint_infos_outflow),
            (
                make_constraint_name(ENERGY_CAPACITY_UP, H),
                make_constraint_name(ENERGY_CAPACITY_DOWN, H),
            ),
            (
                make_variable_name(SPILLAGE, H),
                make_variable_name(ACTIVE_POWER_OUT, H),
                make_variable_name(ENERGY_UP, H),
                make_variable_name(ACTIVE_POWER_IN, H),
                make_variable_name(ENERGY_DOWN, H),
            ),
            (
                UpdateRef{H}(INFLOW, forecast_label_in),
                UpdateRef{H}(OUTFLOW, forecast_label_out),
            ),
        )
    else
        energy_balance_hydro!(
            psi_container,
            get_initial_conditions(psi_container, key),
            (constraint_infos, constraint_infos_outflow),
            (
                make_constraint_name(ENERGY_CAPACITY_UP, H),
                make_constraint_name(ENERGY_CAPACITY_DOWN, H),
            ),
            (
                make_variable_name(SPILLAGE, H),
                make_variable_name(ACTIVE_POWER_OUT, H),
                make_variable_name(ENERGY_UP, H),
                make_variable_name(ACTIVE_POWER_IN, H),
                make_variable_name(ENERGY_DOWN, H),
            ),
        )
    end
    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    device_formulation::Type{<:AbstractHydroUnitCommitment},
) where {H <: PSY.HydroGen}
    status_init(psi_container, devices)
    output_init(psi_container, devices)
    duration_init(psi_container, devices)

    return
end

function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    device_formulation::Type{D},
) where {H <: PSY.HydroGen, D <: AbstractHydroDispatchFormulation}
    output_init(psi_container, devices)

    return
end

########################## Addition to the nodal balances #################################

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractPowerModel},
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return NodalExpressionSpec(
        "get_max_active_power",
        REACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_reactive_power(x) : x -> PSY.get_reactive_power(x),
        1.0,
        T,
    )
end

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractActivePowerModel},
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return NodalExpressionSpec(
        "max_active_power",
        ACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_active_power(x) : x -> PSY.get_active_power(x),
        1.0,
        T,
    )
end

function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.HydroEnergyReservoir},
    device_formulation::Type{D},
    system_formulation::Type{<:PM.AbstractPowerModel},
) where {D <: AbstractHydroFormulation}
    add_to_cost!(
        psi_container,
        devices,
        make_variable_name(ACTIVE_POWER, PSY.HydroEnergyReservoir),
        :fixed,
        -1.0,
    )

    return
end

function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    device_formulation::Type{D},
    system_formulation::Type{<:PM.AbstractPowerModel},
) where {D <: AbstractHydroFormulation, H <: PSY.HydroGen}
    return
end

##################################### Water/Energy Budget Constraint ############################
function energy_budget_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, <:AbstractHydroFormulation},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::IntegralLimitFF,
) where {H <: PSY.HydroGen}
    return
end

"""
This function define the budget constraint for the
active power budget formulation.

`` sum(P[t]) <= Budget ``
"""
function energy_budget_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, <:AbstractHydroFormulation},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {H <: PSY.HydroGen}
    forecast_label = "hydro_budget"
    constraint_data = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        constraint_d =
            DeviceTimeSeriesConstraintInfo(d, x -> PSY.get_storage_capacity(x), ts_vector)
        constraint_data[ix] = constraint_d
    end

    if model_has_parameters(psi_container)
        device_energy_budget_param_ub(
            psi_container,
            constraint_data,
            make_constraint_name(ENERGY_BUDGET, H),
            UpdateRef{H}(ENERGY_BUDGET, forecast_label),
            make_variable_name(ACTIVE_POWER, H),
        )
    else
        device_energy_budget_ub(
            psi_container,
            constraint_data,
            make_constraint_name(ENERGY_BUDGET),
            make_variable_name(ACTIVE_POWER, H),
        )
    end
end

"""
This function define the budget constraint (using params)
for the active power budget formulation.
"""
function device_energy_budget_param_ub(
    psi_container::PSIContainer,
    energy_budget_data::Vector{DeviceTimeSeriesConstraintInfo},
    cons_name::Symbol,
    param_reference::UpdateRef,
    var_names::Symbol,
)
    time_steps = model_time_steps(psi_container)
    resolution = model_resolution(psi_container)
    inv_dt = 1.0 / (Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR)
    variable_out = get_variable(psi_container, var_names)
    set_name = [get_component_name(r) for r in energy_budget_data]
    constraint = add_cons_container!(psi_container, cons_name, set_name)
    container = add_param_container!(psi_container, param_reference, set_name, 1)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)
    for constraint_info in energy_budget_data
        name = get_component_name(constraint_info)
        multiplier[name, 1] = constraint_info.multiplier * inv_dt
        param[name, 1] =
            PJ.add_parameter(psi_container.JuMPmodel, sum(constraint_info.timeseries))
        constraint[name] = JuMP.@constraint(
            psi_container.JuMPmodel,
            sum([variable_out[name, t] for t in time_steps]) <= multiplier[name, 1] * param[name, 1]
        )
    end

    return
end

"""
This function define the budget constraint
for the active power budget formulation.
"""
function device_energy_budget_ub(
    psi_container::PSIContainer,
    energy_budget_constraints::Vector{DeviceTimeSeriesConstraintInfo},
    cons_name::Symbol,
    var_names::Symbol,
)
    time_steps = model_time_steps(psi_container)
    variable_out = get_variable(psi_container, var_names)
    names = [get_component_name(x) for x in energy_budget_constraints]
    constraint = add_cons_container!(psi_container, cons_name, names)

    for constraint_info in energy_budget_constraints
        name = get_component_name(constraint_info)
        resolution = model_resolution(psi_container)
        inv_dt = 1.0 / (Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR)
        forecast = constraint_info.timeseries
        multiplier = constraint_info.multiplier * inv_dt
        constraint[name] = JuMP.@constraint(
            psi_container.JuMPmodel,
            sum([variable_out[name, t] for t in time_steps]) <= multiplier * sum(forecast)
        )
    end

    return
end

##################################### Hydro generation cost ############################
function AddCostSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: PSY.HydroDispatch, U <: AbstractHydroFormulation}
    # Hydro Generators currently have no OperationalCost
    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        fixed_cost = x -> 1.0,
        multiplier = OBJECTIVE_FUNCTION_NEGATIVE,
    )
end

############################
function AddCostSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: PSY.HydroGen, U <: AbstractHydroFormulation}
    # Hydro Generators currently have no OperationalCost
    cost_function = x -> isnothing(x) ? 1.0 : PSY.get_variable(x)
    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        fixed_cost = PSY.get_fixed,
        variable_cost = cost_function,
        multiplier = OBJECTIVE_FUNCTION_POSITIVE,
    )
end


text:
SNN = SpikingNeuralNetworks
param = SNN.IZParameter_more(;a = 0.02, b = 0.2, c = -65.0, d = 8.0, C=100.0, vr=-65.0, k=4.0, vt=-50.0, vPeak=20.0)

ext: jl

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
#-----------------------------------------------------------------------
#
#   OUTLIER DETECTION FUNCTIONS
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    combX( X::RealMtx )

The `combX` function constructs and returns a vector of ``4×4`` matrices
consisting of all choices of 4 column vectors from the ``4×n_{\\rm e}`` matrix
`X`.

"""
function combX( X::RealMtx )
    tpfl=typeof(X[1,1])
    np = size(X)[2]

    if np==4
        return [X]
    elseif np>4
        Z = [zeros(tpfl,4) for _ in 1:np]   # Create containers

        for i=1:np
            Z[i] = copy(X[:,i])
        end

        W = collect(combinations(Z,4))      # Combinations
        nc = length(W)
        Y = [zeros(tpfl,4,4) for _ in 1:nc]

        for i=1:nc
            Y[i][:,1] = W[i][1]
            Y[i][:,2] = W[i][2]
            Y[i][:,3] = W[i][3]
            Y[i][:,4] = W[i][4]
        end

        return Y
    else
        return [zeros(tpfl,4,4)]
    end
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    odetc( XC::Vector , ξ::Real=1e1 )

The `odetc` function takes a vector of vectors as input and returns a
vector of the vectors that differ by less than a threshold of `ξ`.

"""
function odetc( XC::Vector , ξ::Real=1e1 )
    tpfl    = typeof(XC[1][1])
    ns      = length(XC)

    Xc  = zeros(tpfl,4,ns)
    ΔX  = zeros(tpfl,4,ns)
    Xa  = zeros(tpfl,4)
    V   = zeros(tpfl,ns)
    BV  = zeros(Bool,ns)
    σ   = zeros(tpfl,ns)

    for i=1:ns
        Xc[:,i] = XC[i][:]
    end

    for i=1:4
        Xa[i] = median(Xc[i,:])
    end

    for i=1:ns
        for j=1:4
            ΔX[j,i] = (Xc[j,i] - Xa[j])
        end
    end

    for i=1:ns
        V[i] = dot( ΔX[:,i] , ΔX[:,i] )
    end

    IV1 = partialsortperm( V , 1:ns )
    nr = 0

    if ns > 1
        for i=1:ns
            σ[i]    = sqrt(sum(V[IV1][1:i])/i)
            if i==1
                BV[i]   =  true
                nr = 1
            elseif i>1
                if σ[i] < ξ
                    BV[i]   =  true
                    nr += 1
                end
            end
        end
    end
    if nr > 1
        IV2 = partialsortperm( BV , 1:ns , rev=true )
        return (XC[IV1[IV2]][1:nr],IV1[IV2])
    else
        return (XC[IV1][1:1],nr*ones(Int,ns))
    end
end     #---------------------------------------------------------------


text:
"""
    EDMD(X, basis; alg, p, t, C)
    EDMD(X, Y, basis; alg, p, t, C)

Approximates a 'NonlinearKoopman' with the `AbstractKoopmanAlgorithm` 'alg' from the data matrices `X` or `X` and `Y`, respectively.
If only `X` is given, the data is split into `X[:, 1:end-1]` and `X[:, 2:end]`.

Additional keyworded arguments include `p` for the parameter of the basis and `t` for an array of time points.
`C` is the matrix representing the mapping from the Koopman space into the state space.

# Example

```julia
@parameters p[1] t
@variables u[1:2](t)
h = Operation[u; sin.(u); cos(p[1]*t)]
basis = Basis(h, u, parameters = p, iv = t)
koopman = EDMD(X, basis, p = [2.0], t = collect(0:0.2:10.0), C = Float64[1 0 0 0 0; 0 1 0 0 0])
```
"""
function EDMD(X::AbstractArray, Ψ::AbstractBasis; p::AbstractArray = [], t::AbstractVector = [], C::AbstractArray = [], alg::AbstractKoopmanAlgorithm = DMDPINV())
    return EDMD(X[:, 1:end-1], X[:, 2:end], Ψ, p = p, t = t, C = C, alg = alg)
end

function EDMD(X::AbstractArray, Y::AbstractArray, Ψ::AbstractBasis; p::AbstractArray = [], t::AbstractVector = [], C::AbstractArray = [], alg::AbstractKoopmanAlgorithm = DMDPINV())
    @assert size(X)[2] .== size(Y)[2] "Provide consistent dimensions for data"
    @assert size(Y)[1] .<= size(Y)[2] "Provide consistent dimensions for data"

    # Based upon William et.al., A Data-Driven Approximation of the Koopman operator

    # Number of states and measurements
    N,M = size(X)

    # Compute the transformed data
    Ψ₀ = Ψ(X, p, t)
    Ψ₁ = Ψ(Y, p, t)

    A = alg(Ψ₀, Ψ₁)

    # Transform back to states
    if isempty(C)
        C = X*pinv(Ψ₀)
    end

    # TODO Maybe reduce the observable space here
    return NonlinearKoopman(A, [], C , Ψ, Ψ₁*Ψ₀', Ψ₀*Ψ₀', true)
end

function gEDMD(X::AbstractArray, DX::AbstractArray, Ψ::AbstractBasis; p::AbstractArray = [], t::AbstractVector = [], C::AbstractArray = [], alg::AbstractKoopmanAlgorithm = DMDPINV())
    @assert size(X)[2] .== size(DX)[2] "Provide consistent dimensions for data"
    @assert size(DX)[1] .<= size(DX)[2] "Provide consistent dimensions for data"

    # Based upon William et.al., A Data-Driven Approximation of the Koopman operator

    # Number of states and measurements
    N,M = size(X)

    # Compute the transformed data
    Ψ₀ = Ψ(X, p, t)

    # The jacobian to get d/dt(Ψ) = d/dx(Ψ) dx/dt
    ∇ = jacobian(Ψ)

    Ψ₁ = hcat([∇(X[:,i], p, isempty(t) ? zero(eltype(X)) : t[i])*DX[:, i] for i in 1:size(X, 2)]...)
    #return Ψ₁

    A = alg(Ψ₀, Ψ₁)

    # Transform back to states
    if isempty(C)
        C = X*pinv(Ψ₀)
    end

    # TODO Maybe reduce the observable space here
    return NonlinearKoopman(A, [], C , Ψ, Ψ₁*Ψ₀', Ψ₀*Ψ₀', false)
end

"""
    gEDMD(X, Y, basis; alg, p, t, C)
    gEDMD(t, X, basis; dt, p, C, alg, fdm, itp)

Approximates a 'NonlinearKoopman' with the `AbstractKoopmanAlgorithm` 'alg' from the data matrices `X` and `Y`.
`X` should contain the state trajectory and `Y` the differential state trajectory.

If no measurements of the differential state are available, `gEDMD` can be called with measurement time points `t` as the first argument.
It will then create an interpolation using the interpolation method from `DataInterpolations.jl` defined in `itp`. The trajectory will then be resampled
to equidistant measurements over time corresponding to the mean of `diff(t)` or `dt`, if given.
The differential state measurements will be computed via 'FiniteDifferences.jl', given a `FiniteDifferenceMethod` in `fdm`.

# Example

```julia
koopman = gEDMD(X, Y, basis)

fdm = backward_fdm(5,1)
itp = CubicSpline
koopman = gEDMD(t, X, basis, fdm = fdm, itp = itp)
```
"""
function gEDMD(t::AbstractVector, X::AbstractArray, Ψ::DataDrivenDiffEq.AbstractBasis; dt::Real = 0.0, p::AbstractArray = [], C::AbstractArray = [], alg::DataDrivenDiffEq.AbstractKoopmanAlgorithm = DMDPINV(), fdm::FiniteDifferences.FiniteDifferenceMethod = backward_fdm(5, 1), itp = CubicSpline)
    @assert size(X, 2) == length(t) "Sample size must match."
    @assert test_comp = begin

        if itp ∈ [LinearInterpolation, QuadraticInterpolation, QuadraticSpline] && !isa(fdm, FiniteDifferences.Backward{typeof(fdm.grid), typeof(fdm.coefs)})
            false
        else
            true
        end
    end "LinearInterpolation and QuadraticInterpolation need to a backward finite difference method."

    Δt = iszero(dt) ? mean(diff(t)) : dt
    t̂ = t[1]:Δt:t[end]

    X̂ = zeros(eltype(X), size(X, 1), length(t̂))
    Y = similar(X̂)
    for (i, xi) in enumerate(eachrow(X))
        itp_ = itp(xi, t)
        dx(t) = FiniteDifferences.fdm(fdm, itp_, t)
        X̂[i, :] .= itp_.(t̂)
        Y[i, :] .= dx.(t̂)
    end
    return gEDMD(X̂, Y, Ψ, alg = alg, p = p, t = collect(t̂), C = C)
end


text:
function test_KNNUCB()
    d = 2
    A = 2
    traj = Trajectory(d)
    pol = KNNUCB(A, traj, t-> 1.0, 1.0)
    ContextualBandits.add_memory!(traj, 3)

    X = [
            0.0 0.1; # a = 1
            0.0 0.1; # a = 2
            0.0 0.2 # a = 1
        ]
    for i in 1:A
        a = choose(X[i, :], pol)
        @test a == i
        ContextualBandits.extend!(traj, X[i, :], a, Float64(i), 0.0)
        ContextualBandits.update!(pol, X[i, :], a, Float64(i))
    end

    ContextualBandits.extend!(traj, X[3, :], 1, 2.0, 0.0)
    ContextualBandits.update!(pol, X[3, :], 1, 2.0)
    x0 = [0.0, 0.0]
    dst, r = ContextualBandits.sort_nearest(pol, x0, 1)
    @test dst == [0.1, 0.2]
    @test r == [1.0, 2.0]
end

test_KNNUCB()


text:
# ------------------------------------------------------------------------------------------
# ## How can we calculate $\pi$?
#
# Given a square of length $2r$, the square's area is
#
# $$A_{square} = (2r)^2 = 4r^2$$
#
# whereas the area of a circle with radius $r$ is
# $$A_{circle} = \pi r^2$$
#
#
#
# Therefore the ratio of the area of the circle to that of the square above is
#
# $$\frac{A_{circle}}{A_{square}} = \frac{\pi r^2}{4r^2} = \frac{\pi}{4}$$
#
# and we can define $\pi$ as
#
# $$\pi = 4\frac{A_{circle}}{A_{square}}$$
#
# This suggests a way to calculate $\pi$: if we have a square and the largest circle that
# fits inside that square, we can determine the ratio of areas of a circle and a square. We
# can calculate this ratio using a monte carlo simulation. We select random points inside a
# square, and we keep track of how often those points also fall inside the circle that fits
# perfectly inside that square.
#
# Given a large enough sampling points, $\frac{A_{circle}}{A_{square}}$ will be equal to the
# fraction of randomly chosen points inside the square that also fall inside the circle.
# Then we can figure out $\pi$!
#
# #### Pseudo-code
#
# Given the above, our algorithm for determining $\pi$ looks like this:
#
# 1. For each of $N$ iterations,
#     1. Select a random point inside a square of area $4r^2$ as Cartesian, $(x, y)$,
# coordinates.
#     1. Determine if the point also falls inside the circle embedded within this square of
# area $\pi r^2$.
#     1. Keep track of whether or not this point fell inside the circle. At the end of $N$
# iterations, you want to know $M$ -- the number of the $N$ random points that fell inside
# the circle!
# 1. Calculate $\pi$ as $4\frac{M}{N}$
#
# #### Exercise
#
# Write a function that calculates $\pi$ using Julia.
#
# The algorithm above should work for any value of $r$ that you choose to use. Make sure you
# make $N$ big enough that the value of $\pi$ is correct to at least a couple numbers after
# the decimal point!
#
# *Hint*:
#
# This will probably be easier if you center your circle and square at the coordinate (0, 0)
# and use a radius of 1. For example, to choose random coordinates within your square at
# position (x, y), you may want to choose x and y so that they are each a value between -1
# and +1. Then any point within a distance of 1 from (0, 0) will fall inside the circle!
#
#
#
#
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
global charNum = parse(Int, ARGS[1])

##
cd(@__DIR__)
using Pkg
Pkg.activate(".")
Pkg.instantiate()

##
using LinearAlgebra
using StatsBase
using CSV
using DataFrames
using Pipe
using ProgressBars
using ProgressMeter
using StatsPlots
using Base.Threads
using RCall
using StatsFuns
using MCPhylo
using Random
Random.seed!(9011452010850473240)

##

wals = CSV.read("../data/charMtx.csv", DataFrame)
d = CSV.read("../data/fpairMtx.csv", DataFrame)

famFreqs = combine(groupby(wals, :glot_fam), nrow)
sort!(famFreqs, :nrow, rev=true)

families = filter(x -> x.nrow>1, famFreqs).glot_fam
isolates = filter(x -> x.nrow==1, famFreqs).glot_fam

isoDict = @pipe wals |>
      filter(x -> x.glot_fam ∈ isolates, _) |>
      zip(_.glot_fam, _.longname) |>
      Dict

##
fm2trees = Dict()
@showprogress for fm in families
      fm2trees[fm] = MCPhylo.ParseNewick("../data/posteriorTrees/$fm.posterior.tree")
end

for fm in isolates
      fm2trees[fm] = repeat([Node(isoDict[fm])], 1000)
end

##

function renumberTree(tree::GeneralNode)
      tips = [nd for nd in post_order(tree) if nd.nchild==0]
      nonTips = [nd for nd in post_order(tree) if nd.nchild > 0]
      sort!(tips, lt= (x,y) -> x.name < y.name)
      for (i,nd) in enumerate(tips)
            nd.num = i
      end
      for (i,nd) in enumerate(nonTips)
            nd.num = i+length(tips)
      end
      tree
end

##

for fm in keys(fm2trees)
      fm2trees[fm] = renumberTree.(fm2trees[fm])
end

lineages = vcat(families, isolates)

nLineages = length(lineages)

ttrees = [[fm2trees[fm][i] for fm in lineages] for i in 1:1000]

taxa = Vector{String}[]
for t in ttrees[1]
      push!(taxa, sort([nd.name for nd in pre_order(t) if nd.nchild==0]))
end

nnodes = [length(post_order(t)) for t in ttrees[1]]

states = ["a", "b", "c", "d"]

nsites = 1

nbase = 4

rates = [1.0]

##

data = zeros(nbase, nsites, maximum(nnodes), length(ttrees[1]))
char = Dict(zip(d.taxon, d[:,charNum+1]))
for i in 1:length(taxa)
      for l in taxa[i]
            s = char[l]
            t_ind = find_by_name(ttrees[1][i],l).num
            if s in states
                  data[:,1,t_ind, i] = (states .== s)
            else
                  data[:,1,t_ind, i] .= 1
            end
      end
end

##

uSim = read(
      "modelFitting/output/universal_$(lpad(charNum, 2, "0")).jls",
      ModelChains,
)

samples = vcat([uSim.value[:,:,k] for k in 1:2]...)

nSamples = size(samples, 1)

yPred = zeros(nSamples, nLineages)

@threads for i in ProgressBar(1:nSamples)
      s_rates = samples[i,2:(end-1)]
      ind = samples[i,1]
      treeI = ceil(Int, 1000*invlogit(ind))
      fullSrates = exp.(s_rates)
      y_ = zeros(nLineages)
      for fmI = 1:nLineages
            y_[fmI] = logpdf(
                  PhyloDist(ttrees[treeI][fmI], fullSrates, [1.0], freeK),
                  data[:, :, :, fmI],
            )
      end
      yPred[i,:] = y_

end

looU = convert(
      Dict,
      R"
      library(loo)
      loo($yPred)
      ",
)

##

lSim = read(
      "modelFitting/output/lineage_$(lpad(charNum, 2, "0")).jls",
      ModelChains,
)

samples = vcat([lSim.value[:,:,k] for k in 1:2]...)
samples = samples[sample(1:size(samples,1),1000),:]

nSamples = size(samples, 1)

yPred = zeros(nSamples, nLineages)

@threads for i = ProgressBar(1:nSamples)
      s_rates_ = samples[i,2:(end-1)]
      s_rates = reshape(s_rates_, :, nLineages)
      ind = samples[i,1]
      treeI = ceil(Int, 1000 * invlogit(ind))
      fullSrates = exp.(s_rates)
      y_ = zeros(nLineages)
      for fmI = 1:nLineages
            y_[fmI] = logpdf(
                  PhyloDist(ttrees[treeI][fmI], fullSrates[:, fmI], [1.0], freeK),
                  data[:, :, :, fmI],
            )
      end
      yPred[i, :] = y_
end

looL = convert(
      Dict,
      R"
      library(loo)
      loo($yPred)
      ",
)
##

kys = ["looic", "se_looic", "elpd_loo", "se_elpd_loo"]

results = permutedims(
      hcat(
            vcat(
                  [
                        [[k, looU[k], "u"] for k in kys],
                        [[k, looL[k], "l"] for k in kys],
                  ]...,
            )...,
      ),
)

resultsDF = DataFrame(results, :auto)
rename!(resultsDF, [:parameter, :value, :model])
resultsDF[:,:char] .= charNum

try mkdir("looFamResults")
catch e
end

CSV.write("looFamResults/$(lpad(charNum, 2, "0")).csv", resultsDF)

max_stars_repo_name: erathorn/phylogeneticTypology

max_issues_repo_name: erathorn/phylogeneticTypology

max_forks_repo_name: erathorn/phylogeneticTypology


text:
# ============================================================================ #
# Rectangle
# ============================================================================ #

abstract type AbstractRectangle end

struct Rectangle{PointType} <: AbstractRectangle
    p1::PointType # lower_left
    p2::PointType # upper_right
end
function Rectangle(p1::T, p2::T) where {T <: AbstractPoint}
    return Rectangle{T}(p1, p2)
end
function Rectangle(x1::T, y1::T, x2::T, y2::T) where {T <: Number}
    return Rectangle( Point(x1, y1), Point(x2, y2) )
end

mutable struct MutableRectangle{PointType} <: AbstractRectangle
    p1::PointType
    p2::PointType
end

mutable struct BoundingBox{PointType}
    rectangle::MutableRectangle{PointType}
    midpoint::PointType
end

# ---------------------------------------------------------------------------- #
# Methods
# ---------------------------------------------------------------------------- #
Base.:(==)(r1::AbstractRectangle, r2::AbstractRectangle) = r1.p1 == r2.p1 && r1.p2 == r2.p2

function Base.:+(r::T, p::AbstractArray) where {T <: AbstractRectangle}
    q1 = r.p1 + p
    q2 = r.p2 + p
    return T(q1, q2)
end
function Base.:-(r::T, p::AbstractArray) where {T <: AbstractRectangle}
    q1 = r.p1 - p
    q2 = r.p2 - p
    return T(q1, q2)
end
function Base.:+(r::T) where {T <: AbstractRectangle}
    return r
end
function Base.:-(r::T) where {T <: AbstractRectangle}
    return T(-r.p1, -r.p2)
end
function Base.:+(p::AbstractArray, r::T) where {T <: AbstractRectangle}
    return r + p
end
function Base.:-(p::AbstractArray, r::T) where {T <: AbstractRectangle}
    return -(r - p)
end

function Base.:*(f::N, r::T) where {N <: Number, T <: AbstractRectangle}
    return T(f * r.p1, f * r.p2)
end
Base.:*(point::P, f::N) where {N <: Number, P <: AbstractPoint} = f * point;

function shift(rectangle::T, v) where {T <: AbstractRectangle}
    return T( rectangle.p1 + v, rectangle.p2 + v )
end
function scale(rectangle::T, factor) where {T <: AbstractRectangle}
    width  = factor * x_spread(rectangle)
    height = factor * y_spread(rectangle)
    mid    = mid_point(rectangle)

    PointType = typeof(rectangle.p1)
    p1 = PointType( mid.x - width/2, mid.y - height/2 )
    p2 = PointType( mid.x + width/2, mid.y + height/2 )

    return T(p1, p2)
end
function rotate(rectangle::T, rot_mat) where {T <: AbstractRectangle}
    mid = mid_point(rectangle)

    vm2 = rectangle.p2 - mid
    r

    x = rot_mat[1,1] * p.x + rot_mat[1,2] * p.y
    y = rot_mat[2,1] * p.x + rot_mat[2,2] * p.y
    return T(x, y)
end

function Base.getindex(r::T, index::Int) where {T <: AbstractRectangle}
    if index == 1
        return r.p1
    elseif index == 2
        return r.p2
    else
        return NaN
    end
end

function area(rectangle::AbstractRectangle)
    dx = rectangle.p2.x - rectangle.p1.x
    dy = rectangle.p2.y - rectangle.p1.y
    return dx * dy
end

function is_in_rectangle(rectangle::AbstractRectangle, point::AbstractPoint)
    bx = rectangle.p1.x <= point.x <= rectangle.p2.x
    if bx
        by = rectangle.p1.y <= point.y <= rectangle.p2.y
        return by
    else
        return false
    end
end

function centroid(rectangle::AbstractRectangle)
    return (rectangle.p1 + rectangle.p2)./2
end

function mid_point(rectangle::T) where {T <: AbstractRectangle}
    mid_x = 0.5 * (rectangle.p1.x + rectangle.p2.x)
    mid_y = 0.5 * (rectangle.p1.y + rectangle.p2.y)

    return Point(mid_x, mid_y)
end

function x_spread(rectangle::T) where {T <: AbstractRectangle}
    return rectangle.upper_right.x - rectangle.lower_left.x
end

function y_spread(rectangle::T) where {T <: AbstractRectangle}
    return rectangle.upper_right.y - rectangle.lower_left.y
end

function in_circle(rectangle::AbstractRectangle)
    mx = 0.5 * (rectangle.p1.x + rectangle.p2.x)
    my = 0.5 * (rectangle.p1.y + rectangle.p2.y)
    dx = rectangle.p2.x - rectangle.p1.x
    dy = rectangle.p2.y - rectangle.p1.y
    r = 0.5 * min(dx, dy)
    return Circle(Point(mx, my), r)
end

function out_circle(rectangle::AbstractRectangle)
    mx = 0.5 * (rectangle.p1.x + rectangle.p2.x)
    my = 0.5 * (rectangle.p1.y + rectangle.p2.y)
    dx = rectangle.p2.x - rectangle.p1.x
    dy = rectangle.p2.y - rectangle.p1.y
    R = sqrt(dx^2 + dy^2)
    return Circle(Point(mx, my), R)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using RLToolkit
using Test

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


text:
using SimpleSolve
using Base.Test
using StaticArrays

include("nlsolve.jl")
include("fixedpoint.jl")
include("seqsolve.jl")


text:
@testset "Lennard-Jones potential" begin
    @testset "Two interacting particles" begin

        m1 = 1.0
        m2 = 1.0
        r1 = 1.3

        t1 = 0.0
        t2 = 1.0
        τ = (t2 - t1) / 1000

        p1 = MassBody(SVector(-r1 / 2, 0.0, 0.0), SVector(0.0, 0.0, 0.0), m1)
        p2 = MassBody(SVector(r1 / 2, 0.0, 0.0), SVector(0.0, 0.0, 0.0), m2)

        σ = 1.0
        ϵ = 1.0
        parameters = LennardJonesParameters(ϵ, σ, Inf)
        system = PotentialNBodySystem([p1, p2], Dict(:lennard_jones => parameters))
        simulation = NBodySimulation(system, (t1, t2))
        sim_result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        r2 = get_position(sim_result, t2, 2) - get_position(sim_result, t2, 1)
        v_expected = sqrt(4ϵ / m1 * ( ((σ / norm(r1))^12 - (σ / norm(r2))^12) - ((σ / norm(r1))^6 - (σ / norm(r2))^6 ) ))
        v_actual = norm(get_velocity(sim_result, t2, 2))

        ε = 0.001 * v_expected
        @test v_expected ≈ v_actual atol = ε

        io = IOBuffer()

        @test sprint(io -> show(io, system)) ==
    "Potentials: \nLennard-Jones:\n\tϵ:1.0\n\tσ:1.0\n\tR:Inf\n"

        @test sprint(io -> show(io, simulation)) ==
    "Timespan: (0.0, 1.0)\nBoundary conditions: InfiniteBox{Float64}([-Inf, Inf, -Inf, Inf, -Inf, Inf])\nPotentials: \nLennard-Jones:\n\tϵ:1.0\n\tσ:1.0\n\tR:Inf\n"

        @test sprint(io -> show(io, sim_result)) ==
    "N: 2\nTimespan: (0.0, 1.0)\nBoundary conditions: InfiniteBox{Float64}([-Inf, Inf, -Inf, Inf, -Inf, Inf])\nPotentials: \nLennard-Jones:\n\tϵ:1.0\n\tσ:1.0\n\tR:Inf\nTime steps: 1000\nt: 0.0, 1.0\n"
    end

    @testset "Three particles of liquid argon and their \"temperature\"" begin
        T = 120.0 # °K
        kb = 8.3144598e-3 # kJ/(K*mol)
        ϵ = T * kb
        σ = 0.34 # nm
        ρ = 1374/1.6747# Da/nm^3
        m = 39.95# Da
        L = 5σ # 10.229σ
        N = 3 # floor(Int, ρ * L^3 / m)
        R = 2.25σ
        v_dev = sqrt(3 * kb * T / m)
        r1 = SVector(L / 3, L / 3, 2 * L / 3)
        r2 = SVector(L / 3, 2 * L / 3, L / 3)
        r3 = SVector(2 * L / 3, L / 3, L / 3)
        v1 = SVector(0, 0, -v_dev)
        v2 = SVector(0, -v_dev, 0)
        v3 = SVector(-v_dev, 0, 0)
        p1 = MassBody(r1, v1, m)
        p2 = MassBody(r2, v2, m)
        p3 = MassBody(r3, v3, m)

        τ = 0.5e-3 # ps or 1e-12 s
        t1 = 0.0
        t2 = 100τ

        parameters = LennardJonesParameters(ϵ, σ, R)
        lj_system = PotentialNBodySystem([p1, p2, p3], Dict(:lennard_jones => parameters));
        simulation = NBodySimulation(lj_system, (t1, t2), PeriodicBoundaryConditions(L), kb);
        result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        T1 = temperature(result, t1)
        ε = 1e-6
        @test T1 ≈ 120.0 atol = ε

        e_kin_1 = m * (dot(v1, v1) + dot(v2, v2) + dot(v3, v3)) / 2
        @test e_kin_1 == kinetic_energy(result, t1)

        e_tot_1 = total_energy(result, t1)
        ε = 0.1 * e_tot_1
        e_tot_2 = total_energy(result, t2)
        @test e_tot_1 ≈ e_tot_2 atol = ε

        simulation = NBodySimulation(lj_system, (t1, t2), CubicPeriodicBoundaryConditions(L));
        result = run_simulation(simulation, VelocityVerlet(), dt=τ)
        e_tot_1 = total_energy(result, t1)
        ε = 0.1 * e_tot_1
        e_tot_2 = total_energy(result, t2)
        @test e_tot_1 ≈ e_tot_2 atol = ε

        io = IOBuffer()
        pdb_data = sprint(io -> NBodySimulator.write_pdb_data(io, result))
        splitted_data = split(pdb_data, '\n')

        hetatm_count = 0
        timestep_count = 0

        for s in splitted_data
            if length(s) >= 10 && s[1:10] == "REMARK 250"
                timestep_count += 1
            elseif length(s) >= 6 && s[1:6] == "HETATM"
                hetatm_count += 1
            end
        end

        molecule_number = hetatm_count / timestep_count

        @test length(result.solution.t) == timestep_count
        @test molecule_number == length(lj_system.bodies)

    end

    @testset "Testing RDF and MSD calculation" begin
        T = 120.0 # °K
        kb = 8.3144598e-3 # kJ/(K*mol)
        ϵ = T * kb
        σ = 0.34 # nm
        ρ = 1374/1.6747# Da/nm^3
        m = 39.95# Da
        L = 5σ # 10.229σ
        N = 125 # floor(Int, ρ * L^3 / m)
        R = 0.5L
        τ = 0.5e-3 # ps or 1e-12 s
        t1 = 0.0
        t2 = 400τ
        v_dev = sqrt(kb * T / m)
        bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)

        parameters = LennardJonesParameters(ϵ, σ, R)
        lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => parameters));
        simulation = NBodySimulation(lj_system, (t1, t2), CubicPeriodicBoundaryConditions(L), kb);
        result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        (ts, mean_square_displacement) = msd(result)
        @test mean_square_displacement[1] < mean_square_displacement[end]

        (rs, grs) = rdf(result)
        (val, ind) = findmax(grs)
        @test (rs[ind] / σ) ≈ 1.0 atol = 1.0
    end

    @testset "Constructing electorstatic potential parameters entity" begin
        default_potential = LennardJonesParameters()
        @test 1.0 == default_potential.ϵ
        @test 1.0 == default_potential.σ
        @test 2.5 == default_potential.R
        @test 1.0 == default_potential.σ2
        @test 6.25 == default_potential.R2

        io = IOBuffer()
        potential1 = LennardJonesParameters()
        potential2 = LennardJonesParameters(2, 5, 10)
        @test sprint(io -> show(io, potential1)) == "Lennard-Jones:\n\tϵ:1.0\n\tσ:1.0\n\tR:2.5\n"
        @test sprint(io -> show(io, potential2)) == "Lennard-Jones:\n\tϵ:2\n\tσ:5\n\tR:10\n"
    end

    @testset "Periodic boundar conditions" begin
        pbc = PeriodicBoundaryConditions(15e-6)
        boundary = [0.0, 15e-6, 0.0, 15e-6, 0.0, 15e-6]
        ind = 1
        for b in pbc
            b = boundary[ind]
            ind += 1
        end
        for i = 1:6
            @test boundary[i] == pbc[i]
        end
    end
end


text:
struct AlignParam
   mismatches::Int      # Allowable mismatches
   kmer_size::Int       # Minimum number of matches to extend past an edge
   seed_try::Int        # Starting number of seeds
   seed_tolerate::Int   # Allow at most _ hits for a valid seed
   seed_min_qual::Int   # Absolute lowest quality score to allow for worst base in seed
   seed_length::Int     # Seed Size
   seed_buffer::Int     # Ignore first and last _ bases
   seed_inc::Int        # Incrementation for subsequent seed searches
   seed_pair_range::Int # Seeds for paired-end reads must match within _ nt of each other
   score_range::Float64 # Identity scoring range to be considered repetitive
   score_min::Float64   # Minimum percent identity score for a valid alignment
   is_stranded::Bool    # Is input data + strand only?
   is_paired::Bool      # Paired end data?
   is_pair_rc::Bool     # Is mate 2 the reverse complement of mate 1
   is_trans_ok::Bool    # Do we allow edges from one gene to another
   is_circ_ok::Bool     # Do we allow back edges
end

AlignParam( ispaired = false ) = ispaired ? AlignParam( 3, 9, 2, 4, 4, 18, 5, 18, 2500, 0.05, 0.7,
                                                        false, true, true, false, true ) :
                                            AlignParam( 3, 9, 2, 4, 4, 18, 5, 18, 1000, 0.05, 0.7,
                                                        false, false, true, false, true )

# load from command line args
@inline function AlignParam( args::Dict{S,Any}, ispaired=false; kmer=9 ) where S <: AbstractString
   aln = AlignParam( args["mismatches"], kmer,
                           args["seed-try"],
                           args["seed-tol"], 4,
                           args["seed-len"],
                           args["seed-buf"],
                           args["seed-inc"],
                           args["pair-range"], 0.05,
                           args["score-min"],
                           args["stranded"], ispaired,
                          !args["pair-same-strand"], false,
                           args["circ"] )
   aln
end

abstract type UngappedAlignment end

const MatchesInt    = UInt8
const MismatchInt   = UInt8
const MismatchFloat = Float32

mutable struct SGAlignScore
   matches::MatchesInt
   mismatches::MismatchInt
   mistolerance::MismatchFloat
end

@inline Base.zero(t::Type{SGAlignScore}) = SGAlignScore(0x0000,zero(MismatchInt),zero(MismatchFloat))
@inline Base.one(t::Type{SGAlignScore}) = SGAlignScore(0x0001,one(MismatchInt),one(MismatchFloat))
@inline match(t::Type{SGAlignScore}, val::I) where {I <: Integer} = SGAlignScore(convert(MatchesInt, val),zero(MismatchInt),zero(MismatchFloat))

struct SGAlignNode
   gene::NodeInt
   node::NodeInt
   score::SGAlignScore
end

Base.hash( sgn::SGAlignNode ) = hash( sgn.gene, hash(sgn.node) )
Base.hash( sgn::SGAlignNode, h::UInt64 ) = hash( sgn.gene, hash(sgn.node, h) )
Base.isequal( a::SGAlignNode, b::SGAlignNode ) = a.gene == b.gene && a.node == b.node
Base.:(==)( a::SGAlignNode, b::SGAlignNode ) = a.gene == b.gene && a.node == b.node

const SGAlignPath = Vector{SGAlignNode}

@inline SGAlignNode(sgn::SGNode) = SGAlignNode(sgn.gene, sgn.node, zero(SGAlignScore))
@inline Base.push!( v::Vector{SGAlignNode}, sgn::SGNode ) = push!( v, SGAlignNode(sgn) )

const ReadLengthInt = UInt8

mutable struct SGAlignment <: UngappedAlignment
   offset::UInt32
   offsetread::ReadLengthInt
   path::Vector{SGAlignNode}
   strand::Bool
   isvalid::Bool
end

SGAlignment() = SGAlignment(zero(UInt32), zero(ReadLengthInt), SGAlignNode[], true, false)

const SGAlignVec = Nullable{Vector{SGAlignment}}

const DEF_ALIGN = SGAlignment(zero(UInt32), zero(ReadLengthInt), SGAlignNode[], true, false)

@inbounds function matches( v::Vector{SGAlignNode} )
   matches = zero(MatchesInt)
   for i in 1:length(v)
      matches += v[i].score.matches
   end
   matches
end
@inbounds function mistolerance( v::Vector{SGAlignNode} )
   mistol = zero(MismatchFloat)
   for i in 1:length(v)
      mistol += v[i].score.mistolerance
   end
   mistol
end
@inbounds function score( v::Vector{SGAlignNode} )
   score = zero(MismatchFloat)
   for i in 1:length(v)
      score += v[i].score.matches - v[i].score.mistolerance
   end
   score
end
@inline matches( aln::A ) where {A <: UngappedAlignment} = matches( aln.path )
@inline mistolerance( aln::A ) where {A <: UngappedAlignment} = mistolerance( aln.path )
@inline score( aln::A ) where {A <: UngappedAlignment} = score( aln.path )
@inline identity( align::A, readlen::I ) where {A <: UngappedAlignment, I <: Integer} = convert(Float32, @fastmath score( align ) / readlen)
@inline isvalid( align::A ) where {A <: UngappedAlignment} = align.isvalid && length(align.path) >= 1 && matches(align) > mistolerance(align) ? true : false

Base.:>( a::SGAlignment, b::SGAlignment ) = >( score(a), score(b) )
Base.:<( a::SGAlignment, b::SGAlignment ) = <( score(a), score(b) )
Base.:>=( a::SGAlignment, b::SGAlignment ) = >=( score(a), score(b) )
Base.:<=( a::SGAlignment, b::SGAlignment ) = <=( score(a), score(b) )
Base.isless( a::SGAlignment, b::SGAlignment ) = a < b

# add prob of being accurate base to mismatch, rather than integer.
@inline phred_to_prob( phred::UInt8 ) = convert(Float32, @fastmath 1-10^(-convert(Int8,phred)/10))

@inline function Base.empty!( align::SGAlignment )
   align.offset = zero(UInt32)
   align.offsetread = zero(ReadLengthInt)
   empty!( align.path )
   align.strand = true
   align.isvalid = false
end

@inline function remove_invalid!( vec::Vector{SGAlignment} )
   i = 1
   while i <= length(vec)
      @inbounds if !isvalid( vec[i] )
         deleteat!( vec, i )
      else
         i += 1
      end
   end
end

@inline function seed_locate( p::AlignParam, index::FMIndex, read::FASTQRecord; offset_left::Bool=true, both_strands::Bool=true )
   def_sa = 2:1
   readlen = convert(ReadLengthInt, length(read.sequence))
   ctry = 1
   if offset_left
      increment = p.seed_inc
      increment_sm = 1
      curpos = p.seed_buffer
   else
      increment = -p.seed_inc
      increment_sm = -1
      curpos = readlen - p.seed_length - p.seed_buffer
   end
   maxright = readlen - p.seed_length - p.seed_buffer
   while( ctry <= p.seed_try && p.seed_buffer <= curpos <= maxright )
      if minimum( read.quality[curpos:(curpos+p.seed_length-1)] ) <= p.seed_min_qual
         curpos += increment_sm
         continue
      end
      curseq = read.sequence[curpos:(curpos+p.seed_length-1)]
      sa = FMIndexes.sa_range( curseq, index )
      cnt = length(sa)
      ctry += 1
      if cnt == 0 || cnt > p.seed_tolerate
         if both_strands
            reverse_complement!(curseq)
            rev_sa = FMIndexes.sa_range( curseq, index )
            rev_cnt = length(rev_sa)
            if 0 < rev_cnt <= p.seed_tolerate
               rev_pos = readlen - (curpos+p.seed_length-1)+1
               return rev_sa,rev_pos,false
            end
         end
      else
         return sa,curpos,true
      end
      curpos += increment
   end
   def_sa,curpos,true
end

@inline function splice_by_score!( arr::Vector{A}, threshold, buffer ) where A <: UngappedAlignment
   i = 1
   while i <= length( arr )
      @inbounds if threshold - score( arr[i] ) > buffer
         splice!( arr, i )
         i -= 1
      end
      i += 1
   end
end

@inline function splice_by_identity!( arr::Vector{A}, threshold, buffer, readlen ) where A <: UngappedAlignment
   i = 1
   while i <= length( arr )
      @inbounds if threshold - identity( arr[i], readlen ) > buffer
         splice!( arr, i )
         i -= 1
      end
      i += 1
   end
end

@inline function _ungapped_align( p::AlignParam, lib::GraphLib, read::FASTQRecord, indx::Int, readloc::Int;
                                  ispos::Bool=true, geneind::NodeInt=convert(NodeInt,searchsortedlast( lib.offset, indx )) )

   sgidx       = indx - lib.offset[geneind]
   offset_node = convert(NodeInt,searchsortedlast(lib.graphs[geneind].nodeoffset,convert(CoordInt, sgidx + 1)))

   align = ungapped_fwd_extend( p, lib, geneind, sgidx + 1,
                                read, readloc + 1, ispos=ispos,
                                nodeidx=offset_node )

   align = ungapped_rev_extend( p, lib, geneind, sgidx,
                                read, readloc, ispos=ispos, align=align,
                                nodeidx=offset_node )
   align
end

@inline function ungapped_align( p::AlignParam, lib::GraphLib, read::FASTQRecord;
                         ispos::Bool=true, anchor_left::Bool=true )

   seed,readloc,pos = seed_locate( p, lib.index, read, offset_left=anchor_left, both_strands=!p.is_stranded )
   readlen = convert(ReadLengthInt, length(read.sequence))

   if !pos
      reverse_complement!( read )
      ispos = false
   end

   res      = Nullable{Vector{SGAlignment}}()
   maxscore = zero(Float32)

   for sidx in 1:length(seed)
      s = FMIndexes.sa_value( seed[sidx], lib.index ) + 1

      align = _ungapped_align( p, lib, read, s, readloc, ispos=ispos )
      scvar = identity(align, readlen)

      if isvalid(align)
         if isnull( res )
            res = Nullable(SGAlignment[ align ])
            maxscore = scvar
         else
            # new best score
            if scvar > maxscore
               # better than threshold for all of the previous scores
               if scvar - maxscore > p.score_range
                  length( res.value ) >= 1 && empty!( res.value )
               else # keep at least one previous score
                  splice_by_identity!( res.value, scvar, p.score_range, readlen )
               end
               maxscore = scvar
               push!( res.value, align )
            else # new score is lower than previously seen
               if maxscore - scvar <= p.score_range # but tolerable
                  push!( res.value, align )
               end
            end # end score vs maxscore
         end # end isnull
      end # end isvalid

   end
   res
end

# This is the main ungapped alignment extension function in the --> direction
# Returns: SGAlignment
@inbounds function ungapped_fwd_extend( p::AlignParam, lib::GraphLib, geneind::NodeInt, sgidx::Int,
                                read::FASTQRecord, ridx::Int; ispos::Bool=true,
                                align::SGAlignment=SGAlignment(sgidx, ridx, SGAlignNode[],ispos,false),
                                nodeidx::NodeInt=convert(NodeInt,searchsortedlast(lib.graphs[geneind].nodeoffset,CoordInt(sgidx))) )
   const sg       = lib.graphs[geneind]
   const readlen  = convert(ReadLengthInt, length(read.sequence))

   passed_edges  = Nullable{Vector{UInt8}}() # don't allocate array unless needed
   mistol        = mistolerance(align)

   push!( align.path, SGAlignNode( geneind, nodeidx, zero(SGAlignScore) ) ) # starting node

   @inline function spliced_fwd_extend( geneind::CoordInt, edgeind::UInt32,
                                        ridx::Int, right_kmer_ind::Int, align::SGAlignment )
      # Choose extending node through intersection of lib.edges.left ∩ lib.edges.right
      # returns right nodes with matching genes
      isassigned( lib.edges.right, right_kmer_ind ) || return align

      best = align
      const left_kmer_ind = kmer_index(lib.graphs[geneind].edgeleft[edgeind])
      isassigned( lib.edges.left, left_kmer_ind ) || return align
      @inbounds const right_nodes = lib.edges.left[ left_kmer_ind ] ∩ lib.edges.right[ right_kmer_ind ]

      for rn in right_nodes
         rn.gene == align.path[1].gene || continue
         const rn_offset = lib.graphs[rn.gene].nodeoffset[rn.node]

         res_align::SGAlignment = ungapped_fwd_extend( p, lib, rn.gene, Int(rn_offset), read, ridx,
                                                       align=deepcopy(align), nodeidx=rn.node )
         #println(res_align)
         best = res_align > best ? res_align : best
      end
      if score(best) >= score(align) + p.kmer_size
         best.isvalid = true
      end
      best
   end

   while( mistol <= p.mismatches && ridx <= readlen && sgidx <= length(sg.seq) )
      if sgidx == sg.nodeoffset[nodeidx] + sg.nodelen[nodeidx] # hit next edge
         const curedge = nodeidx + 0x01
         if     sg.edgetype[curedge] == EDGETYPE_LR &&
                sg.nodelen[curedge]  >= p.kmer_size && # 'LR' && nodelen >= K
                readlen - ridx + 1   >= p.kmer_size
               # check edgeright[curnode+1] == read[ nextkmer ]
               # move forward K, continue
               # This is a std exon-exon junction, if the read is an inclusion read
               # then we simply jump ahead, otherwise lets try to find a compatible right
               # node
               const rkmer_ind = kmer_index(read.sequence[ridx:(ridx+p.kmer_size-1)])
               rkmer_ind == 0 && break
               if kmer_index(sg.edgeright[curedge]) == rkmer_ind
                  ridx    += p.kmer_size - 1
                  sgidx   += p.kmer_size - 1
                  nodeidx += 0x01
                  push!( align.path, SGAlignNode( geneind, nodeidx, match(SGAlignScore, p.kmer_size - 1) ) )
                  align.isvalid = true
               else
                  align = spliced_fwd_extend( geneind, curedge, ridx, rkmer_ind, align )
                  break
               end
         elseif sg.edgetype[curedge] == EDGETYPE_LR ||
                sg.edgetype[curedge] == EDGETYPE_LL # 'LR' || 'LL'
               # if length(edges) > 0, push! edgematch then set to 0
               # push! edge,  (here we try to extend fwd, but keep track
               # of potential edges we pass along the way.  When the alignment
               # dies if we have not had sufficient matches we then explore
               # those edges.
               if isnull(passed_edges) # now we have to malloc
                  passed_edges = Nullable(Vector{UInt8}())
               end
               push!(get(passed_edges), length(align.path))
               nodeidx += 0x01
               push!( align.path, SGAlignNode( geneind, nodeidx, zero(SGAlignScore) ) )
         elseif sg.edgetype[curedge] == EDGETYPE_LS # 'LS'
               if readlen - ridx + 1   >= p.kmer_size
                  # obligate spliced_extension
                  const rkmer_ind = kmer_index(read.sequence[ridx:(ridx+p.kmer_size-1)])
                  rkmer_ind == 0 && break
                  align = spliced_fwd_extend( geneind, curedge, ridx, rkmer_ind, align )
               end
               break
         elseif sg.edgetype[curedge] == EDGETYPE_SR
               # end of alignment
               break # ?
         else #'RR' || 'SL' || 'RS'
               # ignore 'RR' and 'SL'
               nodeidx += 0x01
               nodeidx <= length(sg.nodelen) && push!( align.path, SGAlignNode( geneind, nodeidx, zero(SGAlignScore) ) )
         end
         # ignore 'N'
      end
      if read.sequence[ridx] == sg.seq[sgidx]
         # match
         align.path[end].score.matches += 0x01
      else
         # mismatch
         const prob  = phred_to_prob( read.quality[ridx] )
         @fastmath align.path[end].score.mistolerance += prob
         @fastmath mistol += prob
         align.path[end].score.mismatches += 0x01
      end
      ridx  += 1
      sgidx += 1
   end

   # if edgemat < K, spliced_extension for each in length(edges)
   if !isnull(passed_edges)
         # go back.
      for c in length(get(passed_edges)):-1:1  #most recent edge first
         const cval = get(passed_edges)[c]
         if matches(align.path[cval+1:end]) >= p.kmer_size + p.mismatches
            align.isvalid = true
            break
         else
            const cur_edge = align.path[cval].node + one(UInt32)
            splice!( align.path, (cval+1):length(align.path) ) # extension into current node failed

            const cur_ridx = ridx - (sgidx - sg.nodeoffset[ cur_edge ])
            (cur_ridx + p.kmer_size - 1) <= readlen || continue

            const rkmer_ind  = kmer_index(read.sequence[cur_ridx:(cur_ridx+p.kmer_size-1)])
            rkmer_ind == 0 && continue
            res_align = spliced_fwd_extend( geneind, cur_edge,
                                            cur_ridx, rkmer_ind, align )
            if length(res_align.path) > length(align.path) # we added a valid node
               align = res_align > align ? res_align : align
            end
         end
      end
   end

   # clean up any bad trailing nodes
   while length(align.path) >= 1 && (align.path[end].score.mismatches >= align.path[end].score.matches ||
                                     Int(align.path[end].score.matches) - Int(align.path[end].score.mismatches) <= 1)
      pop!( align.path )
      length(align.path) == 0 && (align.isvalid = false)
   end

   if identity(align, readlen) >= p.score_min
      align.isvalid = true
   end
   align
end

# This is the main ungapped alignment extension function in the <-- direction
# Returns: SGAlignment
@inbounds function ungapped_rev_extend( p::AlignParam, lib::GraphLib, geneind::NodeInt, sgidx::Int,
                                read::FASTQRecord, ridx::Int; ispos::Bool=true,
                                align::SGAlignment=SGAlignment(sgidx, ridx, SGAlignNode[], ispos, false),
                                nodeidx::NodeInt=convert(NodeInt,searchsortedlast(lib.graphs[geneind].nodeoffset,CoordInt(sgidx))) )
   const sg       = lib.graphs[geneind]
   const readlen  = convert(ReadLengthInt,length(read.sequence))

   passed_edges  = Nullable{Vector{UInt8}}() # don't allocate array unless needed
   mistol        = mistolerance(align)

   if length(align.path) <= 0 || align.path[1].gene != geneind || align.path[1].node != nodeidx
      unshift!( align.path, SGAlignNode( geneind, nodeidx, zero(SGAlignScore) ) ) # starting node if not already there
   end

   @inline function spliced_rev_extend( geneind::NodeInt, edgeind::CoordInt,
                                        ridx::Int, left_kmer_index::Int, align::SGAlignment )
      # Choose extending node through intersection of lib.edges.left ∩ lib.edges.right
      # returns right nodes with matching genes
      isassigned( lib.edges.left, left_kmer_index ) || return align

      best = align
      const right_kmer_ind = kmer_index(lib.graphs[geneind].edgeright[edgeind])
      isassigned( lib.edges.right, right_kmer_ind ) || return align
      const left_nodes = lib.edges.right[ right_kmer_ind ] ∩ lib.edges.left[ left_kmer_index ]

      # do a test for trans-splicing, and reset left_nodes
      for rn in left_nodes
         rn.gene == align.path[1].gene || continue
         const rn_offset = lib.graphs[rn.gene].nodeoffset[rn.node-1] + lib.graphs[rn.gene].nodelen[rn.node-1] - 1

         res_align::SGAlignment = ungapped_rev_extend( p, lib, rn.gene, Int(rn_offset), read, ridx,
                                                       align=deepcopy(align), nodeidx=rn.node-0x01 )
         best = res_align > best ? res_align : best
      end
      if score(best) >= score(align) + p.kmer_size
         best.isvalid = true
      end
      best
   end

   while( mistol <= p.mismatches && ridx > 0 && sgidx > 0 )
      if sgidx == sg.nodeoffset[nodeidx] - 1 # hit right edge
         const leftnode = nodeidx - 0x01
         if     sg.edgetype[nodeidx] == EDGETYPE_LR &&
                sg.nodelen[leftnode] >= p.kmer_size && # 'LR' && nodelen >= K
                                ridx >= p.kmer_size

               const lkmer_ind  = kmer_index(read.sequence[(ridx-p.kmer_size+1):ridx])
               lkmer_ind == 0 && break
               if kmer_index(sg.edgeleft[nodeidx]) == lkmer_ind
                  ridx    -= p.kmer_size - 1
                  sgidx   -= p.kmer_size - 1
                  nodeidx -= 0x01
                  unshift!( align.path, SGAlignNode( geneind, nodeidx, match(SGAlignScore, p.kmer_size - 1) ) )
                  align.isvalid = true
               else
                  align = spliced_rev_extend( geneind, nodeidx, ridx, lkmer_ind, align )
                  break
               end
         elseif sg.edgetype[nodeidx] == EDGETYPE_LR ||
                sg.edgetype[nodeidx] == EDGETYPE_RR # 'LR' || 'RR'

               if isnull(passed_edges) # now we have to malloc
                  passed_edges = Nullable(Vector{UInt8}())
               end
               push!(get(passed_edges), length(align.path)-1)
               nodeidx -= 0x01
               unshift!( align.path, SGAlignNode( geneind, nodeidx, zero(SGAlignScore) ) )
         elseif sg.edgetype[nodeidx] == EDGETYPE_SR # 'SR'
               if ridx >= p.kmer_size
                  # obligate spliced_extension
                  const lkmer_ind  = kmer_index(read.sequence[(ridx-p.kmer_size+1):ridx])
                  lkmer_ind == 0 && break
                  align = spliced_rev_extend( geneind, nodeidx, ridx, lkmer_ind, align )
               end
               break
         elseif sg.edgetype[nodeidx] == EDGETYPE_LS
               # end of alignment
               break # ?
         else #'LL' || 'RS' || 'SL'
               # ignore 'LL' and 'RS'
               nodeidx -= 0x01
               nodeidx > 0 && unshift!( align.path, SGAlignNode( geneind, nodeidx, zero(SGAlignScore) ) )
         end
         # ignore 'N'
      end
      if read.sequence[ridx] == sg.seq[sgidx]
         # match
         align.path[1].score.matches += 0x01
         # mismatch
      else
         const prob  = phred_to_prob( read.quality[ridx] )
         @fastmath align.path[1].score.mistolerance += prob
         @fastmath mistol += prob
         align.path[1].score.mismatches += 0x01
      end
      ridx  -= 1
      sgidx -= 1
   end

   cur_ridx  = ridx + 1
   cur_sgidx = sgidx + 1
   # if passed_extend < K, spliced_extension for each in length(edges)
   if !isnull(passed_edges)
         # go back.
      for c in length(get(passed_edges)):-1:1  #most recent edge first
         const cval = get(passed_edges)[c]
         if matches(align.path[1:(end-(cval+1))]) >= p.kmer_size + p.mismatches
            align.isvalid = true
            break
         else
            const cur_edge = align.path[end-cval].node
            splice!( align.path, 1:(length(align.path)-(cval+1)) ) # extension into current node failed
            cur_ridx  = ridx + (sg.nodeoffset[ cur_edge ] - sgidx)
            cur_sgidx = sg.nodeoffset[ cur_edge ]
            (cur_ridx - p.kmer_size) > 0 || continue

            const lkmer_ind  = kmer_index(read.sequence[(cur_ridx-p.kmer_size):(cur_ridx-1)])
            lkmer_ind == 0 && continue
            res_align = spliced_rev_extend( geneind, cur_edge,
                                            cur_ridx-1, lkmer_ind, align )
            if length(res_align.path) > length(align.path) # we added a valid node
               align = res_align > align ? res_align : align
            end
         end
      end
   end

   if identity(align, readlen) >= p.score_min
      align.isvalid = true
   end
   if sgidx < align.offset
      align.offset = convert(UInt32, cur_sgidx)
   end
   if cur_ridx < align.offsetread
      align.offsetread = convert(ReadLengthInt, cur_ridx)
   end

   # clean up any bad left trailing nodes
   while length(align.path) >= 1 && (align.path[1].score.mismatches >= align.path[1].score.matches ||
                                     Int(align.path[1].score.matches) - Int(align.path[1].score.mismatches) <= 1)
      offset_change = align.path[1].score.matches + align.path[1].score.mismatches
      shift!( align.path )
      if length(align.path) >= 1
         align.offset      = sg.nodeoffset[ align.path[1].node ]
         align.offsetread += offset_change
      else
         align.isvalid = false
      end
   end

   align
end


text:
module ConwayTest

using Base.Test

using CellularAutomata
export blinker_test, turbine_test

# Function to test Conway's blinker oscillator and toad oscillator
function blinker_test()
  init = Array{Int}(zeros(10,10))

  # 'Toad Oscillator'
  init[6, 4] = 1
  init[6, 5] = 1
  init[6, 6] = 1

  init[7, 5] = 1
  init[7, 6] = 1
  init[7, 7] = 1

  ca = CA2d([3], [2,3], init, 9);

  @test ca.cells[:, :, 3] == ca.cells[:, :, 1]
  @test ca.cells[:, :, 5] == ca.cells[:, :, 3]
  @test ca.cells[:, :, 4] == ca.cells[:, :, 2]
  @test ca.cells[:, :, 8] == ca.cells[:, :, 6]

  # 'Blinker Oscillator'
  init[6, 4] = 0
  init[6, 5] = 0
  init[6, 6] = 0

  ca = CA2d([3], [2,3], init, 9)

  @test ca.cells[:, :, 3] == ca.cells[:, :, 1]
  @test ca.cells[:, :, 5] == ca.cells[:, :, 3]
  @test ca.cells[:, :, 4] == ca.cells[:, :, 2]
  @test ca.cells[:, :, 8] == ca.cells[:, :, 6]
end

# Function to test Conway's turbine automota which repeats every ninth step
function turbine_test()
  init = Array{Int}(zeros(20,20))
  init[4, 4:9] = 1
  init[5, 4:9] = 1

  init[4:9, 11] = 1
  init[4:9, 12] = 1

  init[11, 7:12] = 1
  init[12, 7:12] = 1

  init[7:12, 4] = 1
  init[7:12, 5] = 1

  ca = CA2d([3], [2,3], init, 18)

  for i in 9:18
    @test ca.cells[:, :, i] == ca.cells[:, :, i-8]
  end
end
end


text:
import ImageFiltering: imfilter, KernelFactors.gaussian

"""
    GaussianBlur <: ColorOperation

Description
--------------

Blurs an image using a Gaussian filter.

Usage
--------------

    GaussianBlur(k, [σ])

Arguments
--------------

- **`k`** : `Integer` or `AbstractVector` of `Integer` that denote
    the kernel size. It must be an odd positive number.
- **`σ`** : Optional. `Real` or `AbstractVector` of `Real` that denote the
    standard deviation. It must be a positive number.
    Defaults to `0.3 * ((k - 1) / 2 - 1) + 0.8`.

Examples
--------------

```
using Augmentor
img = testpattern()

# use exactly k=3 and σ=1.0
augment(img, GaussianBlur(3, 1.0))

# pick k and σ randomly from the specified ranges
augment(img, GaussianBlur(3:2:7, 1.0:0.1:2.0))
```
"""
struct GaussianBlur{K <: AbstractVector, S <: AbstractVector} <: ColorOperation
    k::K
    σ::S

    function GaussianBlur(k::K, σ::S) where {K <: AbstractVector{<:Integer},
                                             S <: AbstractVector{<:Real}}
        minimum(k) > 0 || throw(ArgumentError("Kernel size must be positive: $(k)"))
        minimum(σ) > 0 || throw(ArgumentError("σ must be positive: $(σ)"))
        new{K, S}(k, σ)
    end
end

# The default value for σ is taken from Albumentations
GaussianBlur(k) = GaussianBlur(k, 0.3 * ((k - 1) / 2 - 1) + 0.8)
GaussianBlur(k, σ) = GaussianBlur(vectorize(k), vectorize(σ))

randparam(op::GaussianBlur, img) = (safe_rand(op.k), safe_rand(op.σ))

@inline supports_eager(::Type{<:GaussianBlur}) = true

function applyeager(op::GaussianBlur, img::AbstractArray, (k, σ))
    n = ndims(img)
    kernel = gaussian(ntuple(_->σ, n), ntuple(_->k, n))
    return imfilter(img, kernel)
end

function showconstruction(io::IO, op::GaussianBlur)
    print(io, typeof(op).name.name, '(', op.k, ", ", op.σ,')')
end

function Base.show(io::IO, op::GaussianBlur)
    if get(io, :compact, false)
        print(io, "GaussianBlur with k=$(op.k) and σ=$(op.σ)")
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end


text:
# functions to download the data from Schneider's database

const SCHNEIDER_API = "https://data.exchange.se.com/explore/"
const DATASET = "dataset/microgrid-energy-management-benchmark-time-series"

function download_site_csv_from_schneider(siteid::Int,
                           path_to_data_folder::String;
                           compressed::Bool = true,
                           file_size = nothing,
                           periods::Union{Nothing, AbstractArray{Int}} = nothing,
                           progress = true, # true/false or Progress or ParallelProgress
                           api_key = get(ENV, "SCHNEIDER_API_KEY", nothing))

    if progress == true
        file_size = get_file_size(siteid, compressed)
        progress = Progress(file_size; desc = "Downloading $siteid.csv.gz ")
    elseif progress == false
        file_size = nothing
        progress = nothing
    end

    isnothing(api_key) && error("""you did not provide your api key
        please set it with the keyword `api_key`
        or with the ENV veriable: `ENV["SCHNEIDER_API_KEY"] = *your api key*`""")

    headers = Dict("Authorization" => "Apikey "*api_key)

    if compressed
        headers["Accept-Encoding"] = "gzip"
    end

    url = SCHNEIDER_API*DATASET*"/download/?format=csv&refine.site_id=$(siteid)"

    if !isnothing(periods)
        for p in periods
            url *= "&refine.period_id=$p"
        end
        url *= "&disjunctive.period_id=true"
    end

    url *= "&use_labels_for_header=true&csv_separator=%3B"

    file_extension = compressed ? ".csv.gz" : ".csv"
    _download(url,
              joinpath(path_to_data_folder, "$(siteid)"*file_extension);
              headers=headers,
              file_size=file_size,
              progress=progress,
              )
end


text:
using Base.Test
using Compat

import AffineTransforms

# Pure translations
@testset "Pure $(nd)D translations" for nd = 1:4
    t1 = rand(nd)
    t2 = rand(nd)
    ident = eye(nd)

    tfm1 = AffineTransforms.AffineTransform(ident, t1)
    x = zeros(nd)
    @test AffineTransforms.tformfwd(tfm1, x) == t1
    @test_throws DimensionMismatch AffineTransforms.tformfwd(tfm1, zeros(nd+1))
    @test AffineTransforms.tforminv(tfm1, x) == -t1
    x = rand(nd)
    @test AffineTransforms.tformfwd(tfm1, x) == x+t1
    @test AffineTransforms.tforminv(tfm1, x) == x-t1
    @test tfm1*x == x+t1
    @test tfm1\x == x-t1
    X = reshape(x, nd, 1)
    @test all(AffineTransforms.tformfwd(tfm1, X) .== x+t1)
    @test all(AffineTransforms.tforminv(tfm1, X) .== x-t1)

    tfm2 = AffineTransforms.tformtranslate(t2)
    tfmP1 = tfm1*tfm2
    @test tfmP1.offset == t1+t2
    @test AffineTransforms.tformfwd(tfmP1, x) ≈ x+t2+t1

    if 2 <= nd <= 3
        # Tuple interface
        r = AffineTransforms.tformfwd(tfm1, x...)
        @test [r...] == tfm1*x
        r = AffineTransforms.tforminv(tfm1, x...)
        @test [r...] == tfm1\x
    end
end

@testset "Pure rotations" begin
    @test AffineTransforms.tformrotate(0).scalefwd ≈ eye(2)
    @test AffineTransforms.tformrotate(pi/2)*[1,0] ≈ [0,1]
    @test AffineTransforms.tformrotate(pi/2).scalefwd ≈ [0 -1; 1 0]
    @test AffineTransforms.tformrotate(pi)*[1,0] ≈ [-1,0]
    @test AffineTransforms.tformrotate(3pi/2)*[1,0] ≈ [0,-1]
    @test AffineTransforms.tformrotate([1,0,0],pi/2).scalefwd ≈ [1 0 0; 0 0 -1; 0 1 0]
    @test AffineTransforms.tformrotate([1,0,0],pi/2)*[1,0,0] ≈ [1,0,0]
    @test AffineTransforms.tformrotate([0,1,0],pi/2)*[1,0,0] ≈ [0,0,-1]
    @test AffineTransforms.tformrotate([0,0,1],pi/2)*[1,0,0] ≈ [0,1,0]
    s2 = 1/sqrt(2)
    @test AffineTransforms.tformrotate([0,1,0],pi/4).scalefwd ≈ [s2 0 s2; 0 1 0; -s2 0 s2]
    @test AffineTransforms.tformrotate([0,0,1],pi/4).scalefwd ≈ [s2 -s2 0; s2 s2 0; 0 0 1]

    theta1 = pi/2
    @testset "> 2d" begin
        tfm1 = AffineTransforms.tformrotate(theta1)
        tfm2 = AffineTransforms.tformrotate(theta1)
        x = rand(2)
        @test tfm1*(tfm2*x) ≈ -x
        @test (tfm1*tfm2)*x ≈ -x
        @test tfm1\(tfm2*x) ≈ x
        @test (tfm1\tfm2)*x ≈ x
        @test (tfm1*tfm1*tfm1*tfm1)*x ≈ x
        r = AffineTransforms.tformfwd(tfm1, x...)
        @test [r...] == tfm1*x
        r = AffineTransforms.tforminv(tfm1, x...)
        @test [r...] == tfm1\x
    end

    @testset "> 3d" begin
        uaxis = rand(3); uaxis = uaxis/norm(uaxis)
        tfm1 = AffineTransforms.tformrotate(uaxis*theta1)
        tfm2 = AffineTransforms.tformrotate(uaxis*theta1)
        x = rand(3)
        xparallel = dot(x,uaxis)*uaxis
        xperp = x - xparallel
        @test tfm1*(tfm2*x) ≈ xparallel-xperp
        @test (tfm1*tfm2)*x ≈ xparallel-xperp
        @test tfm1\(tfm2*x) ≈ x
        @test (tfm1\tfm2)*x ≈ x
        @test (tfm1*tfm1*tfm1*tfm1)*x ≈ x
        r = AffineTransforms.tformfwd(tfm1, x...)
        @test [r...] == tfm1*x
        r = AffineTransforms.tforminv(tfm1, x...)
        @test [r...] == tfm1\x
    end
end

@testset "Random $(nd)D transformations" for nd = 1:4
    A = eye(nd) + 0.1*randn(nd,nd)
    b = randn(nd)
    tfm = AffineTransforms.AffineTransform(A, b)
    x = rand(nd)
    @test tfm*x ≈ A*x+b
    @test tfm\x ≈ A\(x-b)
    @test (tfm\tfm)*x ≈ x
    # Scaling transformations
    s = rand(nd)
    y = tfm*x
    @test AffineTransforms.scale(s, tfm)*x ≈ s.*y
    @test AffineTransforms.scale(tfm, s)*x ≈ tfm*(s.*x)
    if 2 <= nd <= 3
        r = AffineTransforms.tformfwd(tfm, x...)
        @test [r...] == tfm*x
        r = AffineTransforms.tforminv(tfm, x...)
        @test [r...] == tfm\x
    end
end

@testset "Pure rotations in 2d and 3d" begin
    for i = 1:20
        th = pi*(rand()-0.5)
        R = AffineTransforms.rotation2(th)
        @test th ≈ AffineTransforms.rotationparameters(R)[1]
    end
    for i = 1:100
        ax = randn(3)
        n = norm(ax)
        if n == 0
            continue
        end
        ax = ax / n
        th = pi*rand()
        rv = th*ax
        R = AffineTransforms.rotation3(rv)
        @test rv ≈ AffineTransforms.rotationparameters(R)
    end
end

@testset "Transformation of arrays" begin
    using Interpolations, Images # Images for padarray
    import AffineTransforms: center

    padwith0(A) = parent(padarray(A, Fill(0, (2,2))))

    for (IT,GT) in ((BSpline(Constant()), OnCell()),
                    (BSpline(Linear()), OnGrid()),
                    (BSpline(Quadratic(Flat())), OnCell()))
        A = padwith0([1 1; 1 2])
        itp = interpolate(A, IT, GT)
        tfm = AffineTransforms.tformtranslate([1,0])
        tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
        @test @inferred(AffineTransforms.transform(tA))[3:4,3:4] ≈ [1 2; 0 0]
        @test tA[3,3] ≈ 1
        @test tA[3,4] ≈ 2
        tfm = AffineTransforms.tformtranslate([0,1])
        tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
        @test AffineTransforms.transform(tA)[3:4,3:4] ≈ [1 0; 2 0]
        tfm = AffineTransforms.tformtranslate([-1,0])
        tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
        @test AffineTransforms.transform(tA)[3:4,3:4] ≈ [0 0; 1 1]
        tfm = AffineTransforms.tformtranslate([0,-1])
        tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
        @test AffineTransforms.transform(tA)[3:4,3:4] ≈ [0 1; 0 1]

        A = padwith0([2 1; 1 1])
        itp = interpolate(A, IT, GT)
        tfm = AffineTransforms.tformrotate(pi/2)
        tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
        @test AffineTransforms.transform(tA)[3:4,3:4] ≈ [1 2; 1 1]
        tfm = AffineTransforms.tformrotate(-pi/2)
        tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
        @test AffineTransforms.transform(tA)[3:4,3:4] ≈ [1 1; 2 1]

        # Check that getindex and transform yield the same results
        A = rand(6,6)
        itp = interpolate(A, IT, GT)
        tfm = AffineTransforms.tformrotate(pi/7)*AffineTransforms.tformtranslate(rand(2))
        tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
        dest = AffineTransforms.transform(tA)
        tfm_recentered = AffineTransforms.AffineTransform(tfm.scalefwd, tfm.offset + center(A) - tfm.scalefwd*center(dest))
        tA_recentered = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm_recentered)
        nbad = 0
        for j = 1:size(dest,2), i = 1:size(dest,1)
            val = tA_recentered[i,j]
            # try
            #     @test isfinite(dest[i,j]) == isfinite(val)
            # catch err
            #     @show dest tA_recentered
            #     rethrow(err)
            # end
            nbad += isfinite(dest[i,j]) != isfinite(val)
            if isfinite(val) && isfinite(dest[i,j])
                @test abs(val-dest[i,j]) < 1e-12
            end
        end
        @test nbad < 3
        dest = AffineTransforms.transform(tA, origin_src=zeros(2), origin_dest=zeros(2))
        nbad = 0
        for j = 1:size(dest,2), i = 1:size(dest,1)
            val = tA[i,j]
            # try
            #     @test isfinite(dest[i,j]) == isfinite(val)
            # catch err
            #     @show dest tA
            #     rethrow(err)
            # end
            nbad += isfinite(dest[i,j]) != isfinite(val)
            if isfinite(val) && isfinite(dest[i,j])
                @test abs(val-dest[i,j]) < 1e-12
            end
        end
        @test nbad < 3
    end

    A = Float64[1 2 3 4; 5 6 7 8]
    tfm = AffineTransforms.tformeye(2)
    dest = zeros(2,2)

    itp = interpolate(A, BSpline(Linear()), OnGrid())
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    @test AffineTransforms.transform!(dest, tA) ≈ [2 3; 6 7]
    dest3 = zeros(3,3)
    @test all(Compat.isapprox.(AffineTransforms.transform!(dest3, tA), [NaN NaN NaN; 3.5 4.5 5.5; NaN NaN NaN], nans=true))

    itp = interpolate(A, BSpline(Constant()), OnCell())
    tA = AffineTransforms.TransformedArray(extrapolate(itp, NaN), tfm)
    @test AffineTransforms.transform!(dest, tA) == [2 3; 6 7]

    # Transforms with non-real numbers
    using DualNumbers
    a = AffineTransforms.tformtranslate([0,dual(1,0)])
    A = reshape(1:9, 3, 3)
    At = AffineTransforms.transform(A, a)
    @test At[:,1:2] == A[:,2:3]
end

using Unitful: @u_str
@testset "With separates types for linear transform and translation" begin
    aff = AffineTransforms.AffineTransform([0 -1; 1 0], [2, 3]u"nm")
    @test aff * ([1, 2]u"nm") == [0, 4]u"nm"
end

nothing


text:
@testset "constructors" begin
    # inner constructors
    let
        dp = reinterpret(Core.LLVMPtr{Int,AS.Generic}, C_NULL)
        CuDeviceArray{Int,1,AS.Generic}((1,), dp)
    end

    # outer constructors
    for I in [Int32,Int64]
        a = I(1)
        b = I(2)

        dp = reinterpret(CUDA.LLVMPtr{I,AS.Generic}, C_NULL)

        # not parameterized
        CuDeviceArray(b, dp)
        CuDeviceArray((b,), dp)
        CuDeviceArray((b,a), dp)

        # partially parameterized
        CuDeviceArray{I}(b, dp)
        CuDeviceArray{I}((b,), dp)
        CuDeviceArray{I}((a,b), dp)
        CuDeviceArray{I,1}(b, dp)
        CuDeviceArray{I,1}((b,), dp)
        @test_throws MethodError CuDeviceArray{I,1}((a,b), dp)
        @test_throws MethodError CuDeviceArray{I,2}(b, dp)
        @test_throws MethodError CuDeviceArray{I,2}((b,), dp)
        CuDeviceArray{I,2}((a,b), dp)

        # fully parameterized
        CuDeviceArray{I,1,AS.Generic}(b, dp)
        CuDeviceArray{I,1,AS.Generic}((b,), dp)
        @test_throws MethodError CuDeviceArray{I,1,AS.Generic}((a,b), dp)
        @test_throws MethodError CuDeviceArray{I,1,AS.Shared}((a,b), dp)
        @test_throws MethodError CuDeviceArray{I,2,AS.Generic}(b, dp)
        @test_throws MethodError CuDeviceArray{I,2,AS.Generic}((b,), dp)
        CuDeviceArray{I,2,AS.Generic}((a,b), dp)

        # type aliases
        CuDeviceVector{I}(b, dp)
        CuDeviceMatrix{I}((a,b), dp)
    end
end

@testset "basics" begin     # argument passing, get and setindex, length
    dims = (16, 16)
    len = prod(dims)

    function kernel(input::CuDeviceArray{Float32}, output::CuDeviceArray{Float32})
        i = (blockIdx().x-1i32) * blockDim().x + threadIdx().x

        if i <= length(input)
            output[i] = Float64(input[i])   # force conversion upon setindex!
        end

        return
    end

    input = round.(rand(Float32, dims) * 100)

    input_dev = CuArray(input)
    output_dev = CuArray(input)

    @cuda threads=len kernel(input_dev, output_dev)
    output = Array(output_dev)
    @test input ≈ output
end

@testset "iteration" begin     # argument passing, get and setindex, length
    dims = (16, 16)
    function kernel(input::CuDeviceArray{T}, output::CuDeviceArray{T}) where {T}
        acc = zero(T)
        for elem in input
            acc += elem
        end
        output[1] = acc
        return
    end

    input = round.(rand(Float32, dims) * 100)

    input_dev = CuArray(input)
    output_dev = CuArray(Float32[0])

    @cuda kernel(input_dev, output_dev)
    output = Array(output_dev)
    @test sum(input) ≈ output[1]
end

@testset "bounds checking" begin
    @testset "#313" begin
        function kernel(dest)
            dest[1] = 1
            nothing
        end
        tt = Tuple{SubArray{Float64,2,CuDeviceArray{Float64,2,AS.Global},
                            Tuple{UnitRange{Int64},UnitRange{Int64}},false}}

        ir = sprint(io->CUDA.code_llvm(io, kernel, tt))
        @test !occursin("jl_invoke", ir)
        CUDA.code_ptx(devnull, kernel, tt)
    end
end

@testset "views" begin
    function kernel(array)
        i = (blockIdx().x-1i32) * blockDim().x + threadIdx().x

        _sub = view(array, 2:length(array)-1)
        if i <= length(_sub)
            _sub[i] = i
        end

        return
    end

    array = zeros(Int64, 100)
    array_dev = CuArray(array)

    sub = view(array, 2:length(array)-1)
    for i in 1:length(sub)
        sub[i] = i
    end

    @cuda threads=100 kernel(array_dev)
    @test array == Array(array_dev)
end

@testset "non-Int index to unsafe_load" begin
    function load_index(a)
        return a[UInt64(1)]
    end

    a = [1]
    p = pointer(a)
    dp = reinterpret(Core.LLVMPtr{eltype(p), AS.Generic}, p)
    da = CUDA.CuDeviceArray(1, dp)
    load_index(da)
end

function kernel_shmem_reinterpet_equal_size!(y)
  a = CuDynamicSharedArray(Float32, (blockDim().x,))
  b = reinterpret(UInt32, a)
  a[threadIdx().x] = threadIdx().x
  b[threadIdx().x] += 1
  y[threadIdx().x] = a[threadIdx().x]
  return
end

function shmem_reinterpet_equal_size()
  threads = 4
  y = CUDA.zeros(threads)
  shmem = sizeof(Float32) * threads
  @cuda(
    threads = threads,
    blocks = 1,
    shmem = shmem,
    kernel_shmem_reinterpet_equal_size!(y)
  )
  return y
end

@testset "reinterpret shmem: equal size" begin
  gpu = shmem_reinterpet_equal_size()
  a = zeros(Float32, length(gpu))
  b = reinterpret(UInt32, a)
  a .= 1:length(b)
  b .+= 1
  @test collect(gpu) == a
end

function kernel_shmem_reinterpet_smaller_size!(y)
  a = CuDynamicSharedArray(UInt128, (blockDim().x,))
  i = threadIdx().x
  p = i + i * im
  q = i - i * im
  b = reinterpret(typeof(p), a)
  b[1 + 2 * (threadIdx().x - 1i32)] = p
  b[2 + 2 * (threadIdx().x - 1i32)] = q
  y[threadIdx().x] = a[threadIdx().x]
  return
end

function shmem_reinterpet_smaller_size()
  threads = 4
  y = CUDA.zeros(UInt128, threads)
  shmem = sizeof(UInt128) * threads
  @cuda(
    threads = threads,
    blocks = 1,
    shmem = shmem,
    kernel_shmem_reinterpet_smaller_size!(y)
  )
  return y
end

@testset "reinterpret shmem: smaller size" begin
  gpu = shmem_reinterpet_smaller_size()
  n = length(gpu)
  a = zeros(UInt128, n)
  p(i) = Int32(i) + Int32(i) * im
  q(i) = Int32(i) - Int32(i) * im
  b = reinterpret(typeof(p(0)), a)
  b[1:2:end] .= p.(1:n)
  b[2:2:end] .= q.(1:n)
  @test collect(gpu) == a
end

function kernel_shmem_reinterpet_larger_size!(y)
  a = CuDynamicSharedArray(Float32, (4 * blockDim().x,))
  b = reinterpret(UInt128, a)
  a[1 + 4 * (threadIdx().x - 1i32)] = threadIdx().x
  a[2 + 4 * (threadIdx().x - 1i32)] = threadIdx().x * 2
  a[3 + 4 * (threadIdx().x - 1i32)] = threadIdx().x * 3
  a[4 + 4 * (threadIdx().x - 1i32)] = threadIdx().x * 4
  y[threadIdx().x] = b[threadIdx().x]
  return
end

function shmem_reinterpet_larger_size()
  threads = 4
  y = CUDA.zeros(UInt128, threads)
  shmem = sizeof(UInt128) * threads
  @cuda(
    threads = threads,
    blocks = 1,
    shmem = shmem,
    kernel_shmem_reinterpet_larger_size!(y)
  )
  return y
end

@testset "reinterpret shmem: larger size" begin
  gpu = shmem_reinterpet_larger_size()
  n = length(gpu)
  b = zeros(UInt128, n)
  a = reinterpret(Float32, b)
  a[1:4:end] .= 1:n
  a[2:4:end] .= (1:n) .* 2
  a[3:4:end] .= (1:n) .* 3
  a[4:4:end] .= (1:n) .* 4
  @test collect(gpu) == b
end