text:
"""
This function returns the product of marginal cost and quantity provided.
"""
function calculate_cost(marginal_cost::Float64, Q::JuMP.VariableRef)
    return marginal_cost * Q
end

"""
This function returns the product of market price and quantity provided.
"""
function calculate_revenue(price::Float64, Q::JuMP.VariableRef)
    return price * Q
end

"""
This function returns the product of marginal cost and quantity provided.
"""
function calculate_cost(marginal_cost::Float64, Q::JuMP.GenericAffExpr{Float64, JuMP.VariableRef})
    return marginal_cost * Q
end

"""
This function returns the product of market price and quantity provided.
"""
function calculate_revenue(price::Float64, Q::JuMP.GenericAffExpr{Float64, JuMP.VariableRef})
    return price * Q
end

"""
This function does nothing if the product is not either energy, reserve up or reserve down.
"""
function add_to_generationlimits!(min_gen::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                 max_gen::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                 Q::JuMP.VariableRef,
                                 product::T) where T <: Product
return
end

"""
This function adds to the maximum generation expression if the
product is Reserve Up.
Returns nothing.
"""
function add_to_generationlimits!(min_gen::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                 max_gen::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                 Q::JuMP.VariableRef,
                                 product::OperatingReserve{ReserveUpEMIS})
    JuMP.add_to_expression!(max_gen, Q)
    return
end

"""
This function adds to the minimum generation expression if the
product is Reserve Down.
Returns nothing.
"""
function add_to_generationlimits!(min_gen::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                 max_gen::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                 Q::JuMP.VariableRef,
                                 product::OperatingReserve{ReserveDownEMIS})
    JuMP.add_to_expression!(min_gen, -Q)
    return
end

"""
This function adds to both minimum and maximum generation expressions if the
product is Energy.
Returns nothing.
"""
function add_to_generationlimits!(min_gen::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                 max_gen::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                 Q::JuMP.VariableRef,
                                 product::Energy)
    JuMP.add_to_expression!(max_gen, Q)
    JuMP.add_to_expression!(min_gen, Q)
    return
end

"""
This function calculates operating market profits for a generator
based on expected market prices stored in the investor struct.
Returns an axis array of profits for each operating market product.
"""
function calculate_operating_profit(gen::P,
                                    scenario_name::String,
                                    all_prices::MarketPrices,
                                    start_year::Int64,
                                    end_year::Int64,
                                    rep_hour_weight::Vector{Float64},
                                    solver::JuMP.MOI.OptimizerWithAttributes) where P <: GeneratorEMIS{<: BuildPhase}

    num_hours = length(rep_hour_weight)

    products = find_operating_products(get_products(gen))

    market_prices = Dict(get_name(prod) => get_prices(all_prices, prod)[scenario_name][get_zone(get_tech(gen)), :, :] for prod in products)

    max_perc = Dict(zip(get_name.(products), get_maxperc.(products, scenario_name, end_year, num_hours)))

    m = JuMP.Model(solver)

    JuMP.@variable(m, Q[i in get_name.(products),
                   y in start_year:end_year,
                   h in 1:num_hours] >= 0)

    JuMP.@variable(m, temp == 0)

    JuMP.@expression(m, profit[i in 1:length(products), y in start_year:end_year], 0 + temp)

    JuMP.@expression(m, min_gen[y in start_year:end_year, h in 1:num_hours], 0 + temp)

    JuMP.@expression(m, max_gen[y in start_year:end_year, h in 1:num_hours], 0 + temp)

    for y in start_year:end_year
        for h in 1:num_hours
            for (i, product) in enumerate(products)

                hourly_profit = rep_hour_weight[h] *
                         (calculate_revenue(market_prices[get_name(product)][y, h],
                                       Q[get_name(product), y, h]) -
                          calculate_cost(get_marginal_cost(product),
                                   Q[get_name(product), y, h]))

                JuMP.add_to_expression!(profit[i, y], hourly_profit)

                add_to_generationlimits!(min_gen[y, h],
                                        max_gen[y, h],
                                        Q[get_name(product), y, h],
                                        product)
                JuMP.@constraint(m, Q[get_name(product), y, h] <=
                                max_perc[get_name(product)][y, h] * get_maxcap(gen))
            end
            JuMP.@constraint(m, min_gen[y,h] >= get_mincap(gen) * 0)
            JuMP.@constraint(m, max_gen[y,h] <= get_maxcap(gen))
        end
    end

    JuMP.@objective(m, Max, sum(profit))

    JuMP.optimize!(m)

    profit_array = AxisArrays.AxisArray(value.(profit),
                   AxisArrays.Axis{:prod}(get_name.(products)),
                   AxisArrays.Axis{:year}(start_year:end_year))

    energy_production = AxisArrays.AxisArray([sum(value(Q[:Energy, y, h]) * rep_hour_weight[h] for h in 1:num_hours) for y in start_year:end_year],
                        AxisArrays.Axis{:year}(start_year:end_year))

    return profit_array, energy_production

end

"""
This function does nothing if the product is not of Energy, ReserveUp or ReserveDown type
"""
function storage_product_constraints(product::T,
                                    Q::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                    p_out::JuMP.VariableRef,
                                    p_in::JuMP.VariableRef,
                                    p_ru::JuMP.VariableRef,
                                    p_rd::JuMP.VariableRef
                                    ) where T <: OperatingProduct
    return
end

"""
This function calculates the net energy provided by the storage device.
"""
function storage_product_constraints(product::Energy,
                                    Q::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                    p_out::JuMP.VariableRef,
                                    p_in::JuMP.VariableRef,
                                    p_ru::JuMP.VariableRef,
                                    p_rd::JuMP.VariableRef
                                    )
    JuMP.add_to_expression!(Q, p_out - p_in)
    return
end

"""
This function calculates the net reserve up provided by the storage device.
"""
function storage_product_constraints(product::ReserveUpEMIS,
                                    Q::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                    p_out::JuMP.VariableRef,
                                    p_in::JuMP.VariableRef,
                                    p_ru::JuMP.VariableRef,
                                    p_rd::JuMP.VariableRef
                                    )
    JuMP.add_to_expression!(Q, p_ru)
    return
end

"""
This function calculates the net reserve down provided by the storage device.
"""
function storage_product_constraints(product::ReserveDownEMIS,
                                    Q::JuMP.GenericAffExpr{Float64, JuMP.VariableRef},
                                    p_out::JuMP.VariableRef,
                                    p_in::JuMP.VariableRef,
                                    p_ru::JuMP.VariableRef,
                                    p_rd::JuMP.VariableRef
                                    )
    JuMP.add_to_expression!(Q, p_rd)
    return
end

"""
This function calculates operating market profits for a storage device
based on expected market prices stored in the investor struct.
Returns an axis array of profits for each operating market product.
"""
function calculate_operating_profit(storage::P,
                                    scenario_name::String,
                                    all_prices::MarketPrices,
                                    start_year::Int64,
                                    end_year::Int64,
                                    rep_hour_weight::Vector{Float64},
                                    solver::JuMP.MOI.OptimizerWithAttributes ) where P <: StorageEMIS{<: BuildPhase}

    num_hours = length(rep_hour_weight)
    tech = get_tech(storage)

    products = find_operating_products(get_products(storage))
    energy_product = find_energy_product(products)

    capacity_factors = get_capacity_factors(energy_product[1])[scenario_name]

    market_prices = Dict(get_name(prod) => get_prices(all_prices, prod)[scenario_name][get_zone(tech), :, :] for prod in products)

    max_gen = get_maxcap(storage)                                   # maximum generation limit
    min_gen = get_mincap(storage)                                    # minimum generation limit

    input_power_limits = get_input_active_power_limits(tech)
    min_input = input_power_limits[:min]                            # minimum charging limit
    max_input = input_power_limits[:max]                            # maximum charging limit

    max_reserve_up = get_project_max_reserveup(storage)              # maximum reserve up limit
    max_reserve_down = get_project_max_reservedown(storage)           # maximum reserve down limit

    storage_capacity = get_storage_capacity(tech)
    min_storage = storage_capacity[:min]                            # minimum storage capacity
    max_storage = storage_capacity[:max]                            # maximum storage capacity

    efficiency = get_efficiency(tech)
    efficiency_in = efficiency[:in]                                 # charging efficiency
    efficiency_out = efficiency[:out]                               # discharging efficiency

    initstorage = get_soc(tech)
    m = JuMP.Model(solver)

    JuMP.@variable(m, p_e[p in start_year:end_year, t in 1:num_hours] >= 0) # Unit energy production [MW]
    JuMP.@variable(m, p_in[p in start_year:end_year, t in 1:num_hours] >= 0) # Storage charging [MW]

    JuMP.@variable(m, p_ru[p in start_year:end_year, t in 1:num_hours] >= 0) # Unit energy production [MW]
    JuMP.@variable(m, p_rd[p in start_year:end_year, t in 1:num_hours] >= 0) # Storage charging [MW]

    JuMP.@variable(m, p_in_ru[p in start_year:end_year, t in 1:num_hours] >= 0) # reserve up provided by storage charging [MW]
    JuMP.@variable(m, p_in_rd[p in start_year:end_year, t in 1:num_hours] >= 0) # reserve down provided by storage charging [MW]

    JuMP.@variable(m, p_out_ru[p in start_year:end_year, t in 1:num_hours] >= 0) # reserve up provided by storage discharging [MW]
    JuMP.@variable(m, p_out_rd[p in start_year:end_year, t in 1:num_hours] >= 0) # reserve down provided by storage discharging [MW]

    JuMP.@variable(m, temp == 0)

    # Quantity of products provided
    JuMP.@expression(m, Q[i in get_name.(products), p in start_year:end_year, t in 1:num_hours], 0 + temp)

    # Profit for each product
    JuMP.@expression(m, profit[i in 1:length(products), p in start_year:end_year], 0 + temp)

    # Storage level evolution
    JuMP.@expression(m, storage_level[p in start_year:end_year, t in 1:num_hours],
                    p_in[p, t] * efficiency_in - p_e[p, t] / efficiency_out)

    for p in start_year:end_year
        for t in 1:num_hours

            for (i, product) in enumerate(products)

                storage_product_constraints(product,
                                            Q[get_name(product), p, t],
                                            p_e[p, t],
                                            p_in[p, t],
                                            p_ru[p, t],
                                            p_rd[p, t])

                hourly_profit = rep_hour_weight[t] *
                                     (calculate_revenue(market_prices[get_name(product)][p, t],
                                                   Q[get_name(product), p, t]) -
                                      calculate_cost(get_marginal_cost(product),
                                               Q[get_name(product), p, t]))

                            JuMP.add_to_expression!(profit[i, p], hourly_profit)

            end

            if t == 1
                JuMP.add_to_expression!(storage_level[p, t], initstorage)
            else
                JuMP.add_to_expression!(storage_level[p, t], storage_level[p, t - 1])
            end

            # Storage technical constraints
            JuMP.@constraint(m, p_e[p, t] <= capacity_factors[p, t] * max_gen) # Capacity factor constraint

            JuMP.@constraint(m, p_e[p,t] - p_out_rd[p,t] >= min_gen * 0) # Minimum output dispatch
            JuMP.@constraint(m, p_e[p,t] + p_out_ru[p,t] <= max_gen) # Maximum output dispatch

            JuMP.@constraint(m, p_in[p,t] - p_in_ru[p,t] >=  min_input * 0) # Minimum input dispatch
            JuMP.@constraint(m, p_in[p,t] + p_in_rd[p,t] <= max_input) # Maximum input dispatch

            JuMP.@constraint(m, p_ru[p, t] <= p_out_ru[p, t] + p_in_ru[p, t]) # Total storage reserve up
            JuMP.@constraint(m, p_rd[p, t] <= p_out_rd[p, t] + p_in_rd[p, t]) # Total storage reserve down

            JuMP.@constraint(m, storage_level[p, t] >= min_storage) # Minimum storage level
            JuMP.@constraint(m, storage_level[p, t] <= max_storage) # Maximum storage level

        end
    end

    JuMP.@objective(m, Max, sum(profit))

    JuMP.optimize!(m)

    profit_array = AxisArrays.AxisArray(value.(profit),
                   AxisArrays.Axis{:prod}(get_name.(products)),
                   AxisArrays.Axis{:year}(start_year:end_year))

    energy_production = AxisArrays.AxisArray([sum(value(Q[:Energy, p, t]) * rep_hour_weight[t] for t in 1:num_hours) for p in start_year:end_year],
                        AxisArrays.Axis{:year}(start_year:end_year))

    return profit_array, energy_production

end

"""
This function updates the expected profit array of the project
if the product is an operating market product.
Returns nothing.
"""
function update_profit!(project::P,
                        scenario_name::String,
                        product::T,
                        operating_profit::AxisArrays.AxisArray{Float64, 2},
                        capacity_revenues::Vector{Float64},
                        capacity_forward_years::Int64,
                        rec_revenues::Vector{Float64},
                        price_years::NamedTuple{(:start_year, :end_year),
                                     Tuple{Int64, Int64}},
                        update_years::NamedTuple{(:start_year, :end_year),
                                      Tuple{Int64, Int64}},
                        iteration_year::Int64) where {P <: Project{<: BuildPhase}, T <: OperatingProduct}

    profit = [operating_profit[get_name(product), y]
            for y in price_years[:start_year]:price_years[:end_year]]

    set_scenario_profit!(get_finance_data(project),
                scenario_name,
                get_name(product),
                iteration_year,
                update_years[:start_year],
                update_years[:end_year],
                profit)

    return
end


text:
"""
    FlipSignBridge{T, S1, S2}

Bridge constrained variables in `S1` into constrained variables in `S2` by
multiplying the variables by `-1` and taking the point reflection of the set
across the origin. The flipped `MOI.VectorOfVariables`-in-`S` constraint is
stored in the `flipped_constraint` field by convention.
"""
abstract type FlipSignBridge{T,S1<:MOI.AbstractSet,S2<:MOI.AbstractSet} <:
              AbstractBridge end

function supports_constrained_variable(
    ::Type{<:FlipSignBridge{T,S1}},
    ::Type{S1},
) where {T,S1<:MOI.AbstractVectorSet}
    return true
end
function MOIB.added_constrained_variable_types(
    ::Type{<:FlipSignBridge{T,S1,S2}},
) where {T,S1,S2}
    return [(S2,)]
end
function MOIB.added_constraint_types(::Type{<:FlipSignBridge})
    return Tuple{DataType,DataType}[]
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::FlipSignBridge, ::MOI.NumberOfVariables)
    return length(bridge.flipped_variables)
end
function MOI.get(bridge::FlipSignBridge, ::MOI.ListOfVariableIndices)
    return bridge.flipped_variables
end
function MOI.get(
    ::FlipSignBridge{T,S1,S2},
    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,S2},
) where {T,S1,S2<:MOI.AbstractVectorSet}
    return 1
end
function MOI.get(
    bridge::FlipSignBridge{T,S1,S2},
    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S2},
) where {T,S1,S2<:MOI.AbstractVectorSet}
    return [bridge.flipped_constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::FlipSignBridge)
    return MOI.delete(model, bridge.flipped_variables)
end

function MOI.delete(
    model::MOI.ModelLike,
    bridge::FlipSignBridge,
    i::IndexInVector,
)
    MOI.delete(model, bridge.flipped_variables[i.value])
    return deleteat!(bridge.flipped_variables, i.value)
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::FlipSignBridge{T,S1},
) where {T,S1<:MOI.AbstractVectorSet}
    return S1(length(bridge.flipped_variables))
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintDual},
    bridge::FlipSignBridge,
)
    return -MOI.get(model, attr, bridge.flipped_constraint)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.VariablePrimal,MOI.VariablePrimalStart},
    bridge::FlipSignBridge,
    i::IndexInVector,
)
    return -MOI.get(model, attr, bridge.flipped_variables[i.value])
end

function MOIB.bridged_function(
    bridge::FlipSignBridge{T},
    i::IndexInVector,
) where {T}
    func = MOI.SingleVariable(bridge.flipped_variables[i.value])
    return MOIU.operate(-, T, func)
end
function unbridged_map(
    bridge::FlipSignBridge{T},
    vi::MOI.VariableIndex,
    i::IndexInVector,
) where {T}
    func = MOIU.operate(-, T, MOI.SingleVariable(vi))
    return (bridge.flipped_variables[i.value] => func,)
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimalStart,
    ::Type{<:FlipSignBridge},
)
    return MOI.supports(model, attr, MOI.VariableIndex)
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimalStart,
    bridge::FlipSignBridge,
    value,
    i::IndexInVector,
)
    return MOI.set(model, attr, bridge.flipped_variables[i.value], -value)
end

"""
    NonposToNonnegBridge{T, F<:MOI.AbstractVectorFunction, G<:MOI.AbstractVectorFunction} <:
        FlipSignBridge{T, MOI.Nonpositives, MOI.Nonnegatives, F, G}

Transforms constrained variables in `Nonpositives` into constrained variables in
`Nonnegatives`.
"""
struct NonposToNonnegBridge{T} <:
       FlipSignBridge{T,MOI.Nonpositives,MOI.Nonnegatives}
    flipped_variables::Vector{MOI.VariableIndex}
    flipped_constraint::MOI.ConstraintIndex{
        MOI.VectorOfVariables,
        MOI.Nonnegatives,
    }
end
function bridge_constrained_variable(
    ::Type{NonposToNonnegBridge{T}},
    model::MOI.ModelLike,
    set::MOI.Nonpositives,
) where {T}
    flipped_variables, flipped_constraint =
        MOI.add_constrained_variables(model, MOI.Nonnegatives(set.dimension))
    return NonposToNonnegBridge{T}(flipped_variables, flipped_constraint)
end

ext: jl


text:
# Object to store representation of prior
"""
    ConvexPriorClass

Abstract type representing convex classes of probability distributions ``\\mathcal{G}``.
"""
abstract type ConvexPriorClass end

"""
	GaussianMixturePriorClass(σ_prior, grid, solver)

Class of distributions ``G \\in \\mathcal{G}`` that can be written as the convolution of a Normal
distribution with standard deviation `σ_prior` and a discrete distribution ``\\pi``
supported on `grid`.

```math
G = \\mathcal{N}(0,\\sigma_{\\text{prior}}^2) \\star \\pi
```

The `solver` object is a `MathOptInterface` compatible optimizer such as `Mosek.Optimizer` that will be
used to solve the convex modulus of continuity (or worst-case-bias estimation) problem.
"""
mutable struct GaussianMixturePriorClass{T<:Real,
                                         VT<:AbstractVector{T}} <: ConvexPriorClass
    σ_prior::Float64
    grid::VT
    solver
    solver_params
end

GaussianMixturePriorClass(σ_prior, grid, solver) = GaussianMixturePriorClass(σ_prior, grid, solver, [])
#GaussianMixturePriorClass(σ_prior, grid) = GaussianMixturePriorClass(σ_prior, grid, Gurobi.Optimizer)

length(gmix_class::GaussianMixturePriorClass) = length(gmix_class.grid)
location(gmix_class::GaussianMixturePriorClass) = gmix_class.grid

function add_prior_variables!(model, gmix_class::GaussianMixturePriorClass; var_name = "πs") # adds constraints
    n_priors = length(gmix_class)
    tmp_vars = @variable(model, [i=1:n_priors])
    model[Symbol(var_name)] = tmp_vars
    set_lower_bound.(tmp_vars, 0.0)
    con = @constraint(model, sum(tmp_vars) == 1.0)
    tmp_vars
end

function marginalize(gmix_class::GaussianMixturePriorClass,
                     z::DiscretizedStandardNormalSamples,
                     param_vec) #-> A*param_vec
    σ_prior = gmix_class.σ_prior
    grid = gmix_class.grid
    A_mat = hcat([ pdf(marginalize(Normal(μ, σ_prior), z)) for μ in grid]...)
    A_mat*param_vec
end

function linear_functional(gmix_class::GaussianMixturePriorClass,
                           target::EBayesTarget,
                           param_vec)
    σ_prior = gmix_class.σ_prior
    grid = gmix_class.grid
    v = [target(Normal(μ, σ_prior)) for μ in grid]
    dot(v, param_vec)
end

function (gmix_class::GaussianMixturePriorClass)(param_vec)
    σ_prior = gmix_class.σ_prior
    grid = gmix_class.grid
    #TODO: check param_vec is probability vector
    MixtureModel(Normal, [(μ, σ_prior) for μ in grid], param_vec)
end

# TODO: more rigorous projection onto simplex?
function fix_πs(gmix_class::GaussianMixturePriorClass, πs)
    πs = max.(πs, 0.0)
    πs = πs ./ sum(πs)
end

fix_πs(::ConvexPriorClass, πs) = πs

# In fact all computations with marginalized object.

# param -> A*x

# finite-dim representation, here: probs
# want to map: probs ->

function worst_case_bias(general_Q,
                  Z::DiscretizedStandardNormalSamples,
                  gmix_class::ConvexPriorClass,
                  target::LinearEBayesTarget;
                  boundary_mass = Inf)
     grid = Z.mhist.grid
     tmpcalib = DiscretizedAffineEstimator(MCEBHistogram(grid), general_Q)
     worst_case_bias(tmpcalib, Z, gmix_class, target; boundary_mass=boundary_mass)
end

"""
    worst_case_bias(Q,
                    Zs_discr::DiscretizedStandardNormalSamples,
                    prior_class::ConvexPriorClass,
                    target::LinearEBayesTarget)

Computes the worst-case positive and negative bias of the linear estimator `Q` for
estimating the linear functional `target` over `prior_class` based on
[`DiscretizedStandardNormalSamples`](@ref) `Zs_discr`.
"""
function worst_case_bias(Q::DiscretizedAffineEstimator,
                  Z::DiscretizedStandardNormalSamples,
                  prior_class::ConvexPriorClass,
                  target::LinearEBayesTarget;
                  boundary_mass = Inf)

    model = Model(optimizer_with_attributes(prior_class.solver, prior_class.solver_params...))
    πs = add_prior_variables!(model, prior_class; var_name = "πs")
    fs = marginalize(prior_class, Z, πs)
    L = linear_functional(prior_class, target, πs)

    if (!isnothing(Z.f_max))
        idx_upper = idx_enforce_upper_bound(prior_class, Z)#
        @constraint(model, f1_upper, fs[idx_upper] .<= Z.f_max[idx_upper])
    end

    if (!isnothing(Z.f_min))
        idx_nonzero = idx_enforce_lower_bound(prior_class, Z) #
        @constraint(model, f1_lower, fs[idx_nonzero] .>= Z.f_min[idx_nonzero])
    end

    if boundary_mass < Inf
        @constraint(model, fs[1] <= boundary_mass)
        @constraint(model, fs[end] <= boundary_mass)
    end
    #    @constraint(jm, f3 .- f_marginal .<= C*h_marginal_grid)
    #    @constraint(jm, f3 .- f_marginal .>= -C*h_marginal_grid)
    #end

    @objective(model, Max, Q.Qo + dot(Q.Q,fs) - L)
    optimize!(model)
    max_bias = objective_value(model)

    # TODO let prior_class function do this by itself
    max_g = prior_class(fix_πs(prior_class, JuMP.value.(πs)))

    @objective(model, Min, Q.Qo + dot(Q.Q,fs) - L)
    optimize!(model)
    min_bias = objective_value(model)
    min_g = prior_class(fix_πs(prior_class, JuMP.value.(πs)))

    max_abs_bias = max(abs(max_bias), abs(min_bias))

    (max_abs_bias = max_abs_bias,
     max_squared_bias = abs2(max_abs_bias),
     max_bias = max_bias, min_bias = min_bias,
     max_g = max_g, min_g = min_g,
     model = model)
end

# also define this for a single prior
function worst_case_bias(Q::DiscretizedAffineEstimator,
                  Z::DiscretizedStandardNormalSamples,
                  prior::Distribution,
                  target::LinearEBayesTarget)

    marginal_distr = marginalize(prior, Z)
    #TODO: define expectation operator for this type of thing?
    expected_target = Q.Qo + dot(pdf(marginal_distr), Q.Q)
    L = target(prior)

    bias = expected_target - L
    (max_abs_bias = abs(bias),
     max_squared_bias = abs2(bias),
     bias = bias)
end

#function initialize_modulus_problem( )

#nd

#function optimize!(modulus_problem, Z, gmix_class, target, δ)

#end

# helper functions for working with model.
function get_δ(model; recalculate_δ = false)
    if recalculate_δ
        norm(JuMP.value.(model[:Δf]))
    else
      JuMP.value(model[:δ_up])
    end
end

# TODO: figure out why this errors
function idx_enforce_lower_bound(prior_class::ConvexPriorClass, Z::DiscretizedStandardNormalSamples)
    1:length(Z.f_min)#1:length(Z.f_min)
end

function idx_enforce_lower_bound(prior_class::GaussianMixturePriorClass, Z::DiscretizedStandardNormalSamples)
    findall( Z.f_min .> 0.0)
end

# TODO: figure out why this errors
function idx_enforce_upper_bound(prior_class::ConvexPriorClass, Z::DiscretizedStandardNormalSamples)
    findall( Z.f_max .< 1/sqrt(2π))#1:length(Z.f_max)
end

function idx_enforce_upper_bound(prior_class::GaussianMixturePriorClass, Z::DiscretizedStandardNormalSamples)
    findall( Z.f_max .< 1/sqrt(2π))
end

function initialize_modulus_problem(Z::DiscretizedStandardNormalSamples,
                                    prior_class::ConvexPriorClass,
                                    target::EBayesTarget)

    model = Model(optimizer_with_attributes(prior_class.solver, prior_class.solver_params...))
    πs1 = add_prior_variables!(model, prior_class; var_name = "πs1")
    πs2 = add_prior_variables!(model, prior_class; var_name = "πs2")

    fs1 = marginalize(prior_class, Z, πs1)
    fs2 = marginalize(prior_class, Z, πs2)

    model[:fs1] = fs1
    model[:fs2] = fs2

    if (!isnothing(Z.f_max))
        idx_upper = idx_enforce_upper_bound(prior_class, Z)#
        @constraint(model, f1_upper, fs1[idx_upper] .<= Z.f_max[idx_upper])
        @constraint(model, f2_upper, fs2[idx_upper] .<= Z.f_max[idx_upper])
    end

    if (!isnothing(Z.f_min))
        idx_nonzero = idx_enforce_lower_bound(prior_class, Z) #
        @constraint(model, f1_lower, fs1[idx_nonzero] .>= Z.f_min[idx_nonzero])
        @constraint(model, f2_lower, fs2[idx_nonzero] .>= Z.f_min[idx_nonzero])
    end

    L1 = linear_functional(prior_class, target, πs1)
    L2 = linear_functional(prior_class, target, πs2)

    model[:L1] = L1
    model[:L2] = L2

    f̄s_sqrt = sqrt.(Z.var_proxy)
    model[:f̄s_sqrt] = f̄s_sqrt

    @variable(model, δ_up)
    model[:δ_up] = δ_up

    Δf = (fs1 .- fs2)./f̄s_sqrt
    model[:Δf] = Δf
    @constraint(model, pseudo_chisq_constraint,
               [δ_up; Δf] in SecondOrderCone())

    @objective(model, Max, L2 -L1)
    optimize!(model)
    δ = get_δ(model; recalculate_δ = true)
    model[:δ_max] = δ
    @constraint(model, bound_delta, δ_up == δ)
    model
end

function modulus_at_δ!(model, δ)
    set_normalized_rhs(model[:bound_delta], δ)
    optimize!(model)
    model
end

function δ_path(model, δs)
    max_biases = Vector{Float64}(undef, length(δs))
    unit_variances   = Vector{Float64}(undef, length(δs))
    for (i,δ) in enumerate(δs)
        model = modulus_at_δ!(model, δ)
        tmp_bias, tmp_var = get_bias_var(model)
        max_biases[i] = tmp_bias
        unit_variances[i] = tmp_var
    end
    max_biases, unit_variances
end

default_δ_min(n, C∞_density) = sqrt( log(n)/n)*C∞_density

#function solve_modulus_problem

"""
    DeltaTuner

Abstract type used to represent ways of picking
``\\delta`` at which to solve the modulus problem, cf.
Manuscript. Different choices of ``\\delta`` correspond
to different choices of the Bias-Variance tradeoff with
every choice leading to Pareto-optimal tradeoff.
"""
abstract type DeltaTuner end

"""
    FixedDelta(δ) <: DeltaTuner

The simplest `DeltaTuner`. It just uses `δ` which has been fixed
beforehand.  This is used as a component of more complicated
`DeltaTuner`s which need to run repeatedly for different values
of `δ`.
"""
struct FixedDelta <: DeltaTuner
    δ::Float64
end

function optimal_δ(model::JuMP.Model, objective::FixedDelta)
    objective.δ
end

#struct VarProxyUpperBound <: DeltaTuner
#
#end

abstract type BiasVarAggregate <:DeltaTuner end

function get_bias_var(model::JuMP.Model)
    δ = get_δ(model)
    ω_δ = objective_value(model)
    ω_δ_prime = -JuMP.dual(model[:bound_delta])
    max_bias = (ω_δ - δ*ω_δ_prime)/2
    unit_var_proxy = ω_δ_prime^2
    max_bias, unit_var_proxy
end

function (bv::BiasVarAggregate)(model::JuMP.Model)
    bv(get_bias_var(model)...)
end

"""
    RMSE(n::Integer, δ_min::Float64) <: DeltaTuner

A `DeltaTuner` that chooses the `δ ≧ δ_min` the optimizes
the worst-case (root) mean squared error.  Here `n` is the sample
size used for estimation.
"""
struct RMSE <: BiasVarAggregate
    n::Integer
    δ_min::Float64
end

(rmse::RMSE)(bias, unit_var_proxy) =  sqrt(bias^2 + unit_var_proxy/rmse.n)

"""
    HalfCIWidth(n::Integer, α::Float64 δ_min::Float64) <: DeltaTuner

A `DeltaTuner` that chooses the `δ ≧ δ_min` the optimizes
the worst-case (root) mean squared error.  Here `n` is the sample
size used for estimation.
"""
struct HalfCIWidth <: BiasVarAggregate
    n::Integer
    α::Float64
    δ_min::Float64
end

HalfCIWidth(n::Integer, δ_min::Float64) = HalfCIWidth(n, 0.9, δ_min)

function (half_ci::HalfCIWidth)(bias, unit_var_proxy)
    se = sqrt(unit_var_proxy/half_ci.n)
    bias_adjusted_gaussian_ci(se, maxbias=bias, level=half_ci.α)
end

function optimal_δ(model::JuMP.Model, objective::BiasVarAggregate)
    f = δ -> objective(modulus_at_δ!(model, δ))
    δ_max = model[:δ_max]
    δ_min = objective.δ_min
    Optim.optimize(f, δ_min, δ_max; rel_tol=1e-3).minimizer
end

function subfamily_mse(model, n)
    ω_δ = objective_value(model)
    ω_δ_prime = -JuMP.dual(model[:bound_delta])
    δ = JuMP.normalized_rhs(bound_delta)
    max_bias = (ω_δ - δ*ω_δ_prime)/2
    unit_var_proxy = ω_δ_prime^2

    -ω_δ^2/(4 + n*δ^2)
end

Base.@kwdef mutable struct SteinMinimaxEstimator
                Z::DiscretizedStandardNormalSamples
                prior_class::ConvexPriorClass
                target::EBayesTarget
                δ::Float64
                ω_δ::Float64
                ω_δ_prime::Float64
                g1
                g2
                Q::DiscretizedAffineEstimator
                max_bias::Float64
                unit_var_proxy::Float64
                n::Int
                δ_tuner::DeltaTuner
                δ_opt::Float64
                model
end

var(sme::SteinMinimaxEstimator) = sme.unit_var_proxy/sme.n

function worst_case_bias(sme::SteinMinimaxEstimator)
    sme.max_bias
end

"""
    SteinMinimaxEstimator(Zs_discr::DiscretizedStandardNormalSamples,
                          prior_class::ConvexPriorClass,
                          target::LinearEBayesTarget,
                          δ_tuner::DeltaTuner

Computes a linear estimator optimizing a worst-case bias-variance tradeoff (specified by `δ_tuner`)
for estimating a linear `target` over `prior_class` based on [`DiscretizedStandardNormalSamples`](@ref)
`Zs_discr`.
"""
function SteinMinimaxEstimator(Zs_discr::DiscretizedStandardNormalSamples,
                               gmix::ConvexPriorClass,
                               target::LinearEBayesTarget,
                               δ_tuner::DeltaTuner)

   modulus_problem = initialize_modulus_problem(Zs_discr, gmix, target)
   δ_opt = optimal_δ(modulus_problem, δ_tuner)
   modulus_problem = modulus_at_δ!(modulus_problem, δ_opt)
   SteinMinimaxEstimator(Zs_discr, gmix, target, modulus_problem;
                         δ = δ_opt, δ_tuner = δ_tuner)
end

function SteinMinimaxEstimator(Z::DiscretizedStandardNormalSamples,
                               prior_class::ConvexPriorClass,
                               target::LinearEBayesTarget,
                               model::JuMP.Model;
                               δ = get_δ(model),
                               δ_tuner = FixedDelta(δ)
                              )

    ω_δ = objective_value(model)
    ω_δ_prime = -JuMP.dual(model[:bound_delta])

    πs1 = JuMP.value.(model[:πs1])
    πs2 = JuMP.value.(model[:πs2])

    πs1 = fix_πs(prior_class, πs1)
    πs2 = fix_πs(prior_class, πs2)

    g1 = prior_class(πs1)
    g2 = prior_class(πs2)

    L1 = JuMP.value(model[:L1])
    L2 = JuMP.value(model[:L2])
    #construct estimator
    f1 = marginalize(g1, Z)
    f2 = marginalize(g2, Z)

    f̄s = Z.var_proxy
    Q = ω_δ_prime/δ*(pdf(f2) .- pdf(f1))./Z.var_proxy
    Q_0  = (L1+L2)/2 -
           ω_δ_prime/(2*δ)*sum( (pdf(f2) .- pdf(f1)).* (pdf(f2) .+ pdf(f1)) ./ f̄s)

    stein = DiscretizedAffineEstimator(Z, Q, Q_0)

    max_bias = (ω_δ - δ*ω_δ_prime)/2
    unit_var_proxy = ω_δ_prime^2

    n = nobs(Z)
    SteinMinimaxEstimator(Z=Z,
                          prior_class=prior_class,
                          target=target,
                          δ=δ,
                          ω_δ=ω_δ,
                          ω_δ_prime=ω_δ_prime,
                          g1=g1,
                          g2=g2,
                          Q=stein,
                          max_bias=max_bias,
                          unit_var_proxy=unit_var_proxy,
                          n=n,
                          δ_opt = δ,
                          δ_tuner = δ_tuner,
                          model=model)
end

target_bias_std(target::EBayesTarget, model, Zs; kwargs...) = target_bias_std(target, model; kwargs...)

"""
    StatsBase.confint(target::EBayesTarget, estimator, Zs;
					  level = 0.9, clip = true)

Construct confidence intervals for `target` based on samples `Zs`
with the estimator `estimator`. The nominal coverage is set to `level`. If `clip=true` (Default)
and `target` is bounded (say in ``[0,1]``), then the confidence intervals are clipped so
as to  satisfy these bounds.

Some of the estimators have already been fitted on the `Zs`, so that the following
call is also possible:

    StatsBase.confint(target::EBayesTarget, estimator;
                      level = 0.9, clip = true)
"""
confint(target::EBayesTarget, model, Zs; kwargs...) = confint(target, model; kwargs...)

function target_bias_std(target::EBayesTarget,
	                     sme::SteinMinimaxEstimator,
						 Zs=sme.Z;
                         std_type = :empirical)
	Qs = sme.(Zs) # assume Gaussian samples?
    #@show length(Qs)
	estimated_target = mean(Qs)
    if std_type == :empirical
	    estimated_std = std(Qs)/sqrt(nobs(Zs))
    elseif std_type == :var_proxy
        estimated_std = sqrt(sme.unit_var_proxy / nobs(Zs))
    else
        error("only :empirical and :var_proxy sdt_type defined")
    end

	estimated_bias = worst_case_bias(sme)
	(estimated_target = estimated_target,
	 estimated_bias = estimated_bias,
	 estimated_std = estimated_std)
end

function (sme::SteinMinimaxEstimator)(x)
	sme.Q(x)
end

function Distributions.estimate(target::EBayesTarget, fitted_model, args...; kwargs...)
	target_bias_std(target, fitted_model, args...; kwargs...)[:estimated_target]
end
#
#	end

function StatsBase.confint(target::EBayesTarget,
	                       sme, #SteinMinimaxEstimator #default fallback
                           args...;
						   level = 0.9,
						   clip = true,
                           kwargs...)
	res = target_bias_std(target, sme, args...; kwargs...)
	maxbias = abs(res[:estimated_bias])
	q_mult = bias_adjusted_gaussian_ci(res[:estimated_std], maxbias=maxbias , level=level)
    L,U = res[:estimated_target] .+ (-1,1).*q_mult
	L, U = clip ? clamp.((L,U), extrema(target)... ) : (L,U)
    L,U
end

@userplot SteinMinimaxPlot

@recipe function f(h::SteinMinimaxPlot; x_grid = -5:0.01:5,
	                                          ylim_calibrator = nothing,
											  ylim_relative_offset = 1.25,
                                              ylim_panel_e = nothing)
	arg_length = length(h.args)
	arg_length <= 2 || error("at most two plots allowed")
    sme = first(h.args)

	if arg_length == 2
		sme_band = last(h.args)
		layout := @layout [a{0.45w} grid(2,2)]
	else
		layout := @layout [panel1 panel2 panel3]
	end

	def_size = (850,260)

	col_g1 = RGB(77/255,137/255,193/255)
	col_gm1 = RGB(73/255,2/255,53/255)
	col_band = "#F6CD7F"#RGB(243/255,232/255,176/255)
	foreground_color_legend --> :transparent
	background_color_legend --> :transparent

    @series begin
        seriestype := :path
        subplot := 1
		linestyle := :dash
		if !isnothing(ylim_calibrator)
			ylims := ylim_calibrator
		end
        xlims := extrema(x_grid)
	    seriescolor := RGB(205/255,84/255,150/255)	# RGB(68/255,60/255,145/255)
		label := "Minimax"
        sme.Q
    end

    # Important note: In implementation I am doing L(G_2) - L(G_1)
    # So:     G_2 corresponds to G_1 from the paper
    #         G_1 corresponds to G_{-1} from the paper
	@series begin
		seriestype := :path
		subplot := 1
		linestyle := :solid
		seriescolor := RGB(132/255,193/255,216/255)
		xlims:= extrema(x_grid)
	    label := L"Minimax-$\infty$"
        title := "a)"
		legend := :topright
        xguide := L"x"
		sme_band.Q
	end

    # upper histogram
    g1 = sme.g1
    g2 = sme.g2
    g1_xs = pdf.(g1, x_grid)
    g2_xs = pdf.(g2, x_grid)
	ylim_prior = maximum([g1_xs;g2_xs])*ylim_relative_offset
    @series begin
        seriestype := :path
        subplot := 2
        linecolor := [col_g1 col_gm1]
        label :=  [L"g_1" L"g_{-1}"]
        xlims := extrema(x_grid)
		ylims := (0, ylim_prior)
        title := "b)"
        xguide := L"\mu"
		legend := :topright
        x_grid, [g2_xs g1_xs]
    end

	if (arg_length == 2)
		g1_band = sme_band.g1
		g2_band = sme_band.g2
		g1_xs = pdf.(g1_band, x_grid)
		g2_xs = pdf.(g2_band, x_grid)
		@series begin
			seriestype := :path
			subplot := 4
			linecolor := [col_g1 col_gm1]
			label :=  [L"g_1" L"g_{-1}"]
            xlims := extrema(x_grid)
			ylims := (0, maximum([g1_xs;g2_xs])*ylim_relative_offset)
            title := "d)"
			legend := :topright
            xguide := L"\mu"
			x_grid, [g2_xs g1_xs]
		end
	end

	# hardcoding for now to only work with StandardNormal
    Z=sme.Z
    Z_continuous = StandardNormalSample(0.0)

    f1=marginalize(g1, Z_continuous)
    f2=marginalize(g2, Z_continuous)

    f1_xs=pdf.(f1, x_grid)
    f2_xs=pdf.(f2, x_grid)

	ylim_marginal = maximum([f1_xs;f2_xs])*ylim_relative_offset
    @series begin
        seriestype := :path
        subplot := 3
        #bottom_margin := Plots.Measures.Length(:mm,5.0)
        xlims := extrema(x_grid)
		label :=  [L"f_{G_{1}}" L"f_{G_{-1}}"]#[L"\varphi  \star g_1" L"\varphi \star g_{-1} "]
        linecolor --> [col_g1 col_gm1]
		ylims --> (0, ylim_marginal)
        title := "c)"
		legend := :topright
        xguide := L"x"
        x_grid, [f2_xs f1_xs]
    end

	if (arg_length == 2)
		f1_band=marginalize(g1_band, Z_continuous)
		f2_band=marginalize(g2_band, Z_continuous)

		f1_xs=pdf.(f1_band, x_grid)
		f2_xs=pdf.(f2_band, x_grid)

        xs_mid = midpoints(sme_band.Z.mhist.grid)
        h = step(sme_band.Z.mhist.grid)
        ys_lower = sme_band.Z.f_min[2:(end-1)] ./ h
        ys_upper = sme_band.Z.f_max[2:(end-1)] ./ h

		#@series begin
		#	seriestype := :path
		#	subplot := 5
		 #   label :=  [L"f_{G_{1}}" L"f_{G_{-1}}"]#[L"\varphi  \star g_1" L"\varphi \star g_{-1} "]
		#	linecolor := [col_g1 col_gm1]
		#	x_grid, [f2_xs f1_xs]
		#end
        @series begin
            seriestype := :path
             subplot := 5
            # hack for now
             fillalpha := 1
             x := xs_mid
             y := ys_upper
             fillrange := ys_lower
             seriescolor := col_band
             label := L"$\infty$-band"
             ()
        end

        @series begin
        	seriestype := :path
        	subplot := 5
            label :=  L"f_{G_{1}}"#[L"\varphi  \star g_1" L"\varphi \star g_{-1} "]
         	linecolor := col_g1
        	x_grid, f2_xs
        end

        @series begin
          seriestype := :path
          subplot := 5
          if isnothing(ylim_panel_e)
             ylims := (0, ylim_marginal)
          else
             ylims := ylim_panel_e
          end
          xlims := extrema(x_grid)
          title := "e)"
          #title_location := :left
          legend := :topright
          xguide := L"x"
          label := L"f_{G_{-1}}"#[L"\varphi  \star g_1" L"\varphi \star g_{-1} "]
          linecolor := col_gm1
          x_grid, f1_xs
        end
    end
end

##marginal_pdf(gm_class, params) ->  #

#struct LipschitzPriorClass <: ConvexPriorClass
#end

# parametrize by slopes? easier constraints but harder to ....

#function check_bias(ma::MinimaxCalibrator; maximization=true)
#    ds = ma.ds
#    f = ma.f
#    C = ma.C
#    ε = ma.ε_reg
#    m = ma.m
#    target = ma.target

#    Q = ma.Q

#    check_bias(Q, ds, f, m, target; C=C, maximization=maximization)
#end


text:
struct SurvivalCopula{d,CT,VI} <: Copula{d}
    C::CT
    indices::VI
    function SurvivalCopula(C,indices)
        if length(indices) == 0
            return C
        end
        d = length(C)
        @assert all(indices .<= d)
        return new{d,typeof(C),typeof(indices)}(C,indices)
    end
end
function reverse!(u,idx)
    if ndims(u) == 1
        for i in idx
            u[i] = 1 - u[i]
        end
    else
        for i in idx
            u[i,:] .= 1 .- u[i,:]
        end
    end
    return u
end
function reverse(u,idx)
    v = deepcopy(u)
    reverse!(v,idx)
    return v
end
Distributions.cdf(C::SurvivalCopula{d,CT,VI},u) where {d,CT,VI} = Distributions.cdf(C.C,reverse(u,C.indices))
Distributions._logpdf(C::SurvivalCopula{d,CT,VI},u) where {d,CT,VI} = Distributions._logpdf(C.C,reverse(u,C.indices))
function Distributions._rand!(rng::Distributions.AbstractRNG, C::SurvivalCopula{d,CT,VI}, x::AbstractVector{T}) where {d,CT,VI,T}
    Distributions._rand!(rng,C.C,x)
    reverse!(x,C.indices)
end
function Base.rand(rng::Distributions.AbstractRNG,C::SurvivalCopula{d,CT,VI}) where {d,CT,VI}
    x = Base.rand(rnd,C.C)
    reverse!(x,C.indices)
    return x
end
function Distributions.fit(T::Type{CT},u) where {CT <: SurvivalCopula}
    # d = size(u,1)
    d,subCT,indices = T.parameters
    subfit = Distributions.fit(subCT,reverse(u,indices))
    return SurvivalCopula(subfit,indices)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_textarea

"""
    dcc_textarea(;kwargs...)

A Textarea component.
A basic HTML textarea for entering multiline text.
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `value` (String; optional): The value of the textarea
- `autoFocus` (String; optional): The element should be automatically focused after the page loaded.
- `cols` (String | Real; optional): Defines the number of columns in a textarea.
- `disabled` (String | Bool; optional): Indicates whether the user can interact with the element.
- `form` (String; optional): Indicates the form that is the owner of the element.
- `maxLength` (String | Real; optional): Defines the maximum number of characters allowed in the element.
- `minLength` (String | Real; optional): Defines the minimum number of characters allowed in the element.
- `name` (String; optional): Name of the element. For example used by the server to identify the fields in form submits.
- `placeholder` (String; optional): Provides a hint to the user of what can be entered in the field.
- `readOnly` (Bool | a value equal to: 'readOnly', 'readonly', 'READONLY'; optional): Indicates whether the element can be edited.
readOnly is an HTML boolean attribute - it is enabled by a boolean or
'readOnly'. Alternative capitalizations `readonly` & `READONLY`
are also acccepted.
- `required` (a value equal to: 'required', 'REQUIRED' | Bool; optional): Indicates whether this element is required to fill out or not.
required is an HTML boolean attribute - it is enabled by a boolean or
'required'. Alternative capitalizations `REQUIRED`
are also acccepted.
- `rows` (String | Real; optional): Defines the number of rows in a text area.
- `wrap` (String; optional): Indicates whether the text should be wrapped.
- `accessKey` (String; optional): Defines a keyboard shortcut to activate or add focus to the element.
- `className` (String; optional): Often used with CSS to style elements with common properties.
- `contentEditable` (String | Bool; optional): Indicates whether the element's content is editable.
- `contextMenu` (String; optional): Defines the ID of a <menu> element which will serve as the element's context menu.
- `dir` (String; optional): Defines the text direction. Allowed values are ltr (Left-To-Right) or rtl (Right-To-Left)
- `draggable` (a value equal to: 'true', 'false' | Bool; optional): Defines whether the element can be dragged.
- `hidden` (String; optional): Prevents rendering of given element, while keeping child elements, e.g. script elements, active.
- `lang` (String; optional): Defines the language used in the element.
- `spellCheck` (a value equal to: 'true', 'false' | Bool; optional): Indicates whether spell checking is allowed for the element.
- `style` (Dict; optional): Defines CSS styles which will override styles previously set.
- `tabIndex` (String | Real; optional): Overrides the browser's default tab order and follows the one specified instead.
- `title` (String; optional): Text to be displayed in a tooltip when hovering over the element.
- `n_blur` (Real; optional): Number of times the textarea lost focus.
- `n_blur_timestamp` (Real; optional): Last time the textarea lost focus.
- `n_clicks` (Real; optional): Number of times the textarea has been clicked.
- `n_clicks_timestamp` (Real; optional): Last time the textarea was clicked.
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `persistence` (Bool | String | Real; optional): Used to allow user interactions in this component to be persisted when
the component - or the page - is refreshed. If `persisted` is truthy and
hasn't changed from its previous value, a `value` that the user has
changed while using the app will keep that change, as long as
the new `value` also matches what was given originally.
Used in conjunction with `persistence_type`.
- `persisted_props` (Array of a value equal to: 'value's; optional): Properties whose user interactions will persist after refreshing the
component or the page. Since only `value` is allowed this prop can
normally be ignored.
- `persistence_type` (a value equal to: 'local', 'session', 'memory'; optional): Where persisted user changes will be stored:
memory: only kept in memory, reset on page refresh.
local: window.localStorage, data is kept after the browser quit.
session: window.sessionStorage, data is cleared once the browser quit.
"""
function dcc_textarea(; kwargs...)
        available_props = Symbol[:id, :value, :autoFocus, :cols, :disabled, :form, :maxLength, :minLength, :name, :placeholder, :readOnly, :required, :rows, :wrap, :accessKey, :className, :contentEditable, :contextMenu, :dir, :draggable, :hidden, :lang, :spellCheck, :style, :tabIndex, :title, :n_blur, :n_blur_timestamp, :n_clicks, :n_clicks_timestamp, :loading_state, :persistence, :persisted_props, :persistence_type]
        wild_props = Symbol[]
        return Component("dcc_textarea", "Textarea", "dash_core_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
function dijkstra_distances(src::Cell)
    dists = Dict{Cell,Int64}()
    dists[src] = 0
    frontier = [src]
    while !isempty(frontier)
        new_frontier = []
        for cell in frontier
            for next in cell.connections
                if !haskey(dists, next)
                    dists[next] = dists[cell] + 1
                    push!(new_frontier, next)
                end
            end
        end
        frontier = new_frontier
    end
    return dists
end

function dijkstra(src::Cell, target::Cell)
    dists = dijkstra_distances(src)
    current = target
    breadcrumbs = Dict{Cell,Int64}()
    breadcrumbs[current] = dists[current]
    while current != src
        for neighbor in current.connections
            if dists[neighbor] < dists[current]
                breadcrumbs[neighbor] = dists[neighbor]
                current = neighbor
                break
            end
        end
    end

    path = sort(collect(breadcrumbs); by=x -> x[2])
    return path, dists
end

function astar(src::Cell, goal::Cell, heuristic::Function)
    open_set = PriorityQueue()
    enqueue!(open_set, src, 0)
    came_from = Dict{Cell,Cell}()
    gscore = Dict{Cell,Int64}()
    gscore[src] = 0
    while !isempty(open_set)
        current = dequeue!(open_set)
        if current == goal
            break
        end
        for neighbor in current.connections
            tentative_gscore = gscore[current] + 1
            if !haskey(gscore, neighbor) || tentative_gscore < gscore[neighbor]
                came_from[neighbor] = current
                gscore[neighbor] = tentative_gscore
                neighbor_fscore = tentative_gscore + heuristic(neighbor, goal)
                if !haskey(open_set, neighbor)
                    enqueue!(open_set, neighbor, neighbor_fscore)
                end
            end
        end
    end

    return astar_reconstruct_path(came_from, gscore, goal), gscore
end

function astar_reconstruct_path(came_from, gscore, goal)
    path = Dict{Cell,Int64}()
    path[goal] = gscore[goal]
    current = goal
    while haskey(came_from, current)
        current = came_from[current]
        path[current] = gscore[current]
    end
    path = sort(collect(path); by=x -> x[2])
    return path
end


text:
#!/usr/bin/env julia

numeral_lengths =
    Dict(0 => 0,
         1 => 3,
         2 => 3,
         3 => 5,
         4 => 4,
         5 => 4,
         6 => 3,
         7 => 5,
         8 => 5,
         9 => 4,
         10 => 3,
         11 => 6,
         12 => 6,
         13 => 8,
         14 => 8,
         15 => 7,
         16 => 7,
         17 => 9,
         18 => 8,
         19 => 8,
         20 => 6,
         30 => 6,
         40 => 5,
         50 => 5,
         60 => 5,
         70 => 7,
         80 => 6,
         90 => 6,
         100 => 7,
         1000 => 8)

function number_length(n::Int64)
    assert(n < 10_000)
    result = 0
    (dig, n) = divrem(n, 1000)
    result += numeral_lengths[dig]
    if dig > 0
        result += numeral_lengths[1000]
    end
    (dig, n) = divrem(n, 100)
    result += numeral_lengths[dig]
    if dig > 0
        result += numeral_lengths[100]
        result += n > 0 ? 3 : 0
    end
    if n > 19
        (dig, n) = divrem(n, 10)
        result += numeral_lengths[dig*10]
        result += numeral_lengths[n]
    else
        result += numeral_lengths[n]
    end

    return result
end

function main()
    max_len = 0
    num = 0
    for i in 1:1000
        num += number_length(i)
    end
    @printf("The number of letters is %d\n", num)
end

main()


text:
##
## MvNormal
##

using Turing, ContinuousBenchmarks.TuringTools

# Define the target distribution and its gradient
const D = 10

@model target(dim) = begin
   Θ = Vector{Real}(undef, dim)
   θ ~ MvNormal(zeros(D), ones(dim))
end

# Sampling parameter settings
n_samples = 100_000
n_adapts = 2_000

# Sampling
bench_res = @tbenchmark_expr("NUTS(Leapfrog(...))",
                             sample(target(D), HMC(n_samples, 0.1, 5)));

LOG_DATA = build_log_data("MvNormal-Benchmark", bench_res...)
print_log(LOG_DATA)

##
## MvNormal: ForwardDiff vs BackwardDiff (Tracker)
##

using Turing, LinearAlgebra
using BenchmarkTools

@model mdemo(d, N) = begin
    Θ = Vector(undef, N)
   for n=1:N
      Θ[n] ~ d
   end
end

dim2 = 250
A   = rand(Wishart(dim2, Matrix{Float64}(I, dim2, dim2)))
d   = MvNormal(zeros(dim2), A)

# ForwardDiff
Turing.setadbackend(:forward_diff)
@benchmark chain = sample(mdemo(d, 1), HMC(5000, 0.1, 5))

#BackwardDiff
Turing.setadbackend(:reverse_diff)
@benchmark chain = sample(mdemo(d, 1), HMC(5000, 0.1, 5))

# build log and send data back to github.


text:
Random.seed!(7)

# model
T = Float64
integration = :rk4
model = Dynamics.quadrotor
model_d = discretize_model(model,integration)
n = model.n; m = model.m

# cost
Q = (1.0e-2)*Diagonal(I,n)
R = (1.0e-2)*Diagonal(I,m)
Qf = 1000.0*Diagonal(I,n)

# -initial state
x0 = zeros(n)
x0[1:3] = [0.; 0.; 0.]
q0 = [1.;0.;0.;0.]
x0[4:7] = q0

# -final state
xf = copy(x0)
xf[1:3] = [0.;50.;0.] # xyz position
xf[4:7] = q0

# options
verbose=false
opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,cost_tolerance=1.0e-5)
opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,opts_uncon=opts_ilqr,constraint_tolerance=1.0e-3,cost_tolerance=1.0e-5,cost_tolerance_intermediate=1e-4)
opts_altro = ALTROSolverOptions{T}(verbose=verbose,opts_al=opts_al)

N = 101
dt = 0.05
U0 = [0.5*9.81/4.0*ones(m) for k = 1:N-1]
obj = TrajectoryOptimization.LQRObjective(Q,R,Qf,xf,N)

# unconstrained
prob = Problem(model_d, obj, x0=x0, N=N, dt=dt)
initial_controls!(prob, U0)
solve!(prob, opts_ilqr)
@test norm(prob.X[N] - xf) < 5.0e-3

# constrained w/ final position
goal_con = goal_constraint(xf)
con = [goal_con]
prob = Problem(model_d, obj,constraints=Constraints(con,N), x0=x0, N=N, dt=dt)
initial_controls!(prob, U0)
solve!(prob, opts_al)
@test norm(prob.X[N] - xf,Inf) < opts_al.constraint_tolerance
@test max_violation(prob) < opts_al.constraint_tolerance

# constrained w/ final position and control limits
bnd = BoundConstraint(n,m,u_min=0.0,u_max=15.0,trim=true)
con = [bnd,goal_con]
prob = Problem(model_d, obj, constraints=Constraints(con,N), x0=x0, N=N, dt=dt)
initial_controls!(prob, U0)
solve!(prob, opts_al)
@test norm(prob.X[N] - xf) < opts_al.constraint_tolerance
@test max_violation(prob) < opts_al.constraint_tolerance

# constrained w/ final position, control limits, obstacles
r_quad = 1.0
r_sphere = 3.0
spheres = ((0.,10.,0.,r_sphere),(0.,20.,0.,r_sphere),(0.,30.,0.,r_sphere))
n_spheres = 3

function sphere_obs3(c,x,u)
    for i = 1:n_spheres
        c[i] = TrajectoryOptimization.sphere_constraint(x,spheres[i][1],spheres[i][2],spheres[i][3],spheres[i][4]+r_quad)
    end
    return nothing
end

obs = Constraint{Inequality}(sphere_obs3,n,m,n_spheres,:obs)
con = [bnd,obs,goal_con]
prob_con = Constraints(con,N)
prob = Problem(model_d, obj, constraints=Constraints(con,N),x0=x0, N=N, dt=dt)
initial_controls!(prob, U0)
opts_al.constraint_tolerance=1.0e-3
opts_al.constraint_tolerance_intermediate=1.0e-3
solve!(prob, opts_al)
@test norm(prob.X[N] - xf,Inf) < opts_al.constraint_tolerance
@test max_violation(prob) < opts_al.constraint_tolerance

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using Burgers
using Test

@testset "Burgers" begin
    xs, ys = Burgers.get_data(n=1000)

    @test size(xs) == (2, 1024, 1000)
    @test size(ys) == (1, 1024, 1000)

    learner = Burgers.train(epochs=10)
    loss = learner.cbstate.metricsepoch[ValidationPhase()][:Loss].values[end]
    @test loss < 0.1

    # include("deeponet.jl")
end


text:
"""
$(FUNCTIONNAME)(mod::LiquidSaltWater, T::FT,f::FT)

Computes the complex dieletric of liquid salty walter (Ulaby & Long book)
# Arguments
- `mod` an [`LiquidSaltWater`](@ref) type struct, provide Salinity in PSU
- `T`  Temperature in `[⁰K]`
- `f`  Frequency in `[GHz]`
- `S`  Salinity in `[PSU]` comes from the `mod` [`LiquidSaltWater`](@ref) struct

# Examples
```julia-repl
julia> w = CanopyOptics.LiquidSaltWater(S=10.0)     # Create struct for salty seawater
julia> CanopyOptics.dielectric(w,283.0,10.0)
51.79254073931596 + 38.32304044382495im
```
"""
function dielectric(mod::LiquidSaltWater, T::FT,f::FT) where FT<:Real
    @unpack S  = mod
    @assert 0 ≤ S ≤ 45 "Salinity should be ∈ [0,45] PSU"
    @assert 265 ≤ T ≤ 310 "Temperature should be ∈ [265,310] K"
    # Equations were tuned for ⁰C but want K input!
    T = T - FT(273)
    # Compute σ (Eqs 4.21 - 4.21 in Ulaby and Long)
    σ₃₅ = TempPoly(T)
    P   = S * SalPoly1(S) / SalPoly2(S)
    Q   = FT(1) + (α₀PolNom(S)/α₀PolDenom(S) * (T - FT(15)))/(T + α₁(T));

    σ = σ₃₅ * P * Q

    ϵS   = FT(87.85306) * exp(FT(-0.00456992)*T - aU[1]*S - aU[2]*S^2 - aU[3]*S*T);
    ϵOne =  aU[4] * exp(-aU[5]*T - aU[6]*S - aU[7]*S*T);
    # Relaxation time constants
    τ1   = (aU[8]  + aU[9] *S) * exp(aU[10] / (T+aU[11]));
    τ2   = (aU[12] + aU[13]*S) * exp(aU[14] / (T+aU[15]));
    ϵInf =  aU[16] + aU[17]*T  + aU[18]*S;

    # Complex Permittivity Calculation
    ϵ = ((ϵS-ϵOne) / (FT(1) - 1im * 2π * f * τ1)) + ((ϵOne-ϵInf) / (1-1im * 2π * f * τ2)) + (ϵInf) + 1im*((FT(17.9751)*σ) /f);
end

"""
$(FUNCTIONNAME)(mod::LiquidPureWater, T::FT,f::FT)

Computes the complex dieletric of liquid pure walter (Ulaby & Long book)
# Arguments
- `mod` an [`LiquidSaltWater`](@ref) type struct
- `T`  Temperature in `[⁰K]`
- `f`  Frequency in `[GHz]`

# Examples
```julia-repl
julia> w = CanopyOptics.LiquidPureWater()     # Create struct for salty seawater
julia> CanopyOptics.dielectric(w,283.0,10.0)
53.45306674215052 + 38.068642430090044im
```
"""
function dielectric(mod::LiquidPureWater, T::FT,f::FT) where FT<:Real
    @assert 265 ≤ T ≤ 310 "Temperature should be ∈ [265,310] K"
    # Equations were tuned for ⁰C but want K input!
    T = T - FT(273)

    ϵS   =  FT(87.85306) * exp(FT(-0.00456992)*T);
    ϵOne =  aU[4]  * exp(-aU[5]*T);
    τ1   =  aU[8]  * exp(aU[10] / (T+aU[11]));
    τ2   =  aU[12] * exp(aU[14] / (T+aU[15]));
    ϵInf =  aU[16] + aU[17]*T;
    # Complex Permittivity Calculation
    ϵ = ((ϵS-ϵOne) / (FT(1) - 1im * 2π * f * τ1)) + ((ϵOne-ϵInf)./(1-1im * 2π * f * τ2)) + (ϵInf);
end

"""
$(FUNCTIONNAME)(mod::PureIce, T::FT,f::FT)

Computes the complex dieletric of liquid pure walter (Ulaby & Long book)
# Arguments
- `mod` an [`PureIce`](@ref) type struct
- `T`  Temperature in `[⁰K]`
- `f`  Frequency in `[GHz]`

# Examples
```julia-repl
julia> w = CanopyOptics.PureIce()     # Create struct for salty seawater
julia> CanopyOptics.dielectric(w,283.0,10.0)
53.45306674215052 + 38.068642430090044im
```
"""
function dielectric(mod::PureIce, T::FT, f::FT) where FT
    @assert 233 ≤ T ≤ 273.15 "Temperature should be ∈ [233,273.15] K"
    # Some constants first
    θ  = (FT(300)/T) - FT(1);
    B1 = FT(0.0207);
    B2 = FT(1.16e-11);
    b  = FT(335);
    # Eqs 4.22-4.25
    α = (FT(0.00504) + FT(0.0062) * θ) * exp(-FT(22.1) * θ);
    βₘ  = (B1/T) * exp(b/T) / (exp(b/T)-FT(1))^2  + B2 * f^2;
    δβ  = exp(-FT(9.963) + FT(0.0372) * (T-FT(273.16)));

    β = βₘ + δβ;
    ϵ = FT(3.1884) + FT(9.1e-4) *(T-FT(273)) + (α/f + β*f)im;
end

"""
$(FUNCTIONNAME)(mod::SoilMW, T::FT,f::FT)

Computes the complex dieletric of soil (Ulaby & Long book)
# Arguments
- `mod` an [`SoilMW`](@ref) type struct (includes sand_frac, clay_frac, water_frac, ρ)
- `T`  Temperature in `[⁰K]`
- `f`  Frequency in `[GHz]`

# Examples
```julia-repl
julia> w = SoilMW(sand_frac=0.2,clay_frac=0.2, water_frac = 0.3, ρ=1.7 )     # Create struct for salty seawater
julia> dielectric(w,283.0,10.0)
53.45306674215052 + 38.068642430090044im
```
"""
function dielectric(mod::SoilMW{FT}, T::FT,f::FT) where FT
    @unpack sand_frac, clay_frac, mᵥ, ρ = mod
    # T = Tk - FT(273)
    f_hz = f * FT(1e9);    #  transform from GHz to Hz

    α  = FT(0.65);   # eq: 4.68a
    β₁ = FT(1.27) - FT(0.519) * sand_frac - FT(0.152) * clay_frac; # eq: 4.68b
    β₂ = FT(2.06) - FT(0.928) * sand_frac - FT(0.255) * clay_frac; # eq: 4.68c

    ϵ₀ = FT(8.854e-12);

    if f > 1.3
        σ = FT(-1.645) + FT(1.939) * ρ - FT(2.256) * sand_frac + FT(1.594) * clay_frac;  # eq: 4.68d
    elseif 0.3 ≤ f ≤ 1.3
        σ = FT(0.0467) + FT(0.22)  * ρ - FT(0.411) * sand_frac + FT(0.661) * clay_frac;  # eq: 4.70
    else
        σ = FT(0)
    end

    # Dielectric Constant of Pure Water
    ϵw =  dielectric(LiquidPureWater(),T,f)
    # Add conductivity term to ϵ'' (eq. 4.67b)
    ϵw += im*((FT(2.65)-ρ) / FT(2.65 * mᵥ) * σ / (2π * ϵ₀ * f_hz))

    # calculating dielectric constant of soil using eq. 4.66a and 4.66b
    epsr = (FT(1) + FT(0.66) * ρ + mᵥ^β₁ * real(ϵw).^α - mᵥ)^(1/α);
    epsi = mᵥ^β₂ .* imag(ϵw);
    epsr + epsi*im
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 55518988-20c3-11eb-02b0-375c3504f7b1
#this notebook shows Julia type hierarchy

# ╔═╡ dd424284-20c0-11eb-344a-9b14b7570882
begin
	abstract type Asset end;
	abstract type Property <:Asset end;
	abstract type Investment <:Asset end;
	abstract type Cash <: Asset end;
	abstract type Equity <: Investment end;
end

# ╔═╡ 6b99fad0-20c2-11eb-18c3-7511057eef2a
@show subtypes(Asset)

# ╔═╡ ae14a644-20c2-11eb-00ea-5f8ed7d5a6b7
@show supertype(Cash)

# ╔═╡ e9b32cde-20c2-11eb-3903-e765cac8b0f0
@show supertype(Equity)

# ╔═╡ 0a79941e-20c3-11eb-0f77-7b28e2b48560
@show subtypes(Investment)

# ╔═╡ Cell order:
# ╠═55518988-20c3-11eb-02b0-375c3504f7b1
# ╠═dd424284-20c0-11eb-344a-9b14b7570882
# ╠═6b99fad0-20c2-11eb-18c3-7511057eef2a
# ╠═ae14a644-20c2-11eb-00ea-5f8ed7d5a6b7
# ╠═e9b32cde-20c2-11eb-3903-e765cac8b0f0
# ╠═0a79941e-20c3-11eb-0f77-7b28e2b48560


text:
doc"""
    GeneralizedPareto(ξ, σ, μ)

The *Generalized Pareto distribution* with shape parameter `ξ`, scale `σ` and location `μ` has probability density function

$f(x; \xi, \sigma, \mu) = \begin{cases}
        \frac{1}{\sigma}(1 + \xi \frac{x - \mu}{\sigma} )^{-\frac{1}{\xi} - 1} & \text{for } \xi \neq 0 \\
        \frac{1}{\sigma} e^{-\frac{\left( x - \mu \right) }{\sigma}} & \text{for } \xi = 0
    \end{cases}~,
    \quad x \in \begin{cases}
        \left[ \mu, \infty \right] & \text{for } \xi \geq 0 \\
        \left[ \mu, \mu - \sigma / \xi \right] & \text{for } \xi < 0
    \end{cases}$

```julia
GeneralizedPareto()             # Generalized Pareto distribution with unit shape and unit scale, i.e. GeneralizedPareto(1.0, 1.0, 0.0)
GeneralizedPareto(k, s)         # Generalized Pareto distribution with shape k and scale s, i.e. GeneralizedPareto(k, s, 0.0)
GeneralizedPareto(k, s, m)      # Generalized Pareto distribution with shape k, scale s and location m.

params(d)       # Get the parameters, i.e. (k, s, m)
shape(d)        # Get the shape parameter, i.e. k
scale(d)        # Get the scale parameter, i.e. s
location(d)     # Get the location parameter, i.e. m
```

External links

* [Generalized Pareto distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_Pareto_distribution)

"""

immutable GeneralizedPareto <: ContinuousUnivariateDistribution
    ξ::Float64
    σ::Float64
    μ::Float64

    function GeneralizedPareto(ξ::Real, σ::Real, μ::Real)
        @check_args(GeneralizedPareto, σ > zero(σ))
        new(ξ, σ, μ)
    end
    GeneralizedPareto(ξ::Real, σ::Real) = GeneralizedPareto(ξ::Real, σ::Real, 0.0)
    GeneralizedPareto() = new(1.0, 1.0, 0.0)
end

minimum(d::GeneralizedPareto) = d.μ
maximum(d::GeneralizedPareto) = d.ξ < 0.0 ? d.μ - d.σ / d.ξ : Inf

#### Parameters

shape(d::GeneralizedPareto) = d.ξ
scale(d::GeneralizedPareto) = d.σ
location(d::GeneralizedPareto) = d.μ
params(d::GeneralizedPareto) = (d.ξ, d.σ, d.μ)

#### Statistics

median(d::GeneralizedPareto) = d.μ + d.σ * expm1(d.ξ * log(2.0)) / d.ξ

function mean(d::GeneralizedPareto)
    if d.ξ < 1.0
        return d.μ + d.σ / (1.0 - d.ξ)
    else
        return Inf
    end
end

function var(d::GeneralizedPareto)
    if d.ξ < 0.5
        return d.σ^2 / ((1.0 - d.ξ)^2 * (1.0 - 2.0 * d.ξ))
    else
        return Inf
    end
end

function skewness(d::GeneralizedPareto)
    (ξ, σ, μ) = params(d)

    if ξ < (1.0 / 3.0)
        return 2.0 * (1.0 + ξ) * sqrt(1.0 - 2.0 * ξ) / (1.0 - 3.0 * ξ)
    else
        return Inf
    end
end

function kurtosis(d::GeneralizedPareto)
    (ξ, σ, μ) = params(d)

    if ξ < 0.25
        k1 = (1.0 - 2.0 * ξ) * (2.0 * ξ^2 + ξ + 3.0)
        k2 = (1.0 - 3.0 * ξ) * (1.0 - 4.0 * ξ)
        return 3.0 * k1 / k2 - 3.0
    else
        return Inf
    end
end

#### Evaluation

function logpdf(d::GeneralizedPareto, x::Float64)
    (ξ, σ, μ) = params(d)

    # The logpdf is log(0) outside the support range.
    p = -Inf

    if x >= μ
        z = (x - μ) / σ
        if abs(ξ) < eps()
            p = -z - log(σ)
        elseif ξ > 0.0 || (ξ < 0.0 && x < maximum(d))
            p = (-1.0 - 1.0 / ξ) * log1p(z * ξ) - log(σ)
        end
    end

    return p
end

pdf(d::GeneralizedPareto, x::Float64) = exp(logpdf(d, x))

function logccdf(d::GeneralizedPareto, x::Float64)
    (ξ, σ, μ) = params(d)

    # The logccdf is log(0) outside the support range.
    p = -Inf

    if x >= μ
        z = (x - μ) / σ
        if abs(ξ) < eps()
            p = -z
        elseif ξ > 0.0 || (ξ < 0.0 && x < maximum(d))
            p = (-1.0 / ξ) * log1p(z * ξ)
        end
    end

    return p
end

ccdf(d::GeneralizedPareto, x::Float64) = exp(logccdf(d, x))
cdf(d::GeneralizedPareto, x::Float64) = -expm1(logccdf(d, x))

function quantile(d::GeneralizedPareto, p::Float64)
    (ξ, σ, μ) = params(d)

    if p == 0.0
        z = 0.0
    elseif p == 1.0
        z = ξ < 0.0 ? -1.0 / ξ : Inf
    elseif 0.0 < p < 1.0
        if abs(ξ) < eps()
            z = -log1p(-p)
        else
            z = expm1(-ξ * log1p(-p)) / ξ
        end
    else
      z = NaN
    end

    return μ + σ * z
end

#### Sampling

function rand(d::GeneralizedPareto)
    # Generate a Float64 random number uniformly in (0,1].
    u = 1.0 - rand()

    if abs(d.ξ) < eps()
        rd = -log(u)
    else
        rd = expm1(-d.ξ * log(u)) / d.ξ
    end

    return d.μ + d.σ * rd
end


text:
module MultiThreadedCaches

using Base: @lock

export MultiThreadedCache, init_cache!

"""
    MultiThreadedCache{K,V}()
    MultiThreadedCache{K,V}([initial_kv_values])

`MultiThreadedCache{K,V}()` constructs a fast, thread-safe cache.

This cache stores k=>v pairs that cache a deterministic computation. The only API into the
cache is `get!()`: you can look up a key, and if it is not available, you can produce a
value which will be added to the cache.

Accesses to the cache will look first in the per-thread cache, and then fall back to the
shared thread-safe cache. Concurrent misses to the same key in the shared cache will
coordinate, so that only one Task will perform the compuatation for that value, and the
other Task(s) will block.

# Examples:
```julia
julia> cache = MultiThreadedCache{Int, Int}(Dict(1=>2, 2=>3))
MultiThreadedCache{Int64, Int64}(Dict(2 => 3, 1 => 2))

julia> init_cache!(cache)
MultiThreadedCache{Int64, Int64}(Dict(2 => 3, 1 => 2))

julia> get!(cache, 2) do
           2+1
       end
3

julia> get!(cache, 5) do
           5+1
       end
6

julia> get!(cache, 5) do
           5+10
       end
6
```
"""
struct MultiThreadedCache{K,V}
    thread_caches::Vector{Dict{K,V}}
    thread_locks::Vector{ReentrantLock}

    base_cache::Dict{K,V}  # Guarded by: base_cache_lock
    base_cache_lock::ReentrantLock
    base_cache_futures::Dict{K,Channel{V}}  # Guarded by: base_cache_lock

    function MultiThreadedCache{K,V}() where {K,V}
        base_cache = Dict{K,V}()

        return MultiThreadedCache{K,V}(base_cache)
    end

    function MultiThreadedCache{K,V}(base_cache::Dict) where {K,V}
        thread_caches = Dict{K,V}[]
        thread_locks = ReentrantLock[]

        base_cache_lock = ReentrantLock()
        base_cache_futures = Dict{K,Channel{V}}()

        return new(thread_caches, thread_locks, base_cache, base_cache_lock, base_cache_futures)
    end
end

"""
    init_cache!(cache::MultiThreadedCache{K,V})

This function must be called inside a Module's `__init__()` method, to ensure that the
number of threads are set *at runtime, not precompilation time*.

!!! note
NOTE: This function is *not thread safe*, it must not be called concurrently with any other
code that touches the cache. This should only be called once, during module initialization.
"""
function init_cache!(cache::MultiThreadedCache{K,V}) where {K,V}
    # Statically resize the vector, but wait to lazily create the dictionaries when
    # requested, so that the object will be allocated on the thread that will consume it.
    # (This follows the guidance from Julia Base.)
    resize!(cache.thread_caches, Threads.nthreads())
    resize!(cache.thread_locks, Threads.nthreads())
    return cache
end

# Based upon the thread-safe Global RNG implementation in the Random stdlib:
# https://github.com/JuliaLang/julia/blob/e4fcdf5b04fd9751ce48b0afc700330475b42443/stdlib/Random/src/RNGs.jl#L369-L385
# Get or lazily construct the per-thread cache when first requested.
function _thread_cache(mtcache::MultiThreadedCache, tid)
    length(mtcache.thread_caches) >= Threads.nthreads() || _thread_cache_length_assert()
    if @inbounds isassigned(mtcache.thread_caches, tid)
        @inbounds cache = mtcache.thread_caches[tid]
    else
        # We copy the base cache to all the thread caches, so that any precomputed values
        # can be shared without having to wait for a cache miss.
        cache = Base.@lock mtcache.base_cache_lock begin
            copy(mtcache.base_cache)
        end
        @inbounds mtcache.thread_caches[tid] = cache
    end
    return cache
end
@noinline _thread_cache_length_assert() = @assert false "** Must call `init_cache!(cache)` in your Module's __init__()! - length(cache.thread_caches) < Threads.nthreads() "
function _thread_lock(cache::MultiThreadedCache, tid)
    length(cache.thread_locks) >= Threads.nthreads() || _thread_cache_length_assert()
    if @inbounds isassigned(cache.thread_locks, tid)
        @inbounds lock = cache.thread_locks[tid]
    else
        lock = eltype(cache.thread_locks)()
        @inbounds cache.thread_locks[tid] = lock
    end
    return lock
end

const CACHE_MISS = :__MultiThreadedCaches_key_not_found__

function Base.get!(func::Base.Callable, cache::MultiThreadedCache{K,V}, key) where {K,V}
    # If the thread-local cache has the value, we can return immediately.
    # We store tcache in a local variable, so that even if the Task migrates Threads, we are
    # still operating on the same initial cache object.
    tid = Threads.threadid()
    tcache = _thread_cache(cache, tid)
    tlock = _thread_lock(cache, tid)
    # We have to lock during access to the thread-local dict, because it's possible that the
    # Task may migrate to another thread by the end, and we really might be mutating the
    # dict in parallel. But most of the time this lock should have 0 contention, since it's
    # only held during get() and set!().
    Base.@lock tlock begin
        thread_local_cached_value_or_miss = get(tcache, key, CACHE_MISS)
        if thread_local_cached_value_or_miss !== CACHE_MISS
            return thread_local_cached_value_or_miss::V
        end
    end
    # If not, we need to check the base cache.
    # When we're done, call this function to set the result
    @inline function _store_result!(v::V; test_haskey::Bool)
        # Set the value into thread-local cache for the supplied key.
        # Note that we must perform two separate get() and setindex!() calls, for
        # concurrency-safety, in case the dict has been mutated by another task in between.
        # TODO: For 100% concurrency-safety, we maybe want to lock around the get() above
        # and the setindex!() here.. it's probably fine without it, but needs considering.
        # Currently this is relying on get() and setindex!() having no yields.
        Base.@lock tlock begin
            if test_haskey
                if !haskey(tcache, key)
                    tcache[key] = v
                end
            else
                tcache[key] = v
            end
        end
    end

    # Even though we're using Thread-local caches, we still need to lock during
    # construction to prevent multiple tasks redundantly constructing the same object,
    # and potential thread safety violations due to Tasks migrating threads.
    # NOTE that we only grab the lock if the key doesn't exist, so the mutex contention
    # is not on the critical path for most accessses. :)
    is_first_task = false
    local future  # used only if the base_cache doesn't have the key
    # We lock the mutex, but for only a short, *constant time* duration, to grab the
    # future for this key, or to create the future if it doesn't exist.
    @lock cache.base_cache_lock begin
        value_or_miss = get(cache.base_cache, key, CACHE_MISS)
        if value_or_miss !== CACHE_MISS
            return value_or_miss::V
        end
        future = get!(cache.base_cache_futures, key) do
            is_first_task = true
            Channel{V}(1)
        end
    end
    if is_first_task
        v = try
            func()
        catch e
            # In the case of an exception, we abort the current computation of this
            # key/value pair, and throw the exception onto the future, so that all
            # pending tasks will see the exeption as well.
            #
            # NOTE: we could also cache the exception and throw it from now on, but this
            # would make interactive development difficult, since once you fix the
            # error, you'd have to clear out your cache. So instead, we just rethrow the
            # exception and don't cache anything, so that you can fix the exception and
            # continue on. (This means that throwing exceptions remains expensive.)

            # close(::Channel, ::Exception) requires an Exception object, so if the user
            # threw a non-Exception, we convert it to one, here.
            e isa Exception || (e = ErrorException("Non-exception object thrown during get!(): $e"))
            close(future, e)
            # As below, the future isn't needed after this returns (see below).
            @lock cache.base_cache_lock begin
                delete!(cache.base_cache_futures, key)
            end
            rethrow(e)
        end
        # Finally, lock again for a *constant time* to insert the computed value into
        # the shared cache, so that we can free the Channel and future gets can read
        # from the shared base_cache.
        @lock cache.base_cache_lock begin
            cache.base_cache[key] = v
            # We no longer need the Future, since all future requests will see the key
            # in the base_cache. (Other Tasks may still hold a reference, but it will
            # be GC'd once they have all completed.)
            delete!(cache.base_cache_futures, key)
        end
        # Store the result in this thread-local dictionary.
        _store_result!(v, test_haskey=false)
        # Return v to any other Tasks that were blocking on this key.
        put!(future, v)
        return v
    else
        # Block on the future until the first task that asked for this key finishes
        # computing a value for it.
        v = fetch(future)
        # Store the result in our original thread-local cache (if another Task hasn't) set
        # it already.
        _store_result!(v, test_haskey=true)
        return v
    end
end

function Base.show(io::IO, cache::MultiThreadedCache{K,V}) where {K,V}
    # Contention optimization: don't hold the lock while printing, since that could block
    # for an arbitrarily long time. Instead, print the data to an intermediate buffer first.
    # Note that this has the same CPU complexity, since printing is already O(n).
    iobuf = IOBuffer()
    let io = IOContext(iobuf, io)
        Base.@lock cache.base_cache_lock begin
            _oneline_show(io, cache)
        end
    end
    # Now print the data without holding the lock.
    seekstart(iobuf)
    write(io, read(iobuf))
    return nothing
end
_oneline_show(io::IO, cache::MultiThreadedCache{K,V}) where {K,V} =
    print(io, "$(MultiThreadedCache{K,V})(", cache.base_cache, ")")

function Base.show(io::IO, mime::MIME"text/plain", cache::MultiThreadedCache{K,V}) where {K,V}
    # Contention optimization: don't hold the lock while printing. See above for more info.
    iobuf = IOBuffer()
    let io = IOContext(iobuf, io)
        Base.@lock cache.base_cache_lock begin
            if isempty(cache.base_cache)
                _oneline_show(io, cache)
            else
                print(io, "$(MultiThreadedCache): ")
                Base.show(io, mime, cache.base_cache)
            end
        end
    end
    # Now print the data without holding the lock.
    seekstart(iobuf)
    write(io, read(iobuf))
    return nothing
end

end # module


text:
module TestWorkStealing

using Folds
using FoldsThreads
using Test
using Transducers

@testset "many self-steals" begin
    @testset for basesize in [2^6, 2^8, 2^10]
        @testset for nchunks in [2^6, 2^8, 2^10]
            len = basesize * nchunks + 1
            xs = rand(typemin(UInt):typemax(UInt), len)
            @test Folds.reduce(
                xor,
                xs,
                WorkStealingEx(basesize = basesize);
                init = UInt(0),  # TODO: get rid of this
            ) == reduce(xor, xs)
            @test Folds.collect(xs |> Filter(isodd), WorkStealingEx(basesize = basesize)) ==
                  filter(isodd, xs)
        end
    end
end

function random_work(i)
    n = 0
    t0 = time_ns()
    if mod(t0, 32) == 0
        t1 = t0 + 100_000
        while t1 > time_ns()
        end
    end
    return i
end

@testset "random scheduling" begin
    @testset for _trial_ in 1:5
        xs = 1:1000
        @test Folds.collect(xs |> Map(random_work), WorkStealingEx(basesize = 1)) == xs
    end
end

end  # module


text:
"""
    circ_smooth(v, L=501)

Smooth a signal with a gaussian window with support over about `L`. Assumes the
signal is periodic.
"""
function circ_smooth(v, L=501)
    N = length(v)
    win = gaussian(N, 0.15*L/N; zerophase=true)
    win ./= sum(win)
    circ_conv(v, win)
end

"""
Perform circular convolution using the FFT algorithm. The result length will be the
greater of the arguments.
"""
function circ_conv(u, v)
    lu = length(u)
    lv = length(v)
    nfft = nextfastfft(max(lu, lv))
    # always copy `u` because we re-use this buffer for the result
    # TODO: add zero-padding to `pad` function and use here
    upad = similar(u, nfft)
    copyto!(upad, u)
    for i in lu+1:nfft
        @inbounds upad[i] = zero(eltype(u))
    end
    vpad = if lv != nfft
        similar(v, nfft)
        copyto!(vpad, v)
        for i in lv+1:nfft
            @inbounds vpad[i] = zero(eltype(v))
        end
    else
        v
    end
    p = plan_rfft(upad)
    U = p*upad
    V = p*vpad
    @. U *= V
    # I think the arguments for this are reversed in the FFTW docs
    ldiv!(upad, p, U)
    resize!(upad, max(lu,lv))
end

"""
    circ_resample(v, rate)

Performs time-domain resampling of `v` by the ratio `rate`, treating `v` as
periodic. When downsampling the signal is lowpassed before resampling.
"""
function circ_resample(v, rate)
    L = length(v)
    s = rate < 1 ? circ_smooth(v, 2/rate) : v
    map(0:round(Int, L*rate)-1) do i
        cubic(s, i/rate+1, :periodic)
    end
end

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
module FactorOperations

export Factor
export product
export marginalize
export divice

using JoinTreeInference: Node, Potential

type Factor
    nodes::Array{String, 1}
    cardinality::Array{Int, 1}
    stride::Array{Int, 1}
    data::Array{Float64, 1}
end

function Factor(nodes::Array{String, 1}, cardinality::Array{Int, 1}, data::Array{Float64, 1})
    n = length(nodes)
    stride = Array{Int, 1}(n)
    s = 1
    for i in 1:n
        stride[i] = s
        s *= cardinality[i]
    end
    return Factor(nodes, cardinality, stride, data)
end

function product(factors::Array{Factor, 1})
        num_factors = length(factors)
        vars = union([f.nodes for f in factors]...)
        var_ids = Dict(j=>i for (i,j) in enumerate(vars))

        n = length(vars)
        cardinality = Array{Int, 1}(n)
        stride = Array{Int, 1}(n)

        idmap = Array{Array{Int, 1}, 1}(num_factors)
        for (i, f) in enumerate(factors)
            nv = length(f.nodes)
            idmap[i] = fill(0, n)
            for (j, v) in enumerate(f.nodes)
                id = var_ids[v]
                cardinality[id] = f.cardinality[j]
                idmap[i][id] = j
            end
        end

        s = 1
        for i in 1:n
            stride[i] = s
            s *= cardinality[i]
        end

        assignments = fill(1, n)
        ind = fill(1, num_factors)
        num_rows = reduce(*, cardinality)
        data = fill(1.0, num_rows)
        for i in 1:num_rows
            for (j, f) in enumerate(factors)
                data[i] *= f.data[ind[j]]
            end
            for l in 1:n
                assignments[l] += 1
                if assignments[l] > cardinality[l]
                    assignments[l] = 1
                    for j in 1:num_factors
                        id = idmap[j][l]
                        if id > 0
                            ind[j] -= factors[j].stride[id] * (cardinality[l]-1)
                        end
                    end
                else
                    for j in 1:num_factors
                        id = idmap[j][l]
                        if id > 0
                            ind[j] += factors[j].stride[id]
                        end
                    end
                    break
                end
            end
        end
        return Factor(vars, cardinality, stride, data)
end

function marginalize(f::Factor, v::Set{String})
    fn = length(f.nodes)
    n = fn - length(v)

    vars = Array{String, 1}(n)
    cardinality = Array{Int, 1}(n)
    idmap = fill(0, fn)
    isin = fill(false, fn)

    j = 1
    for (i, var) in enumerate(f.nodes)
        if !(var in v)
            vars[j] = var
            cardinality[j] = f.cardinality[i]
            idmap[i] = j
            isin[i] = true
            j += 1
        end
    end

    stride = Array{Int, 1}(n)
    s = 1
    for i in 1:n
        stride[i] = s
        s *= cardinality[i]
    end

    data = fill(0.0, reduce(*, cardinality))
    num_rows = reduce(*, f.cardinality)
    assignments = fill(1, fn)
    ind = 1
    for i in 1:num_rows
        data[ind] += f.data[i]
        for (l, n) in enumerate(f.nodes)
            assignments[l] += 1
            if assignments[l] > f.cardinality[l]
                assignments[l] = 1
                if isin[l]
                    ind -= stride[idmap[l]] * (f.cardinality[l]-1)
                end
            else
                if isin[l]
                    ind += stride[idmap[l]]
                end
                break
            end
        end
    end
    return Factor(vars, cardinality, stride, data)
end

function divide(f1::Factor, f2::Factor)
    d = Dict((j, i) for (i,j) in enumerate(f2.nodes))
    v = Set(f2.nodes)

    fn = length(f1.nodes)
    idmap = fill(0, fn)
    isin = fill(false, fn)
    for (i, n) in enumerate(f1.nodes)
        if n in v
            isin[i] = true
            idmap[i] = d[n]
        end
    end

    assignments = fill(1, fn)
    f = deepcopy(f1)
    num_rows = reduce(*, f.cardinality)
    ind = 1
    for i in 1:num_rows
        f.data[i] /= f2.data[ind]
        for (l, n) in enumerate(f1.nodes)
            assignments[l] += 1
            if assignments[l] > f.cardinality[l]
                assignments[l] = 1
                if isin[l]
                    ind -= f.stride[idmap[l]] * (f.cardinality[l]-1)
                end
            else
                if isin[l]
                    ind += f.stride[idmap[l]]
                end
                break
            end
        end
    end
    return f
end

end

max_stars_repo_name: bhrzslm/uncertainty-reasoning

max_issues_repo_name: bhrzslm/uncertainty-reasoning

max_forks_repo_name: bhrzslm/uncertainty-reasoning


text:
# test index notation using @tensor macro
#-----------------------------------------
A=randn((3,5,4,6))
p=randperm(4)
C1=permutedims(A,p)
@eval @tensor C2[$(p...)] := A[1,2,3,4]
@test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
@test_throws TensorOperations.IndexError begin
    @tensor C[1,2,3,4] := A[1,2,3]
end
@test_throws TensorOperations.IndexError begin
    @tensor C[1,2,3,4] := A[1,2,2,4]
end

B=randn((5,6,3,4))
p=[3,1,4,2]
@tensor C1[3,1,4,2] := A[3,1,4,2] + B[1,2,3,4]
C2=A+permutedims(B,p)
@test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
@test_throws DimensionMismatch begin
    @tensor C[1,2,3,4] := A[1,2,3,4] + B[1,2,3,4]
end

A=randn(50,100,100)
@tensor C1[a] := A[a,b',b']
C2=zeros(50)
for i=1:50
    for j=1:100
        C2[i]+=A[i,j,j]
    end
end
@test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
A=randn(3,20,5,3,20,4,5)
@tensor C1[e,a,d] := A[a,b,c,d,b,e,c]
C2=zeros(4,3,3)
for i1=1:4, i2=1:3, i3=1:3
    for j1=1:20,j2=1:5
        C2[i1,i2,i3]+=A[i2,j1,j2,i3,j1,i1,j2]
    end
end
@test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)

A=randn(3,20,5,3,4)
B=randn(5,6,20,3)
@tensor C1[a,g,e,d,f] := A[a,b,c,d,e]*B[c,f,b,g]
C2=zeros(3,3,4,3,6)
for a=1:3, b=1:20, c=1:5, d=1:3, e=1:4, f=1:6, g=1:3
    C2[a,g,e,d,f] += A[a,b,c,d,e]*B[c,f,b,g]
end
@test_throws TensorOperations.IndexError begin
    @tensor A[a,b,c,d]*B[c,f,b,g]
end

A=randn(5,5,5,5)
B=rand(Complex128,(5,5,5,5))
@tensor C1[1,2,5,6,3,4,7,8] := A[1,2,3,4]*B[5,6,7,8]
C2=reshape(kron(reshape(B,(25,25)),reshape(A,(25,25))),(5,5,5,5,5,5,5,5))
@test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
@test_throws TensorOperations.IndexError begin
    @tensor C[a,b,c,d,e,f,g,i] := A[a,b,c,d]*B[e,f,g,h]
end

Da=10
Db=15
Dc=4
Dd=8
De=6
Df=7
Dg=3
Dh=2
A=rand(Complex128,(Da,Dc,Df,Da,De,Db,Db,Dg))
B=rand(Complex128,(Dc,Dh,Dg,De,Dd))
C=rand(Complex128,(Dd,Dh,Df))
@tensor D1[d,f,h] := A[a,c,f,a,e,b,b,g]*B[c,h,g,e,d]+0.5*C[d,h,f]
D2=zeros(Complex128,(Dd,Df,Dh))
for d=1:Dd, f=1:Df, h=1:Dh
    D2[d,f,h] += 0.5*C[d,h,f]
    for a=1:Da, b=1:Db, c=1:Dc, e=1:De, g=1:Dg
        D2[d,f,h] += A[a,c,f,a,e,b,b,g]*B[c,h,g,e,d]
    end
end
@test vecnorm(D1-D2)<eps()*sqrt(length(D1))*vecnorm(D1+D2)
@test_approx_eq(vecnorm(D1),sqrt(abs(@tensor scalar(D1[d,f,h]*conj(D1[d,f,h])))))

Abig=randn((30,30,30,30))
A=sub(Abig,1+3*(0:9),2+2*(0:6),5+4*(0:6),4+3*(0:8))
p=[3,1,4,2]
Cbig=zeros(Complex128,(50,50,50,50))
C=sub(Cbig,13+(0:6),11+4*(0:9),15+4*(0:8),4+3*(0:6))
Acopy = copy(A)
Ccopy = copy(C)
@tensor C[3,1,4,2] = A[1,2,3,4]
@tensor Ccopy[3,1,4,2] = Acopy[1,2,3,4]
@test vecnorm(C-Ccopy)<eps()*sqrt(length(C))*vecnorm(C+Ccopy)
@test_throws TensorOperations.IndexError begin
    @tensor C[3,1,4,2] = A[1,2,3]
end
@test_throws DimensionMismatch begin
    @tensor C[3,1,4,2] = A[3,1,4,2]
end
@test_throws TensorOperations.IndexError begin
    @tensor C[1,1,2,3] = A[1,2,3,4]
end

Cbig=zeros(Complex128,(50,50,50,50))
C=sub(Cbig,13+(0:6),11+4*(0:9),15+4*(0:8),4+3*(0:6))
Acopy=tensorcopy(A,1:4,p)
Ccopy=copy(C)
alpha=randn()
beta=randn()
@tensor C[3,1,4,2] = beta*C[3,1,4,2] + alpha*A[1,2,3,4]
Ccopy=beta*Ccopy+alpha*Acopy
@test vecnorm(C-Ccopy)<eps()*sqrt(length(C))*vecnorm(C+Ccopy)
@test_throws TensorOperations.IndexError begin
    @tensor C[3,1,4,2] = 0.5*C[3,1,4,2] + 1.2*A[1,2,3]
end
@test_throws DimensionMismatch  begin
    @tensor C[3,1,4,2] = 0.5*C[3,1,4,2] + 1.2*A[3,1,2,4]
end
@test_throws TensorOperations.IndexError  begin
    @tensor C[1,1,2,3] = 0.5*C[1,1,2,3] + 1.2*A[1,2,3,4]
end

Abig=rand((30,30,30,30))
A=sub(Abig,1+3*(0:8),2+2*(0:14),5+4*(0:6),7+2*(0:8))
Bbig=rand(Complex128,(50,50))
B=sub(Bbig,13+(0:14),3+5*(0:6))
Acopy=copy(A)
Bcopy=copy(B)
alpha=randn()
@tensor B[b,c] += alpha*A[a,b,c,a]
for i=1+(0:8)
    Bcopy+=alpha*slice(A,i,:,:,i)
end
@test vecnorm(B-Bcopy)<eps()*vecnorm(B+Bcopy)*sqrt(length(B))
@test_throws TensorOperations.IndexError begin
    @tensor B[b,c] += alpha*A[a,b,c]
end
@test_throws DimensionMismatch begin
    @tensor B[c,b] += alpha*A[a,b,c,a]
end
@test_throws TensorOperations.IndexError begin
    @tensor B[c,b] += alpha*A[a,b,a,a]
end
@test_throws DimensionMismatch begin
    @tensor B[c,b] += alpha*A[a,b,a,c]
end

Abig=rand((30,30,30,30))
A=sub(Abig,1+3*(0:8),2+2*(0:14),5+4*(0:6),7+2*(0:8))
Bbig=rand(Complex128,(50,50,50))
B=sub(Bbig,3+5*(0:6),7+2*(0:7),13+(0:14))
Cbig=rand(Complex64,(40,40,40))
C=sub(Cbig,3+2*(0:8),13+(0:8),7+3*(0:7))
Acopy=copy(A)
Bcopy=copy(B)
Ccopy=copy(C)
alpha=randn()
for d=1+(0:8),a=1+(0:8),e=1+(0:7)
    for b=1+(0:14),c=1+(0:6)
        Ccopy[d,a,e] -=alpha*A[a,b,c,d]*conj(B[c,e,b])
    end
end
@tensor C[d,a,e] -= alpha*A[a,b,c,d]*conj(B[c,e,b])
@test vecnorm(C-Ccopy)<eps(Float32)*vecnorm(C+Ccopy)*sqrt(length(C))
Cbig=rand(Complex64,(40,40,40))
C=sub(Cbig,3+2*(0:8),13+(0:8),7+3*(0:7))
Ccopy=copy(C)

for d=1+(0:8),a=1+(0:8),e=1+(0:7)
    for b=1+(0:14),c=1+(0:6)
        Ccopy[d,a,e] += alpha*A[a,b,c,d]*conj(B[c,e,b])
    end
end
@tensor C[d,a,e] += alpha*A[a,b,c,d]*conj(B[c,e,b])
@test vecnorm(C-Ccopy)<eps(Float32)*vecnorm(C+Ccopy)*sqrt(length(C))
@test_throws TensorOperations.IndexError begin
    @tensor C[d,a,e] += alpha*A[a,b,c,a]*B[c,e,b]
end
@test_throws TensorOperations.IndexError begin
    @tensor C[d,a,e] += alpha*A[a,b,c,d]*B[c,b]
end
@test_throws TensorOperations.IndexError begin
    @tensor C[d,e] += alpha*A[a,b,c,d]*B[c,e,b]
end
@test_throws DimensionMismatch begin
    @tensor C[d,e,a] += alpha*A[a,b,c,d]*B[c,e,b]
end

D=10;
A=randn(D,D,D,D);
B=randn(D,D,D,D);
C=randn(D,D,D,D);

@test_throws TensorOperations.IndexError begin
    @tensor C[a,b,c,d] = conj(A[c,a])*B[c,b]
end


text:
#!/usr/bin/env julia

using OrderedCollections,JSON

# tstep needs to be a divisor of
# 3600 = 2^4 ⋅ 3^2 ⋅ 5^2
tsteps = sort([2^i2*3^i3*5^i5  for i2=0:4 for i3=0:2 for i5=0:2])
filename = basename(@__FILE__)

r = OrderedCollections.OrderedDict()
r[raw"$comment"] = "File created by $filename"
r["type"] = "integer"
r["enum"]  = tsteps

io = open("tstep.schema.json","w")
JSON.print(io, r, 4)
close(io)


text:
# Engine configuration
function load_engine_config(filename::String=DEFAULT_GARAMONDRC_FILE)
    if isfile(filename)
        try
            command = "include(\"$filename\")"
            Garamond.eval(Meta.parse(command))
            @debug "Parsed and evaluated $filename."
        catch
            @warn """There was a problem evaluating $filename.
                     Assuming engine defaults."""
        end
    end
end
# Run the engine config loader
load_engine_config()
# Define the variables
SYMBOLS = Dict(:PDFTOTEXT_PROGRAM => DEFAULT_PDFTOTEXT_PROGRAM,
               :MAX_EDIT_DISTANCE => DEFAULT_MAX_EDIT_DISTANCE,
               :MAX_MATCHES => DEFAULT_MAX_MATCHES,
               :MAX_SUGGESTIONS => DEFAULT_MAX_SUGGESTIONS,
               :RESULT_AGGREGATION_STRATEGY => DEFAULT_RESULT_AGGREGATION_STRATEGY)
for (symbol, default_value) in SYMBOLS
    if !isdefined(Garamond, symbol)
        if default_value isa AbstractString
            statement = "const $symbol=\"$default_value\""
        elseif default_value isa Number
            statement = "const $symbol=$default_value"
        elseif default_value isa Symbol
            statement = "const $symbol=:$default_value"
        elseif default_value isa Type
            statement = "const $symbol=$default_value"
        else
            statement = ""  # silently ignore other types
        end
        Garamond.eval(Meta.parse(statement))
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 62f910de-3462-11ec-0510-b983317ce54a
begin
	using Pkg
	Pkg.activate("../Project.toml")
end

# ╔═╡ a2e307c3-c6c9-4689-9b16-0bd269a62c9e
begin
	using Flux, Flux.Data.MNIST, Statistics
	using Flux: throttle, params
	using Flux: @epochs, mse, onehotbatch, argmax, crossentropy
	using Base.Iterators: repeated, partition
	using Plots
	using Colors
	using Images
end

# ╔═╡ d473bd52-bc87-4c27-89ce-e54a5e63ab8a
using PlutoUI

# ╔═╡ e6a6e6e5-9b31-480e-badb-a8a65be9ad08
md"""
## References:

1) [Chapter 11, Variational Autoencoders by Hugh Murrell](https://nextjournal.com/DeepLearningNotes/Ch11Autoencoders) - trying to reproduce the above in Pluto.

2) [ACDME - Variational Autoencoder](https://docs.juliahub.com/ADCME/b8Ld2/0.7.3/vae/)

3) EECS 505 material

4) Flux Model Zoo (look for `vae_mnist` under `vision`). [Github](https://github.com/FluxML/model-zoo)
"""

# ╔═╡ 7d77c20a-6a38-4aac-aa41-443eeaf34b9d
md"""
## Display and helper functions
"""

# ╔═╡ 9d5cec86-b665-4b24-9660-57f47f0c4f13
my_norm(v) = sqrt(sum(v.^2))

# ╔═╡ a845a832-1088-4de9-992e-31d7c402e15c
function accuracy(lm, range, centroids)
    sum([ argmin(mapslices(my_norm,centroids.-lm(X[:,i]).data,dims=1))[2]-1==Y[i]
        for i in range ]) / length(range)
end

# ╔═╡ b9dfee1a-a5ac-4cb8-a7fb-6931cec23dff
function computeClasses(Y)
    sort(unique(Y))
end

# ╔═╡ 2ec27fff-e3c0-4093-8c67-2ec581a94f16
function computeCentroids(lm, data, classes)
    centroids = zeros((2,length(classes)))
    for k in 1:length(data)
        x = first(data[k])
        y = last(data[k])
        for i in 1:length(classes)
            points = lm(x[:,y.==classes[i]])
            centroids[:,i] += Flux.data.(mapslices(mean,points,dims=2))
        end
    end;
    centroids ./= length(data)
    return centroids
end

# ╔═╡ 52a6b4c8-989c-4828-9f6d-5078c5029f2d
function plotLatentSpace(lm, data, classes, width=0)
    centroids = computeCentroids(lm, data, classes)
    colors = distinguishable_colors(12, [RGB(1,1,1)])[3:12]
    if ( width > 0 )
        p = scatter(xlim=(-width,width),ylim=(-width,width))
    else
        p = scatter()
    end
    for lab = 1:length(classes)
        scatter!([centroids[1,lab]],[centroids[2,lab]],
                label=classes[lab],markersize=[9],markercolor=[colors[lab]])
    end
    for k in 1:length(data)
        x = first(data[k])
        y = last(data[k])
        for i in 1:length(classes)
            points = Flux.data(lm(x[:,y.==classes[i]]))
            scatter!(points[1,:],points[2,:],label="",
               markercolor=[colors[i]])
        end
    end
    for lab = 1:length(classes)
        scatter!([centroids[1,lab]],[centroids[2,lab]],
                label="",markersize=[9],markercolor=[colors[lab]])
    end
    scatter!(centroids[1,:],centroids[2,:],label="centroids")
    return p
end

# ╔═╡ c25f4c25-959d-4077-9ce5-fa0b26641b17
md"""
## Simple Autoencoder
"""

# ╔═╡ 3db0aa21-7a71-4d21-bc7f-f236b6dde600
# begin
# 	# Latent dimensionality, # hidden units.
# 	Dz, Dh = 2, 500

# 	# encoder
# 	g = Chain(Dense(28^2, Dh, tanh), Dense(Dh, Dz))
# 	# decoder
# 	f = Chain(Dense(Dz, Dh, tanh), Dense(Dh, 28^2, σ))
# 	# model
# 	sae = Chain(g,f) |> gpu
# 	# loss
# 	loss(X,Y) = mse(sae(X), X)
# 	# callback
# 	evalcb = throttle(() -> ( p=rand(1:N, M); @show(loss(Float32.(X[:, p]),Y[p]))), 20)
# 	# optimization
# 	opt = ADAM()
# 	# parameters
# 	ps = params(sae);
# end

# ╔═╡ 9b6d280f-0994-4caf-8a07-bdf7f03cfda9
# p=rand(1:N, M)

# ╔═╡ 122fd414-1bbc-45cc-988a-bf4dda9ee5ae
# loss(Float32.(X[:,p]), Y[p])

# ╔═╡ 04ce4d2a-8b37-4e33-af12-15e2f05d9fb5
# trainDataGPU = gpu.(trainData);

# ╔═╡ 72b3993c-55fb-4a34-bf6b-87e538ab8b4c
# trainDataGPU

# ╔═╡ b78786c3-df99-46e6-a574-3c9a101f2bb6
# @epochs 2 Flux.train!(loss, params(sae), trainDataGPU, opt, cb = evalcb)

# ╔═╡ Cell order:
# ╟─e6a6e6e5-9b31-480e-badb-a8a65be9ad08
# ╠═62f910de-3462-11ec-0510-b983317ce54a
# ╠═a2e307c3-c6c9-4689-9b16-0bd269a62c9e
# ╠═d473bd52-bc87-4c27-89ce-e54a5e63ab8a
# ╠═7d77c20a-6a38-4aac-aa41-443eeaf34b9d
# ╠═9d5cec86-b665-4b24-9660-57f47f0c4f13
# ╠═a845a832-1088-4de9-992e-31d7c402e15c
# ╠═b9dfee1a-a5ac-4cb8-a7fb-6931cec23dff
# ╠═2ec27fff-e3c0-4093-8c67-2ec581a94f16
# ╠═52a6b4c8-989c-4828-9f6d-5078c5029f2d
# ╠═c25f4c25-959d-4077-9ce5-fa0b26641b17
# ╠═3db0aa21-7a71-4d21-bc7f-f236b6dde600
# ╠═9b6d280f-0994-4caf-8a07-bdf7f03cfda9
# ╠═122fd414-1bbc-45cc-988a-bf4dda9ee5ae
# ╠═04ce4d2a-8b37-4e33-af12-15e2f05d9fb5
# ╠═72b3993c-55fb-4a34-bf6b-87e538ab8b4c
# ╠═b78786c3-df99-46e6-a574-3c9a101f2bb6

max_stars_repo_name: aniketjivani/variational-autoencoders

max_issues_repo_name: aniketjivani/variational-autoencoders

max_forks_repo_name: aniketjivani/variational-autoencoders


text:
"""
	vec2mat(x::Vector{Vector{Float64}})
	vec2mat(x::Vector{Float64})

Convert vector of vectors into matrix. If input is a single vector do nothing.
"""
vec2mat(x::Vector{Vector{Float64}}) = reduce(hcat, x)
vec2mat(x::Vector{Float64}) = x

"""
    vec2mat!(Σ::ControlAffineSystem)

In-place modification of system trajectory
"""
function vec2mat!(Σ::ControlAffineSystem)
    Σ.xs = vec2mat(Σ.xs)

    return Σ
end