text:
"""
    diff_optimizer(optimizer_constructor)::Optimizer

Creates a `DiffOpt.Optimizer`, which is an MOI layer with an internal optimizer
and other utility methods. Results (primal, dual and slack values) are obtained
by querying the internal optimizer instantiated using the
`optimizer_constructor`. These values are required for find jacobians with respect to problem data.

One define a differentiable model by using any solver of choice. Example:

```julia
julia> import DiffOpt, HiGHS

julia> model = DiffOpt.diff_optimizer(HiGHS.Optimizer)
julia> model.add_variable(x)
julia> model.add_constraint(...)

julia> _backward_quad(model)  # for convex quadratic models

julia> _backward_quad(model)  # for convex conic models
```
"""
function diff_optimizer(optimizer_constructor)::Optimizer
    optimizer = MOI.instantiate(optimizer_constructor, with_bridge_type=Float64)
    # When we do `MOI.copy_to(diff, optimizer)` we need to efficiently `MOI.get`
    # the model information from `optimizer`. However, 1) `optimizer` may not
    # implement some getters or it may be inefficient and 2) the getters may be
    # unimplemented or inefficient through some bridges.
    # For this reason we add a cache layer, the same cache JuMP adds.
    caching_opt = MOIU.CachingOptimizer(
        MOIU.UniversalFallback(MOIU.Model{Float64}()),
        optimizer,
    )
    return Optimizer(caching_opt)
end

mutable struct Optimizer{OT <: MOI.ModelLike} <: MOI.AbstractOptimizer
    optimizer::OT

    program_class::ProgramClassCode

    diff::Union{Nothing,MOI.Bridges.LazyBridgeOptimizer{QPDiff},MOI.Bridges.LazyBridgeOptimizer{ConicDiff}}

    index_map::Union{Nothing,MOI.Utilities.IndexMap}

    # sensitivity input cache using MOI like sparse format
    input_cache::DiffInputCache

    function Optimizer(optimizer::OT) where {OT <: MOI.ModelLike}
        new{OT}(
            optimizer,
            AUTOMATIC,
            nothing,
            nothing,
            DiffInputCache(),
        )
    end
end

function MOI.add_variable(model::Optimizer)
    model.diff = nothing
    vi = MOI.add_variable(model.optimizer)
    return vi
end

function MOI.add_variables(model::Optimizer, N::Int)
    model.diff = nothing
    return VI[MOI.add_variable(model) for i in 1:N]
end

function MOI.add_constraint(model::Optimizer, f::SUPPORTED_SCALAR_FUNCTIONS, s::SUPPORTED_SCALAR_SETS)
    model.diff = nothing
    return MOI.add_constraint(model.optimizer, f, s)
end

function MOI.add_constraint(model::Optimizer, vf::SUPPORTED_VECTOR_FUNCTIONS, s::SUPPORTED_VECTOR_SETS)
    model.diff = nothing
    return MOI.add_constraint(model.optimizer, vf, s)
end

function MOI.add_constraints(model::Optimizer, f::AbstractVector{F}, s::AbstractVector{S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.diff = nothing
    return CI{F, S}[MOI.add_constraint(model, f[i], s[i]) for i in eachindex(f)]
end

function MOI.set(model::Optimizer, attr::MOI.ObjectiveFunction{<: SUPPORTED_OBJECTIVES}, f::SUPPORTED_OBJECTIVES)
    model.diff = nothing
    MOI.set(model.optimizer, attr, f)
end

MOI.supports(model::Optimizer, attr::MOI.ObjectiveSense) = MOI.supports(model.optimizer, attr)

function MOI.set(model::Optimizer, attr::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.diff = nothing
    return MOI.set(model.optimizer, attr, sense)
end

function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute)
    return MOI.get(model.optimizer, attr)
end

function MOI.get(
    model::Optimizer, attr::MOI.ListOfConstraintIndices{F, S}
) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(
    model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}
) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(
    model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S
) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    model.diff = nothing
    return MOI.set(model.optimizer, attr, ci, s)
end

function MOI.get(
    model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}
) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(
    model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F
) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    model.diff = nothing
    return MOI.set(model.optimizer, attr, ci, f)
end

function MOI.get(
    model::Optimizer, attr::MOI.ListOfConstraintIndices{F, S}
) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(
    model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}
) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(
    model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S
) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    model.diff = nothing
    return MOI.set(model.optimizer, attr, ci, s)
end

function MOI.get(
    model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}
) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(
    model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F
) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    model.diff = nothing
    return MOI.set(model.optimizer, attr, ci, f)
end

# `MOI.supports` methods

function MOI.supports(model::Optimizer, attr::MOI.AbstractModelAttribute)
    return MOI.supports(model.optimizer, attr)
end

function MOI.supports(
    model::Optimizer,
    attr::MOI.ObjectiveFunction{<: SUPPORTED_OBJECTIVES}
)
    return MOI.supports(model.optimizer, attr)
end

function MOI.supports_constraint(
    model::Optimizer,
    ::Type{F},
    ::Type{S}
) where {
    F <: SUPPORTED_SCALAR_FUNCTIONS,
    S <: SUPPORTED_SCALAR_SETS,
}
    MOI.supports_constraint(model.optimizer, F, S)
end
function MOI.supports_constraint(
    model::Optimizer,
    ::Type{F},
    ::Type{S}
) where {
    F <: SUPPORTED_VECTOR_FUNCTIONS,
    S <: SUPPORTED_VECTOR_SETS,
}
    MOI.supports_constraint(model.optimizer, F, S)
end
function MOI.supports(
    model::Optimizer, attr::MOI.ConstraintName, ::Type{CI{F, S}}
) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.supports(model.optimizer, attr, CI{F,S})
end
function MOI.supports(
    model::Optimizer, attr::MOI.ConstraintName, ::Type{CI{F, S}}
) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.supports(model.optimizer, attr, CI{F,S})
end

MOI.get(model::Optimizer, attr::MOI.SolveTimeSec) = MOI.get(model.optimizer, attr)

function MOI.empty!(model::Optimizer)
    MOI.empty!(model.optimizer)
    model.diff = nothing
    model.index_map = nothing
    empty!(model.input_cache)
    return
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.optimizer) &&
           model.diff === nothing
end

function MOI.supports_incremental_interface(model::Optimizer)
    if !MOI.supports_incremental_interface(model.optimizer)
        error(
            "DiffOpt requires a solver that " *
            "`MOI.supports_incremental_interface`, which is not the case for " *
            "$(model.optimizer)")
    end
    return true
end

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike)
    model.diff = nothing
    return MOIU.default_copy_to(model.optimizer, src)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    return MOI.get(model.optimizer, MOI.TerminationStatus())
end

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::VI, value::Float64)
    model.diff = nothing
    MOI.set(model.optimizer, MOI.VariablePrimalStart(), vi, value)
end

function MOI.supports(model::Optimizer, attr::MOI.AbstractVariableAttribute,
                      ::Type{MOI.VariableIndex})
    return MOI.supports(model.optimizer, attr, MOI.VariableIndex)
end

function MOI.set(model::Optimizer, attr::MOI.AbstractVariableAttribute, v::MOI.VariableIndex, value)
    MOI.set(model.optimizer, attr, v, value)
end

function MOI.get(model::Optimizer, attr::MOI.AbstractVariableAttribute, vi::MOI.VariableIndex)
    return MOI.get(model.optimizer, attr, vi)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.diff = nothing
    MOI.delete(model.optimizer, ci)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    model.diff = nothing
    MOI.delete(model.optimizer, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.is_valid(model::Optimizer, v::VI)
    return MOI.is_valid(model.optimizer, v::VI)
end

MOI.is_valid(model::Optimizer, con::CI) = MOI.is_valid(model.optimizer, con)

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintBasisStatus, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintBasisStatus(), ci)
end

# helper methods to check if a constraint contains a Variable
function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VariableIndex, S}) where {S <: SUPPORTED_SCALAR_SETS}
    return v == MOI.get(model, MOI.ConstraintFunction(), ci)
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.ScalarAffineFunction{Float64}, S}) where {S <: SUPPORTED_SCALAR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.variable, func.terms)
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorOfVariables, S}) where {S <: SUPPORTED_VECTOR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return v in func.variables
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorAffineFunction{Float64}, S}) where {S <: SUPPORTED_VECTOR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.scalar_term.variable, func.terms)
end

function MOI.delete(model::Optimizer, v::VI)
    model.diff = nothing
    MOI.delete(model.optimizer, v)
end

# for array deletion
function MOI.delete(model::Optimizer, indices::Vector{VI})
    model.diff = nothing
    for i in indices
        MOI.delete(model, i)
    end
end

function MOI.modify(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    chg::MOI.AbstractFunctionModification
)
    model.diff = nothing
    MOI.modify(
        model.optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        chg
    )
end

function MOI.modify(model::Optimizer, ci::CI{F, S}, chg::MOI.AbstractFunctionModification) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.diff = nothing
    MOI.modify(model.optimizer, ci, chg)
end

function MOI.modify(model::Optimizer, ci::CI{F, S}, chg::MOI.AbstractFunctionModification) where {F<:MOI.VectorAffineFunction{Float64}, S <: SUPPORTED_VECTOR_SETS}
    model.diff = nothing
    MOI.modify(model.optimizer, ci, chg)
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    return MOI.get(model.optimizer, MOI.VariableIndex, name)
end

function MOI.set(model::Optimizer, ::MOI.ConstraintName, con::CI, name::String)
    MOI.set(model.optimizer, MOI.ConstraintName(), con, name)
end

function MOI.get(model::Optimizer, ::Type{MOI.ConstraintIndex}, name::String)
    return MOI.get(model.optimizer, MOI.ConstraintIndex, name)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, con::CI)
    return MOI.get(model.optimizer, MOI.ConstraintName(), con)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{F,S})
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{VF, VS}}) where {VF<:SUPPORTED_VECTOR_FUNCTIONS, VS<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{VF,VS})
end

function MOI.set(model::Optimizer, ::MOI.Name, name::String)
    MOI.set(model.optimizer, MOI.Name(), name)
end

function MOI.get(model::Optimizer, ::Type{CI{F, S}}, name::String) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, CI{F,S}, name)
end

function MOI.get(model::Optimizer, ::Type{CI{VF, VS}}, name::String) where {VF<:SUPPORTED_VECTOR_FUNCTIONS, VS<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, CI{VF,VS}, name)
end

function MOI.supports(model::Optimizer, attr::MOI.TimeLimitSec)
    MOI.supports(model.optimizer, attr)
end
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Union{Real, Nothing})
    MOI.set(model.optimizer, MOI.TimeLimitSec(), value)
end
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(model.optimizer, MOI.TimeLimitSec())
end

MOI.supports(model::Optimizer, ::MOI.Silent) = MOI.supports(model.optimizer, MOI.Silent())

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    MOI.set(model.optimizer, MOI.Silent(), value)
end

function MOI.get(model::Optimizer, ::MOI.Silent)
    return MOI.get(model.optimizer, MOI.Silent())
end

function MOI.get(model::Optimizer, ::MOI.SolverName)
    return MOI.get(model.optimizer, MOI.SolverName())
end

function MOI.optimize!(model::Optimizer)
    model.diff = nothing
    MOI.optimize!(model.optimizer)

    # do not fail. interferes with MOI.Tests.linear12test
    if !in(MOI.get(model.optimizer, MOI.TerminationStatus()),  (MOI.LOCALLY_SOLVED, MOI.OPTIMAL))
        @warn "problem status: $(MOI.get(model.optimizer, MOI.TerminationStatus()))"
        return
    end

    return
end

"""
    ProgramClass <: MOI.AbstractOptimizerAttribute

Determines which program class to used from [`ProgramClassCode`](@ref). The
default is `AUTOMATIC`.

One important advantage of setting the class explicitly is that it will allow
necessary bridges to be used. If the class is `AUTOMATIC` then
`DiffOpt.Optimizer` will report that it supports both objective and constraints
of the QP and CP classes. For instance, it will reports that is supports both
quadratic objective and conic constraints. However, at the differentiation
stage, we won't be able to differentiate since QP does not support conic
constraints and CP does not support quadratic objective. On the other hand, if
the `ProgramClass` is set to `CONIC` then `DiffOpt.Optimizer` will report that
it does not support quadratic objective hence it will be bridged to second-order
cone constraints and we will be able to use CP to differentiate.
"""
struct ProgramClass <: MOI.AbstractOptimizerAttribute end

MOI.supports(::Optimizer, ::ProgramClass) = true
MOI.get(model::Optimizer, ::ProgramClass) = model.program_class
function MOI.set(model::Optimizer, ::ProgramClass, class::ProgramClassCode)
    model.program_class = class
end

"""
    ProgramClassUsed <: MOI.AbstractOptimizerAttribute

Program class actually used, same as [`ProgramClass`](@ref) except that it does
not return `AUTOMATIC` but the class automatically chosen instead. This attribute
is read-only, it cannot be set, set [`ProgramClass`](@ref) instead.
"""
struct ProgramClassUsed <: MOI.AbstractOptimizerAttribute end

function MOI.get(model::Optimizer, ::ProgramClassUsed)
    if model.program_class == AUTOMATIC
        if _qp_supported(model.optimizer)
            return QUADRATIC
        else
            return CONIC
        end
    else
        return model.program_class
    end
end

"""
    backward(model::Optimizer)

Wrapper method for the backward pass.
This method will consider as input a currently solved problem and differentials
with respect to the solution set with the [`BackwardInVariablePrimal`](@ref) attribute.
The output problem data differentials can be queried with the
attributes [`BackwardOutObjective`](@ref) and [`BackwardOutConstraint`](@ref).
"""
function backward(model::Optimizer)
    diff = _diff(model)
    for (vi, value) in model.input_cache.dx
        MOI.set(diff, BackwardInVariablePrimal(), model.index_map[vi], value)
    end
    backward(diff)
end

function _copy_forward_in_constraint(diff, index_map, con_map, constraints)
    for (index, value) in constraints
        MOI.set(diff, ForwardInConstraint(), con_map[index], MOI.Utilities.map_indices(index_map, value))
    end
end

"""
    forward(model::Optimizer)

Wrapper method for the forward pass.
This method will consider as input a currently solved problem and
differentials with respect to problem data set with
the [`ForwardInObjective`](@ref) and  [`ForwardInConstraint`](@ref) attributes.
The output solution differentials can be queried with the attribute
[`ForwardOutVariablePrimal`](@ref).
"""
function forward(model::Optimizer)
    diff = _diff(model)
    if model.input_cache.objective !== nothing
        MOI.set(diff, ForwardInObjective(), MOI.Utilities.map_indices(model.index_map, model.input_cache.objective))
    end
    for (F, S) in keys(model.input_cache.scalar_constraints.dict)
        _copy_forward_in_constraint(diff, model.index_map, model.index_map.con_map[F, S], model.input_cache.scalar_constraints[F, S])
    end
    for (F, S) in keys(model.input_cache.vector_constraints.dict)
        _copy_forward_in_constraint(diff, model.index_map, model.index_map.con_map[F, S], model.input_cache.vector_constraints[F, S])
    end
    forward(diff)
end

function _copy_constraint_start(dest, src, index_map::MOIU.DoubleDicts.IndexDoubleDictInner{F, S}, dest_attr, src_attr) where {F, S}
    for ci in MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
        value = MOI.get(src, src_attr, ci)
        MOI.set(dest, dest_attr, index_map[ci], value)
    end
end

function _diff(model::Optimizer)
    if model.diff === nothing
        if MOI.get(model, ProgramClassUsed()) == QUADRATIC
            diff = QPDiff()
        else
            _check_termination_status(model)
            diff = ConicDiff()
        end
        model.diff = MOI.Bridges.LazyBridgeOptimizer(diff)
        # We don't add any variable bridge here because:
        # 1) If `ZerosBridge` is used, `MOI.Bridges.unbridged_function` does not work.
        #    This is in fact expected: since `ZerosBridge` drops the variable, we dont
        #    compute the derivative of the value of this variable as a function of its fixed value.
        #    This could be easily determined as the same as the derivative of the value but
        #    since the variable was also dropped from other constraints, we would ignore its impact on the other constraints.
        # 2) For affine variable bridges, `bridged_function` and `unbridged_function` don't treat the function as a derivative hence they will add constants
        MOI.Bridges.Constraint.add_all_bridges(model.diff, Float64)
        MOI.Bridges.Objective.add_all_bridges(model.diff, Float64)
        model.index_map = MOI.copy_to(model.diff, model.optimizer)
        vis_src = MOI.get(model.optimizer, MOI.ListOfVariableIndices())
        MOI.set(model.diff, MOI.VariablePrimalStart(), getindex.(Ref(model.index_map), vis_src), MOI.get(model.optimizer, MOI.VariablePrimal(), vis_src))
        for (F, S) in MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())
            _copy_constraint_start(model.diff, model.optimizer, model.index_map.con_map[F, S], MOI.ConstraintPrimalStart(), MOI.ConstraintPrimal())
            _copy_constraint_start(model.diff, model.optimizer, model.index_map.con_map[F, S], MOI.ConstraintDualStart(), MOI.ConstraintDual())
        end
    end
    return model.diff
end

function _check_termination_status(model::Optimizer)
    if !in(
        MOI.get(model, MOI.TerminationStatus()), (MOI.LOCALLY_SOLVED, MOI.OPTIMAL)
        )
        error("problem status: ", MOI.get(model.optimizer, MOI.TerminationStatus()))
    end
end

# DiffOpt attributes redirected to `diff`

function _checked_diff(model::Optimizer, attr::MOI.AnyAttribute, call)
    if model.diff === nothing
        error("Cannot get attribute `attr`. First call `DiffOpt.$call`.")
    end
    return model.diff
end

function MOI.get(model::Optimizer, attr::BackwardOutObjective)
    return IndexMappedFunction(
        MOI.get(_checked_diff(model, attr, :backward), attr),
        model.index_map,
    )
end
function MOI.get(model::Optimizer, ::ForwardInObjective)
    return model.input_cache.objective
end
function MOI.set(model::Optimizer, ::ForwardInObjective, objective)
    model.input_cache.objective = objective
    return
end

function MOI.get(model::Optimizer, attr::ForwardOutVariablePrimal, vi::MOI.VariableIndex)
    return MOI.get(_checked_diff(model, attr, :forward), attr, model.index_map[vi])
end
function MOI.get(model::Optimizer, ::BackwardInVariablePrimal, vi::VI)
    return get(model.input_cache.dx, vi, 0.0)
end
function MOI.set(model::Optimizer, ::BackwardInVariablePrimal, vi::VI, val)
    model.input_cache.dx[vi] = val
    return
end

function MOI.get(model::Optimizer, attr::BackwardOutConstraint, ci::MOI.ConstraintIndex)
    return IndexMappedFunction(
        MOI.get(_checked_diff(model, attr, :backward), attr, model.index_map[ci]),
        model.index_map,
    )
end
function MOI.get(model::Optimizer,
    ::ForwardInConstraint, ci::CI{MOI.ScalarAffineFunction{T},S}
) where {T,S}
    return get(model.input_cache.scalar_constraints, ci, zero(MOI.ScalarAffineFunction{T}))
end
function MOI.get(model::Optimizer,
    ::ForwardInConstraint, ci::CI{MOI.VectorAffineFunction{T},S}
) where {T,S}
    func = get(model.input_cache.vector_constraints, ci, nothing)
    if func === nothing
        set = MOI.get(model, MOI.ConstraintSet(), ci)
        dim = MOI.dimension(set)
        return MOI.Utilities.zero_with_output_dimension(MOI.VectorAffineFunction{T}, dim)
    else
        return func
    end
end
function MOI.set(model::Optimizer,
    ::ForwardInConstraint,
    ci::CI{MOI.ScalarAffineFunction{T},S},
    func::MOI.ScalarAffineFunction{T},
) where {T,S}
    model.input_cache.scalar_constraints[ci] = func
    return
end
function MOI.set(model::Optimizer,
    ::ForwardInConstraint,
    ci::CI{MOI.VectorAffineFunction{T},S},
    func::MOI.VectorAffineFunction{T},
) where {T,S}
    model.input_cache.vector_constraints[ci] = func
    return
end


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

mktempdir() do dir

  # Create test file
  filename = joinpath(dir, "file.txt")
  text = "123456"
  write(filename, text)

  # test filesystem truncate (shorten)
  file = Base.Filesystem.open(filename, Base.Filesystem.JL_O_RDWR)
  Base.Filesystem.truncate(file, 2)
  @test length(read(file)) == 2
  close(file)

  # test filesystem truncate (lengthen)
  file = Base.Filesystem.open(filename, Base.Filesystem.JL_O_RDWR)
  Base.Filesystem.truncate(file, 20)
  @test length(read(file)) == 20
  close(file)

  # test filesystem futime
  file = Base.Filesystem.open(filename, Base.Filesystem.JL_O_RDWR)
  Base.Filesystem.futime(file, 1.0, 2.0)
  @test Base.Filesystem.stat(file).mtime == 2.0
  close(file)

  # test filesystem readbytes!
  file = Base.Filesystem.open(filename, Base.Filesystem.JL_O_RDWR)
  res = ones(UInt8, 80)
  Base.Filesystem.readbytes!(file, res)
  @test res == UInt8[text..., (i > 20 for i in (length(text) + 1):length(res))...]
  close(file)

end

ext: jl

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using Literate, MIToS

COOKBOOK = joinpath(@__DIR__, "src", "cookbook")
MD_OUTPUT = joinpath(@__DIR__, "src")
NB_OUTPUT = joinpath(@__DIR__, "src", "cookbook", "notebooks")

for file in [ "01_Change_B_factors.jl",
              "02_Linking_structural_and_evolutionary_information.jl"]
    Literate.markdown(joinpath(COOKBOOK, file),
        MD_OUTPUT, execute=false, documenter=true)
    Literate.notebook(joinpath(COOKBOOK, file),
        NB_OUTPUT, execute=false, documenter=true)
end


text:
function second_order_filter(f, q, fs)
    cft = f ./ fs
    rho = exp.(-pi .* cft ./ q)
    theta = 2pi .* cft .* sqrt.(1 .- 1.0 ./ (4 .* q .^ 2))
    filts = hcat(ones(size(rho)), -2rho .* cos.(theta), rho .^ 2)
end

function freq_resp(filter, f, fs)
    cf = exp(1im * 2 * pi * f / fs)
    denom = filter[5] + filter[4] * cf + cf ^ 2
    mag = (filter[3] + filter[2] * cf + filter[1] * cf ^ 2) / denom
    mag = 20 * log10(abs(mag))
end

function set_gain(filter, desired, f, fs)
    old_gain = 10 ^ (freq_resp(filter, f, fs) / 20)
    filter[1:3] *= desired / old_gain
    filter
end

function design_lyon_filters(fs, ear_q = 8, step_factor = ear_q / 32)
    Eb = 1000.0
    EarZeroOffset = 1.5
    EarSharpness = 5.0
    EarPremphCorner = 300

    # Find top frequency, allowing space for first cascade filter.
    topf = fs / 2.0
    topf_neg_delta = sqrt(topf ^ 2 + Eb ^ 2) / ear_q * step_factor * EarZeroOffset
    topf = topf - topf_neg_delta + sqrt(topf ^ 2 + Eb ^ 2) / ear_q * step_factor

    # Find place where CascadePoleQ < .5
    lowf = Eb / sqrt(4 * ear_q ^ 2 - 1)
    _log_val_low = log(lowf + sqrt(lowf ^ 2 + Eb ^ 2))
    _log_val_top = log(topf + sqrt(Eb ^ 2 + topf ^ 2))
    NumberOfChannels = floor(Int, (ear_q * (-_log_val_low + _log_val_top)) / step_factor)

    # Now make an array of CenterFreqs..... This expression was derived by
    # Mathematica by integrating 1/EarBandwidth(cf) and solving for f as a
    # function of channel number.
    cn = 1:NumberOfChannels
    denom = exp.((cn .* step_factor) ./ ear_q)
    nom_add = topf + sqrt(Eb ^ 2 + topf ^ 2)
    center_freqs = (-1 .* ((exp.((cn .* step_factor) ./ ear_q) .* Eb ^ 2) ./ nom_add) .+ nom_add ./ denom) ./ 2.0

    # OK, now we can figure out the parameters of each stage filter.
    EarBandwidth = sqrt.(center_freqs .^ 2 .+ Eb ^ 2) ./ ear_q
    CascadeZeroCF = center_freqs .+ EarBandwidth .* step_factor .* EarZeroOffset
    CascadeZeroQ = EarSharpness .* CascadeZeroCF ./ EarBandwidth
    CascadePoleCF = center_freqs
    CascadePoleQ = center_freqs ./ EarBandwidth

    # Now lets find some filters.... first the zeros then the poles
    zerofilts = second_order_filter(CascadeZeroCF, CascadeZeroQ, fs)
    polefilts = second_order_filter(CascadePoleCF, CascadePoleQ, fs)
    # filters = np.vstack([zerofilts, polefilts[1:, :]])
    filters = cat(zerofilts, polefilts[:, 2:end], dims = 2)

    # Now we can set the DC gain of each stage.
    dcgain = zeros(NumberOfChannels)

    dcgain[2:end] .= center_freqs[1:end-1] ./ center_freqs[2:end]
    dcgain[1] = dcgain[2]

    for i in 1:NumberOfChannels
        filters[i, :] = set_gain(filters[i, :], dcgain[i], 0, fs)
    end

    # Finally, let's design the front filters.
    front = zeros(2, 5)
    front_0 = [0, 1, -exp(-2 * pi * EarPremphCorner / fs), 0, 0]
    front[1, :] = set_gain(front_0, 1, fs / 4, fs)
    top_poles = second_order_filter(topf, CascadePoleQ[1], fs)
    front_1 = [1, 0, -1, top_poles[2], top_poles[3]]
    front[2, :] = set_gain(front_1, 1, fs / 4, fs)

    # Now, put them all together.
    filters = cat(front, filters, dims = 1)
    return filters, center_freqs
end

function epsilon_from_tau(tau, sample_rate)
    1 - exp(-1.0 / tau / sample_rate)
end


text:
function build(diags::AbstractVector, hasLoop = true; verbose::Int = 0)
    @assert all(d -> (d.id.para == diags[1].id.para), diags) "Parameters of all diagrams shoud be the same!"

    diags = DiagTree.optimize(diags, verbose = verbose)

    if hasLoop
        tree = newExprTree(diags[1].id.para, :none)
    else
        tree = newExprTreeXT(diags[1].id.para, :none)
    end

    verbose > 0 && println("Constructing expression tree...")
    nodes = Dict{Int,Any}()
    for diag in diags
        for d in PostOrderDFS(diag)
            if haskey(nodes, d.hash) == false
                id = d.id
                if isempty(d.subdiagram)
                    K = hasLoop ? id.extK : nothing
                    nodes[d.hash] = addpropagator!(tree, d.name, d.factor; site = collect(id.extT), loop = K, para = id)
                else
                    children = [nodes[sub.hash] for sub in d.subdiagram]
                    nodes[d.hash] = addnode!(tree, operator(d.operator), d.name, children, d.factor, para = id)
                end
            end
        end
    end

    tree.root = [nodes[d.hash] for d in diags]

    initialize!(tree.node)
    return tree
end
function build(diag::Diagram{W}; verbose::Int = 0) where {W}
    diag = build([diag,], verbose = verbose)
    return diag
end

function operator(op::Operator)
    if op isa Sum
        return ADD
    elseif op isa Prod
        return MUL
    else
        error("$op not implemented!")
    end
end

function newExprTree(para, name::Symbol = :none)
    weightType = para.weightType
    Kpool = LoopPool(:K, para.loopDim, para.totalLoopNum, Float64)
    return ExpressionTree(loopBasis = Kpool, nodePara = DiagramId, weight = weightType, name = name)
end

function newExprTreeXT(para, name::Symbol = :none)
    weightType = para.weightType
    return ExpressionTree(loopBasis = nothing, nodePara = DiagramId, weight = weightType, name = name)
end


text:
using Dates: Day, Hour, Minute

utc = FixedTimeZone("UTC", 0)
dst = first(compile("America/Winnipeg", tzdata["northamerica"]))
no_dst = first(compile("America/Regina", tzdata["northamerica"]))
fixed = FixedTimeZone("Fixed", -5 * 3600)

# Test date ranges for dates in which no DST transition occurs.
d = DateTime(2013, 2, 13)
raw_range = collect(d:Hour(1):d + Day(1))

utc_range = ZonedDateTime(d, utc):Hour(1):ZonedDateTime(d, utc) + Day(1)
dst_range = ZonedDateTime(d, dst):Hour(1):ZonedDateTime(d, dst) + Day(1)
no_dst_range = ZonedDateTime(d, no_dst):Hour(1):ZonedDateTime(d, no_dst) + Day(1)
fixed_range = ZonedDateTime(d, fixed):Hour(1):ZonedDateTime(d, fixed) + Day(1)

# Each range should have 25 elements.
@test length(utc_range) == 25
@test length(dst_range) == 25
@test length(no_dst_range) == 25
@test length(fixed_range) == 25

# Values in ZoneDateTime range should match raw DateTime range with appropriate offsets.
@test DateTime.(utc_range, UTC) == raw_range
@test DateTime.(dst_range, UTC) == raw_range .+ Hour(6)
@test DateTime.(no_dst_range, UTC) == raw_range .+ Hour(6)
@test DateTime.(fixed_range, UTC) == raw_range .+ Hour(5)

# Test date ranges for the day of the DST "spring forward".
d = DateTime(2013, 3, 10)
raw_range = collect(d:Hour(1):d + Day(1))

utc_range = ZonedDateTime(d, utc):Hour(1):ZonedDateTime(d, utc) + Day(1)
dst_range = ZonedDateTime(d, dst):Hour(1):ZonedDateTime(d, dst) + Day(1)
no_dst_range = ZonedDateTime(d, no_dst):Hour(1):ZonedDateTime(d, no_dst) + Day(1)
fixed_range = ZonedDateTime(d, fixed):Hour(1):ZonedDateTime(d, fixed) + Day(1)

# The range that observes DST should have only 24 elements because 02:00 is skipped.
@test length(utc_range) == 25
@test length(dst_range) == 24
@test length(no_dst_range) == 25
@test length(fixed_range) == 25

# Values in ZoneDateTime range should match raw DateTime range with appropriate offsets.
@test DateTime.(utc_range, UTC) == raw_range
@test DateTime.(dst_range, UTC) == [
    raw_range[1:2] .+ Hour(6);
    raw_range[4:end] .+ Hour(5)
]
@test DateTime.(no_dst_range, UTC) == raw_range .+ Hour(6)
@test DateTime.(fixed_range, UTC) == raw_range .+ Hour(5)

# Test date ranges for the day of the DST "fall back".
d = DateTime(2013, 11, 3)
raw_range = collect(d:Hour(1):d + Day(1))

utc_range = ZonedDateTime(d, utc):Hour(1):ZonedDateTime(d, utc) + Day(1)
dst_range = ZonedDateTime(d, dst):Hour(1):ZonedDateTime(d, dst) + Day(1)
no_dst_range = ZonedDateTime(d, no_dst):Hour(1):ZonedDateTime(d, no_dst) + Day(1)
fixed_range = ZonedDateTime(d, fixed):Hour(1):ZonedDateTime(d, fixed) + Day(1)

# The range that observes DST should have 26 elements because of the two instances of 01:00.
@test length(utc_range) == 25
@test length(dst_range) == 26
@test length(no_dst_range) == 25
@test length(fixed_range) == 25

# Values in ZoneDateTime range should match raw DateTime range with appropriate offsets.
@test DateTime.(utc_range, UTC) == raw_range
@test DateTime.(dst_range, UTC) == [
    raw_range[1:2] .+ Hour(5);
    raw_range[2:end] .+ Hour(6);
]
@test DateTime.(no_dst_range, UTC) == raw_range .+ Hour(6)
@test DateTime.(fixed_range, UTC) == raw_range .+ Hour(5)

# filter behaviour with a non-existent hour
range = ZonedDateTime(2015, 3, 8, dst):Dates.Hour(1):ZonedDateTime(2015, 3, 10, dst)
@test filter(dt -> Dates.hour(dt) == 2, range) == [
    # ZonedDateTime(2015, 3, 8, 2, dst)  # Non-existent hour
    ZonedDateTime(2015, 3, 9, 2, dst)
]

# filter behaviour with ambiguous hour
range = ZonedDateTime(2015, 11, 1, dst):Dates.Hour(1):ZonedDateTime(2015, 11, 3, dst)
@test filter(dt -> Dates.hour(dt) == 1, range) == [
    ZonedDateTime(2015, 11, 1, 1, dst, 1)
    ZonedDateTime(2015, 11, 1, 1, dst, 2)
    ZonedDateTime(2015, 11, 2, 1, dst)
]

# default step for a ZonedDateTime range
# Currently unable to test the deprecation
# @test_warn(
#     "colon(start::T, stop::T) where T <: ZonedDateTime is deprecated, use start:Day(1):stop instead.",
#     ZonedDateTime(2017, 10, 1, 9, utc):ZonedDateTime(2017, 12, 8, 23, utc)
# )


text:
#
#--------------------------------------#
export hlmz
#--------------------------------------#
"""
 for v,u in H^1 of Omega

 (v,(-ν∇² + k)u)\n
          = ν * a(v,u)
          + k *  (v,u)
"""
function hlmz(u::Array
             ,ν,k,msh::Mesh)

    Hu   = ν .* lapl(u,msh)
    Hu .+= k .* mass(u,msh)

return Hu
end
#--------------------------------------#
# dealiased version
function hlmz(u::Array
             ,ν,k
             ,msh1::Mesh
             ,msh2::Mesh)

    Hu = hlmz(u,ν,k,msh1)

return Hu
end
#--------------------------------------#
#


text:
# SD_Wave Functions for Keysight M32XXA/M33XXA AWGs
export SD_Wave_newFromArrayInteger, SD_Wave_newFromArrayDouble,
    SD_Wave_newFromFile, SD_Wave_delete

## new
"""
This function creates a waveform object from data points contained in an array
in memory.
### AWG Waveform types
- Analog 16 Bits : Analog normalized waveforms (-1..1) defined with doubles
--> WAVE_ANALOG_16 = 0
- Analog 32 Bits : Analog normalized waveforms (-1..1) defined with doubles
--> WAVE_ANALOG_32 = 1
- Analog 16 Bits Dual : Analog normalized waveforms (-1..1) defined with
doubles, with two components (A and B) --> WAVE_ANALOG_16_DUAL = 4
- Analog 32 Bits Dual : Analog normalized waveforms (-1..1) defined with
doubles, with two components (A and B) --> WAVE_ANALOG_32_DUAL = 6
- IQ : Analog normalized waveforms (-1..1) defined with doubles, with two
components (I and Q) --> WAVE_IQ = 2
- IQ Polar : Analog waveforms (-1..1 module, -180..+180 phase) defined with
doubles, with two components (Module and Phase) --> WAVE_IQPOLAR = 3
- Digital : Digital waveforms defined with integers --> WAVE_DIGITAL = 5
"""
## int SD_Wave_newFromArrayInteger(int waveformType, int waveformPoints, int *waveformDataA, int *waveformDataB SD_DEFAULT_NULL);
function SD_Wave_newFromArrayInteger(waveformType::Integer,
    waveformDataA::Vector, waveformDataB::Vector=[])
    c_waveformDataA = Vector{Cint}(waveformDataA)
    c_waveformDataB = Vector{Cint}(waveformDataB)
    waveformPoints = length(c_waveformDataA)
    val = ccall((:SD_Wave_newFromArrayInteger, klib), Cint,
        (Cint, Cint, Ref{Cint}, Ref{Cint}),
        waveformType, waveformPoints, c_waveformDataA, c_waveformDataB)
    return val
end

## int SD_Wave_newFromArrayDouble(int waveformType, int waveformPoints, double *waveformDataA, double *waveformDataB SD_DEFAULT_NULL);
function SD_Wave_newFromArrayDouble(waveformType::Integer,
    waveformDataA::Vector, waveformDataB::Vector=[])
    c_waveformDataA = Vector{Cdouble}(waveformDataA)
    c_waveformDataB = Vector{Cdouble}(waveformDataB)
    waveformPoints = length(c_waveformDataA)
    val = ccall((:SD_Wave_newFromArrayDouble, klib), Cint,
        (Cint, Cint, Ref{Cdouble}, Ref{Cdouble}),
        waveformType, waveformPoints, c_waveformDataA, c_waveformDataB)
    return val
end

## int SD_Wave_newFromFile(const char *waveformFile);
"""
This function creates a waveform object from data points contained in an array
in a file.
"""
SD_Wave_newFromFile(waveformFile::String) =
    ccall((:SD_Wave_newFromFile, klib), Cint, (Cstring,), waveformFile)

## int SD_Wave_delete(int waveformID);
"""
This function removes a waveform created with the function.
"""
SD_Wave_delete(waveformID::Integer) =
    ccall((:SD_Wave_delete, klib), Cint, (Cint,), waveformID)


text:
using Test
import LanguageServer: Document, Range, Position,
                       get_text, get_offset, get_line_offsets, get_position_at

@testset "UTF16 handling" begin
    doc = Document("", "aaa", false)
    @test sizeof(get_text(doc)) == 3
    @test get_offset(doc, 0, 0) == 0
    @test get_position_at(doc, 0) == (0, 0)
    @test get_offset(doc, 0, 1) == 1
    @test get_position_at(doc, 1) == (0, 1)
    @test get_offset(doc, 0, 2) == 2
    @test get_position_at(doc, 2) == (0, 2)
    @test get_offset(doc, 0, 3) == 3
    @test get_position_at(doc, 3) == (0, 3)

    doc = Document("", "ααα", false)
    @test sizeof(get_text(doc)) == 6
    @test get_offset(doc, 0, 0) == 0
    @test get_position_at(doc, 0) == (0, 0)
    @test get_offset(doc, 0, 1) == 1
    @test get_position_at(doc, 1) == (0, 1)
    @test get_offset(doc, 0, 2) == 3
    @test get_position_at(doc, 3) == (0, 2)
    @test get_offset(doc, 0, 3) == 5
    @test get_position_at(doc, 5) == (0, 3)

    doc = Document("", "ࠀࠀࠀ", false) # 0x0800
    @test sizeof(get_text(doc)) == 9
    @test get_offset(doc, 0, 0) == 0
    @test get_position_at(doc, 0) == (0, 0)
    @test get_offset(doc, 0, 1) == 1
    @test get_position_at(doc, 1) == (0, 1)
    @test get_offset(doc, 0, 2) == 4
    @test get_position_at(doc, 4) == (0, 2)
    @test get_offset(doc, 0, 3) == 7
    @test get_position_at(doc, 7) == (0, 3)

    doc = Document("", "𐐀𐐀𐐀", false)
    @test sizeof(get_text(doc)) == 12
    @test get_offset(doc, 0, 0) == 0
    @test get_position_at(doc, 0) == (0, 0)
    @test get_offset(doc, 0, 2) == 1
    @test get_position_at(doc, 1) == (0, 2)
    @test get_offset(doc, 0, 4) == 5
    @test get_position_at(doc, 5) == (0, 4)
    @test get_offset(doc, 0, 6) == 9
    @test get_position_at(doc, 9) == (0, 6)

    doc = Document("", "𐀀𐀀𐀀", false) # 0x010000
    @test sizeof(get_text(doc)) == 12
    @test get_offset(doc, 0, 0) == 0
    @test get_position_at(doc, 0) == (0, 0)
    @test get_offset(doc, 0, 2) == 1
    @test get_position_at(doc, 1) == (0, 2)
    @test get_offset(doc, 0, 4) == 5
    @test get_position_at(doc, 5) == (0, 4)
    @test get_offset(doc, 0, 6) == 9
    @test get_position_at(doc, 9) == (0, 6)
end


text:
######################################################################
# info.jl
#
# These functions used to be in scoring but they are unrelated
######################################################################

export collect_info, optimal_info

"""
`collect_info(em::ErgodicManager, traj::VVF, d_rate::Float64)`

`collect_info(em::,rgodicManager, traj::VVF)`

modifies em.phi according to some submodular.
we don't use the last point in the trajectory

decreases at rate `D/T`
If you spend `h` time there, it is equivalent to `h*D/(h*N) = D/N`

Returns total info picked up (a scalar value).
"""
function collect_info(em::ErgodicManager, traj::VVF; steps=0, d_rate=0.0)
	N = length(traj) - 1
	D = sum(em.phi)

	# if user has not specified decay rate, just set it to normal
	if d_rate == 0.0
		d_rate = D/(N+1)
	end

	total_info = 0.0
	if steps != 0
		N = steps
	end

	for n in 0:N
		# determine the cell that this point is in
		xi,yi = find_cell(em.domain, traj[n+1])

		# if there is enough info, grab it
		info_value = min(em.phi[xi,yi], d_rate)
		em.phi[xi,yi] -= info_value
		total_info += info_value
	end
	return total_info
end

function find_cell(em::ErgodicManager, x::Vector{Float64})
	return find_cell(em.bins, em.cell_size, x)
end
# should be the proper one
function find_cell(d::Domain, x::Vector{Float64})
	x1 = round(Int, x[1] / d.cell_lengths[1], RoundDown) + 1
	x2 = round(Int, x[2] / d.cell_lengths[2], RoundDown) + 1
	if x1 > d.cells[1]; x1 -= 1; end
	if x2 > d.cells[2]; x2 -= 1; end
	if x1 < 1; x1 += 1; end
	if x2 < 1; x2 += 1; end
	return x1, x2
end

"""
`optimal_info(em, N)`

Returns optimal information gathered.
"""
function optimal_info(em::ErgodicManager, N::Int)
	d_rate = sum(em.phi) / (N+1)
	total_info = 0.0

	# for each point in the trajectory...
	for n = 1:(N+1)
		best_ind = indmax(em.phi)
		info_value = min(em.phi[best_ind], d_rate)
		em.phi[best_ind] -= info_value
		total_info += info_value
	end

	return total_info
end

# create trajectory to do the above
export optimal_traj
function optimal_traj(em::ErgodicManager, tm::TrajectoryManager)
	d_rate = sum(em.phi) / (tm.N+1)
	total_info = 0.0
	xd = Array(Vector{Float64}, tm.N+1)
	xd[1] = deepcopy(tm.x0)
	temp_phi = deepcopy(em.phi)
	size_tuple = (x_cells(em), y_cells(em))
	for n = 1:tm.N
		bi = indmax(temp_phi)
		xi, yi = ind2sub(size_tuple, bi)
		xd[n+1] = [(xi-0.5)*x_size(em), (yi-0.5)*y_size(em)]
		temp_phi[bi] -= min(temp_phi[bi], d_rate)
	end
	ud = compute_controls(xd, tm.h)
	return xd,ud
end


text:
import ..UncertainValues:
    AbstractUncertainValue
import ..Resampling:
    resample

# Cosine functions
"""
    Base.cos(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the cosine of the uncertain value `x`, where `x` is in radians. Computes the
element-wise cosine for `n` realizations.
"""
Base.cos(x::AbstractUncertainValue; n::Int = 10000) = cos.(resample(x, n))

"""
    Base.cos(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the cosine of the uncertain value `x`, where `x` is in radians. Computes the
element-wise cosine for `n` realizations.
"""
Base.cos(x::AbstractUncertainValue, n::Int) = cos.(resample(x, n))

"""
    Base.cos(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the cosine of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise cosine for `n` realizations.
"""
Base.cosd(x::AbstractUncertainValue; n::Int = 10000) = cosd.(resample(x, n))

"""
    Base.cos(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the cosine of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise cosine for `n` realizations.
"""
Base.cosd(x::AbstractUncertainValue, n::Int) = cosd.(resample(x, n))

"""
    Base.cos(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic cosine of the uncertain value `x`.
Computes the element-wise hyperbolic cosine for `n` realizations.
"""
Base.cosh(x::AbstractUncertainValue; n::Int = 10000) = cosh.(resample(x, n))

"""
    Base.cos(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic cosine of the uncertain value `x`.
Computes the element-wise hyperbolic cosine for `n` realizations.
"""
Base.cosh(x::AbstractUncertainValue, n::Int) = cosh.(resample(x, n))

# Sine functions
"""
    Base.sin(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the sine of the uncertain value `x`, where `x` is in radians. Computes the
element-wise sine for `n` realizations.
"""
Base.sin(x::AbstractUncertainValue; n::Int = 10000) = sin.(resample(x, n))

"""
    Base.sin(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the sine of the uncertain value `x`, where `x` is in radians. Computes the
element-wise sine for `n` realizations.
"""
Base.sin(x::AbstractUncertainValue, n::Int) = sin.(resample(x, n))

"""
    Base.sind(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the sine of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise sine for `n` realizations.
"""
Base.sind(x::AbstractUncertainValue; n::Int = 10000) = sind.(resample(x, n))

"""
    Base.sind(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the sine of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise sine for `n` realizations.
"""
Base.sind(x::AbstractUncertainValue, n::Int) = sind.(resample(x, n))

"""
    Base.sinh(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic sine of the uncertain value `x`.
Computes the element-wise hyperbolic sine for `n` realizations.
"""
Base.sinh(x::AbstractUncertainValue; n::Int = 10000) = sinh.(resample(x, n))

"""
    Base.sinh(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic sine of the uncertain value `x`.
Computes the element-wise hyperbolic sine for `n` realizations.
"""
Base.sinh(x::AbstractUncertainValue, n::Int) = sinh.(resample(x, n))

# Tan functions
"""
    Base.tan(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the tangent of the uncertain value `x`, where `x` is in radians. Computes the
element-wise tangent for `n` realizations.
"""
Base.tan(x::AbstractUncertainValue; n::Int = 10000) = tan.(resample(x, n))

"""
    Base.tan(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the tangent of the uncertain value `x`, where `x` is in radians. Computes the
element-wise tangent for `n` realizations.
"""
Base.tan(x::AbstractUncertainValue, n::Int) = tan.(resample(x, n))

"""
    Base.tand(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the tangent of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise tangent for `n` realizations.
"""
Base.tand(x::AbstractUncertainValue; n::Int = 10000) = tand.(resample(x, n))

"""
    Base.tand(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the tangent of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise tangent for `n` realizations.
"""
Base.tand(x::AbstractUncertainValue, n::Int) = tand.(resample(x, n))

"""
    Base.tanh(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic tangent of the uncertain value `x`.
Computes the element-wise hyperbolic tangent for `n` realizations.
"""
Base.tanh(x::AbstractUncertainValue; n::Int = 10000) = tanh.(resample(x, n))

"""
    Base.tanh(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic tangent of the uncertain value `x`.
Computes the element-wise hyperbolic tangent for `n` realizations.
"""
Base.tanh(x::AbstractUncertainValue, n::Int) = tanh.(resample(x, n))

# Other

"""
    Base.sincos(x::AbstractUncertainValue; n::Int = 10000)

Simultaneously compute the sine and cosine of the uncertain value `x`, where `x` is in
radians. Computes the element-wise `sincos` for `n` realizations.
"""
Base.sincos(x::AbstractUncertainValue; n::Int = 10000) = sincos.(resample(x, n))

"""
    Base.sincos(x::AbstractUncertainValue, n::Int = 10000)

Simultaneously compute the sine and cosine of the uncertain value `x`, where `x` is in
radians. Computes the element-wise `sincos` for `n` realizations.
"""
Base.sincos(x::AbstractUncertainValue, n::Int) = sincos.(resample(x, n))

"""
    Base.sinc(x::AbstractUncertainValue; n::Int = 10000)

In an element-wise manner for `n` realizations of the uncertain value `x`, compute
``\\sin(\\pi x) / (\\pi x)`` if ``x \\neq 0``, and ``1`` if ``x = 0``.
"""
Base.sinc(x::AbstractUncertainValue; n::Int = 10000) = sinc.(resample(x, n))

"""
    Base.sinc(x::AbstractUncertainValue, n::Int = 10000)

Compute ``\\sin(\\pi x) / (\\pi x)`` if ``x \\neq 0``, and ``1`` if ``x = 0`` element-wise
over `n` realizations of the uncertain value `x`.
"""
Base.sinc(x::AbstractUncertainValue, n::Int) = sinc.(resample(x, n))

"""
    Base.sinpi(x::AbstractUncertainValue; n::Int = 10000)

Compute ``\\sin(\\pi x)`` more accurately than `sin(pi*x)`, especially for large `x`,
in an element-wise over `n` realizations of the uncertain value `x`.
"""
Base.sinpi(x::AbstractUncertainValue; n::Int = 10000) = sinpi.(resample(x, n))

"""
    Base.sinpi(x::AbstractUncertainValue; n::Int = 10000)

Compute ``\\sin(\\pi x)`` more accurately than `sin(pi*x)`, especially for large `x`,
in an element-wise over `n` realizations of the uncertain value `x`.
"""
Base.sinpi(x::AbstractUncertainValue, n::Int) = sinpi.(resample(x, n))

"""
    Base.cosc(x::AbstractUncertainValue; n::Int = 10000)

Compute ``\\cos(\\pi x) / x - \\sin(\\pi x) / (\\pi x^2)`` if ``x \\neq 0``, and ``0`` if
``x = 0``, in an element-wise manner over `n` realizations of the uncertain value `x`.

This is the derivative of `sinc(x)`.
"""
Base.cosc(x::AbstractUncertainValue; n::Int = 10000) = cosc.(resample(x, n))

"""
    Base.cosc(x::AbstractUncertainValue, n::Int = 10000)

Compute ``\\cos(\\pi x) / x - \\sin(\\pi x) / (\\pi x^2)`` if ``x \\neq 0``, and ``0`` if
``x = 0``, in an element-wise manner over `n` realizations of the uncertain value `x`.

This is the derivative of `sinc(x)`.
"""
Base.cosc(x::AbstractUncertainValue, n::Int) = cosc.(resample(x, n))

"""
    Base.cospi(x::AbstractUncertainValue; n::Int = 10000)

Compute ``\\cos(\\pi x)`` more accurately than `cos(pi*x)`, especially for large `x`,
in an element-wise over `n` realizations of the uncertain value `x`.
"""
Base.cospi(x::AbstractUncertainValue; n::Int = 10000) = cospi.(resample(x, n))

"""
    Base.cospi(x::AbstractUncertainValue, n::Int = 10000)

Compute ``\\cos(\\pi x)`` more accurately than `cos(pi*x)`, especially for large `x`,
in an element-wise over `n` realizations of the uncertain value `x`.
"""
Base.cospi(x::AbstractUncertainValue, n::Int) = cospi.(resample(x, n))

############################
# Inverse trig functions
############################
"""
    Base.acos(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse cosine of the uncertain value `x`, where `x` is in radians. Computes the
element-wise inverse cosine for `n` realizations.
"""
Base.acos(x::AbstractUncertainValue; n::Int = 10000) = acos.(resample(x, n))

"""
    Base.acos(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse cosine of the uncertain value `x`, where `x` is in radians. Computes the
element-wise inverse cosine for `n` realizations.
"""
Base.acos(x::AbstractUncertainValue, n::Int) = acos.(resample(x, n))

"""
    Base.acosd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse cosine of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise inverse cosine for `n` realizations.
"""
Base.acosd(x::AbstractUncertainValue; n::Int = 10000) = acosd.(resample(x, n))

"""
    Base.acosd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse cosine of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise inverse cosine for `n` realizations.
"""
Base.acosd(x::AbstractUncertainValue, n::Int) = acosd.(resample(x, n))

"""
    Base.acosh(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic cosine of the uncertain value `x`.
Computes the element-wise inverse hyperbolic cosine for `n` realizations.
"""
Base.acosh(x::AbstractUncertainValue; n::Int = 10000) = acosh.(resample(x, n))

"""
    Base.acosh(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic cosine of the uncertain value `x`.
Computes the element-wise inverse hyperbolic cosine for `n` realizations.
"""
Base.acosh(x::AbstractUncertainValue, n::Int) = acosh.(resample(x, n))

"""
    Base.asin(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse sine of the uncertain value `x`, where `x` is in radians. Computes the
element-wise inverse sine for `n` realizations.
"""
Base.asin(x::AbstractUncertainValue; n::Int = 10000) = asin.(resample(x, n))

"""
    Base.asin(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse sine of the uncertain value `x`, where `x` is in radians. Computes the
element-wise inverse sine for `n` realizations.
"""
Base.asin(x::AbstractUncertainValue, n::Int) = asin.(resample(x, n))

"""
    Base.asind(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse sine of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise inverse sine for `n` realizations.
"""
Base.asind(x::AbstractUncertainValue; n::Int = 10000) = asind.(resample(x, n))

"""
    Base.asind(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse sine of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise inverse sine for `n` realizations.
"""
Base.asind(x::AbstractUncertainValue, n::Int) = asind.(resample(x, n))

"""
    Base.asinh(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic sine of the uncertain value `x`.
Computes the element-wise inverse hyperbolic csine for `n` realizations.
"""
Base.asinh(x::AbstractUncertainValue; n::Int = 10000) = asinh.(resample(x, n))

"""
    Base.asinh(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic sine of the uncertain value `x`.
Computes the element-wise inverse hyperbolic csine for `n` realizations.
"""
Base.asinh(x::AbstractUncertainValue, n::Int) = asinh.(resample(x, n))

"""
    Base.atan(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse tangent of the uncertain value `x`, where `x` is in radians.
Computes the element-wise inverse tangent for `n` realizations.
"""
Base.atan(x::AbstractUncertainValue; n::Int = 10000) = atan.(resample(x, n))

"""
    Base.atan(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse tangent of the uncertain value `x`, where `x` is in radians.
Computes the element-wise inverse tangent for `n` realizations.
"""
Base.atan(x::AbstractUncertainValue, n::Int) = atan.(resample(x, n))

"""
    Base.atand(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse tangent of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise inverse tangent for `n` realizations.
"""
Base.atand(x::AbstractUncertainValue; n::Int = 10000) = atand.(resample(x, n))

"""
    Base.atand(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse tangent of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise inverse tangent for `n` realizations.
"""
Base.atand(x::AbstractUncertainValue, n::Int) = atand.(resample(x, n))

"""
    Base.atanh(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse hypoerbolic tangent of the uncertain value `x`.
Computes the element-wise inverse hypoerbolic tangent for `n` realizations.
"""
Base.atanh(x::AbstractUncertainValue; n::Int = 10000) = atanh.(resample(x, n))

"""
    Base.atanh(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse hypoerbolic tangent of the uncertain value `x`.
Computes the element-wise inverse hypoerbolic tangent for `n` realizations.
"""
Base.atanh(x::AbstractUncertainValue, n::Int) = atanh.(resample(x, n))

"""
    Base.asec(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse secant of the uncertain value `x`, where `x` is in radians.
Computes the element-wise inverse secant for `n` realizations.
"""
Base.asec(x::AbstractUncertainValue; n::Int = 10000) = asec.(resample(x, n))

"""
    Base.asec(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse secant of the uncertain value `x`, where `x` is in radians.
Computes the element-wise inverse secant for `n` realizations.
"""
Base.asec(x::AbstractUncertainValue, n::Int) = asec.(resample(x, n))

"""
    Base.asecd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse secant of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise inverse secant for `n` realizations.
"""
Base.asecd(x::AbstractUncertainValue; n::Int = 10000) = asec.(resample(x, n))

"""
    Base.asecd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse secant of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise inverse secant for `n` realizations.
"""
Base.asecd(x::AbstractUncertainValue, n::Int) = asec.(resample(x, n))

"""
    Base.asech(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic secant of the uncertain value `x`.
Computes the element-wise inverse hyperbolic secant for `n` realizations.
"""
Base.asech(x::AbstractUncertainValue; n::Int = 10000) = asec.(resample(x, n))

"""
    Base.asech(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic secant of the uncertain value `x`.
Computes the element-wise inverse hyperbolic secant for `n` realizations.
"""
Base.asech(x::AbstractUncertainValue, n::Int) = asec.(resample(x, n))

"""
    Base.acsc(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse cosecant of the uncertain value `x`, where `x` is in radians.
Computes the element-wise inverse cosecant for `n` realizations.
"""
Base.acsc(x::AbstractUncertainValue; n::Int = 10000) = acsc.(resample(x, n))

"""
    Base.acsc(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse cosecant of the uncertain value `x`, where `x` is in radians.
Computes the element-wise inverse cosecant for `n` realizations.
"""
Base.acsc(x::AbstractUncertainValue, n::Int) = acsc.(resample(x, n))

"""
    Base.acscd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse cosecant of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise inverse cosecant for `n` realizations.
"""
Base.acscd(x::AbstractUncertainValue; n::Int = 10000) = acsc.(resample(x, n))

"""
    Base.acscd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse cosecant of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise inverse cosecant for `n` realizations.
"""
Base.acscd(x::AbstractUncertainValue, n::Int) = acsc.(resample(x, n))

"""
    Base.acscd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic cosecant of the uncertain value `x`.
Computes the element-wise inverse hypoerbolic cosecant for `n` realizations.
"""
Base.acsch(x::AbstractUncertainValue; n::Int = 10000) = acsc.(resample(x, n))

"""
    Base.acscd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic cosecant of the uncertain value `x`.
Computes the element-wise inverse hypoerbolic cosecant for `n` realizations.
"""
Base.acsch(x::AbstractUncertainValue, n::Int) = acsc.(resample(x, n))

"""
    Base.acot(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse cotangent of the uncertain value `x`, where `x` is in radians.
Computes the element-wise inverse secant for `n` realizations.
"""
Base.acot(x::AbstractUncertainValue; n::Int = 10000) = acot.(resample(x, n))

"""
    Base.acot(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse cotangent of the uncertain value `x`, where `x` is in radians.
Computes the element-wise inverse secant for `n` realizations.
"""
Base.acot(x::AbstractUncertainValue, n::Int) = acot.(resample(x, n))

"""
    Base.acotd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse cotangent of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise inverse secant for `n` realizations.
"""
Base.acotd(x::AbstractUncertainValue; n::Int = 10000) = acot.(resample(x, n))

"""
    Base.acotd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse cotangent of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise inverse secant for `n` realizations.
"""
Base.acotd(x::AbstractUncertainValue, n::Int) = acot.(resample(x, n))

"""
    Base.acoth(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic cotangent of the uncertain value `x`.
Computes the element-wise inverse hyperbolic secant for `n` realizations.
"""
Base.acoth(x::AbstractUncertainValue; n::Int = 10000) = acot.(resample(x, n))

"""
    Base.acoth(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the inverse hyperbolic cotangent of the uncertain value `x`.
Computes the element-wise inverse hyperbolic secant for `n` realizations.
"""
Base.acoth(x::AbstractUncertainValue, n::Int) = acot.(resample(x, n))

############################
# Reciprocal trig functions
############################

"""
    Base.csc(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the cosecant of the uncertain value `x`, where `x` is in radians. Computes the
element-wise cosecant for `n` realizations.
"""
Base.csc(x::AbstractUncertainValue; n::Int = 10000) = csc.(resample(x, n))

"""
    Base.csc(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the cosecant of the uncertain value `x`, where `x` is in radians. Computes the
element-wise cosecant for `n` realizations.
"""
Base.csc(x::AbstractUncertainValue, n::Int) = csc.(resample(x, n))

"""
    Base.cscd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the cosecant of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise cosecant for `n` realizations.
"""
Base.cscd(x::AbstractUncertainValue; n::Int = 10000) = cscd.(resample(x, n))

"""
    Base.cscd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the cosecant of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise cosecant for `n` realizations.
"""
Base.cscd(x::AbstractUncertainValue, n::Int) = cscd.(resample(x, n))

"""
    Base.cscd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic cosecant of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise hyperbolic cosecant for `n` realizations.
"""
Base.csch(x::AbstractUncertainValue; n::Int = 10000) = csch.(resample(x, n))

"""
    Base.cscd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic cosecant of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise hyperbolic cosecant for `n` realizations.
"""
Base.csch(x::AbstractUncertainValue, n::Int) = csch.(resample(x, n))

"""
    Base.sec(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the secant of the uncertain value `x`, where `x` is in radians. Computes the
element-wise secant for `n` realizations.
"""
Base.sec(x::AbstractUncertainValue; n::Int = 10000) = sec.(resample(x, n))

"""
    Base.sec(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the secant of the uncertain value `x`, where `x` is in radians. Computes the
element-wise secant for `n` realizations.
"""
Base.sec(x::AbstractUncertainValue, n::Int) = sec.(resample(x, n))

"""
    Base.secd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the secant of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise cosecant for `n` realizations.
"""
Base.secd(x::AbstractUncertainValue; n::Int = 10000) = secd.(resample(x, n))

"""
    Base.secd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the secant of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise cosecant for `n` realizations.
"""
Base.secd(x::AbstractUncertainValue, n::Int) = secd.(resample(x, n))

"""
    Base.sech(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic secant of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise hyperbolic secant for `n` realizations.
"""
Base.sech(x::AbstractUncertainValue; n::Int = 10000) = sech.(resample(x, n))

"""
    Base.sech(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic secant of the uncertain value `x`, where `x` is in degrees.
Computes the element-wise hyperbolic secant for `n` realizations.
"""
Base.sech(x::AbstractUncertainValue, n::Int) = sech.(resample(x, n))

"""
    Base.cot(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the cotangent of the uncertain value `x`, where `x` is in radians. Computes the
element-wise cotangent for `n` realizations.
"""
Base.cot(x::AbstractUncertainValue; n::Int = 10000) = cot.(resample(x, n))

"""
    Base.cot(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the cotangent of the uncertain value `x`, where `x` is in radians. Computes the
element-wise cotangent for `n` realizations.
"""
Base.cot(x::AbstractUncertainValue, n::Int) = cot.(resample(x, n))

"""
    Base.cotd(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the cotangent of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise cotangent for `n` realizations.
"""
Base.cotd(x::AbstractUncertainValue; n::Int = 10000) = cotd.(resample(x, n))

"""
    Base.cotd(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the cotangent of the uncertain value `x`, where `x` is in degrees. Computes the
element-wise cotangent for `n` realizations.
"""
Base.cotd(x::AbstractUncertainValue, n::Int) = cotd.(resample(x, n))

"""
    Base.coth(x::AbstractUncertainValue; n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic cotangent of the uncertain value `x`, where `x` is in radians.
Computes the element-wise hyperbolic cotangent for `n` realizations.
"""
Base.coth(x::AbstractUncertainValue; n::Int = 10000) = coth.(resample(x, n))

"""
    Base.coth(x::AbstractUncertainValue, n::Int = 10000) -> Vector{Float64}

Compute the hyperbolic cotangent of the uncertain value `x`, where `x` is in radians.
Computes the element-wise hyperbolic cotangent for `n` realizations.
"""
Base.coth(x::AbstractUncertainValue, n::Int) = coth.(resample(x, n))

#####################################################################################
# Special cases
#####################################################################################

trigfuncs = [:(cos), :(cosd), :(cosh), :(sin), :(sind), :(sinh), :(tan), :(tand), :(tanh),
    :(sincos), :(sinc), :(sinpi), :(cosc), :(cospi), :(acos), :(acosd), :(acosh), :(asin),
    :(asind), :(asinh), :(atan), :(atand), :(atanh), :(asec), :(asecd), :(asech), :(acsc),
    :(acscd), :(acsch), :(acot), :(acotd), :(acoth), :(csc), :(cscd), :(csch), :(sec),
    :(secd), :(sech), :(cot), :(cotd), :(coth)]

##################
# `CertainValue`s
#################
import ..UncertainValues: CertainValue

"""
    Base.:-(a::Union{CertainValue, Real}, b::Union{CertainValue, Real})

Subtraction of certain values with themselves or scalars acts as regular subtraction,
but returns the result wrapped in a `CertainValue` instance.
"""

for trigfunc in trigfuncs
    f = Meta.parse("Base.$(trigfunc)")

    regular_func = quote
        """
            $($f)(x::CertainValue)

        Compute `$($trigfunc)(x)`.
        """
        $(f)(x::CertainValue) = x.value
    end

    n_reps_func = quote
        """
            $($f)(x::CertainValue, n::Int)

        Compute `$($trigfunc)(x)` `n` times and return the result(s) as a vector (just
        repeating the value `n` times).
        """
        $(f)(x::CertainValue, n::Int) = [x.value for i = 1:n]
    end

    eval(regular_func)
    eval(n_reps_func)
end


text:
"""
    pdi{T<:Number}(x::Array{T, 1})

This function will range the values of each row, so that the strongest link has
a value of one. This works for deterministic and quantitative networks.

#### References

Poisot, T., Bever, J.D., Nemri, A., Thrall, P.H., Hochberg, M.E., 2011. A
conceptual framework for the evolution of ecological specialisation. Ecol. Lett.
14, 841–851. https://doi.org/10.1111/j.1461-0248.2011.01645.x

Poisot, T., Canard, E., Mouquet, N., Hochberg, M.E., 2012. A comparative study
of ecological specialization estimators. Methods in Ecology and Evolution 3,
537–544. https://doi.org/10.1111/j.2041-210X.2011.00174.x
"""
function pdi(x::Vector{T}) where {T<:Number}
    s = reverse(sort(x))
    s = s ./ maximum(s)
    p = s[1] .- s
    return sum(p[2:end])/(length(s)-1)
end

"""
    specificity(N::DeterministicNetwork)

Measure of specificity in a deterministic network. This returns a value between
0 and 1, where 1 indicates maximal specificity.

#### References

Poisot, T., Bever, J.D., Nemri, A., Thrall, P.H., Hochberg, M.E., 2011. A
conceptual framework for the evolution of ecological specialisation. Ecol. Lett.
14, 841–851. https://doi.org/10.1111/j.1461-0248.2011.01645.x

Poisot, T., Canard, E., Mouquet, N., Hochberg, M.E., 2012. A comparative study
of ecological specialization estimators. Methods in Ecology and Evolution 3,
537–544. https://doi.org/10.1111/j.2041-210X.2011.00174.x
"""
function specificity(N::BinaryNetwork)
    A = map(Int64, N.A)
    return Dict(zip(species(N; dims=1), vec(mapslices(pdi, A; dims=2))))
end

"""
    specificity(N::QuantitativeNetwork)

Measure of specificity in a quantitative network. This returns a value between 0
and 1, where 1 indicates maximal specificity. Note that the PDI is measured
species-wise, and the maximal interaction strength of every species is set to 1.

#### References

Poisot, T., Bever, J.D., Nemri, A., Thrall, P.H., Hochberg, M.E., 2011. A
conceptual framework for the evolution of ecological specialisation. Ecol. Lett.
14, 841–851. https://doi.org/10.1111/j.1461-0248.2011.01645.x

Poisot, T., Canard, E., Mouquet, N., Hochberg, M.E., 2012. A comparative study
of ecological specialization estimators. Methods in Ecology and Evolution 3,
537–544. https://doi.org/10.1111/j.2041-210X.2011.00174.x
"""
function specificity(N::QuantitativeNetwork)
    return Dict(zip(species(N; dims=1), vec(mapslices(pdi, N.A; dims=2))))
end


text:
using MLStyle

module Linq
    select(arr, f) = map(f, arr)
end

macro linq(expr)
    @match expr begin
        :($subject.$method($(args...))) =>
            let method = getfield(Linq, method)
                quote $method($subject, $(args...)) end
            end
        _ => @error "invalid"
    end
end

@test (@linq [1, 2, 3].select(x -> x * 5)) == [5, 10,  15]

@data internal AAA begin
    AAA_1(Int, Int)
    AAA_2(Float32)
end

@test (@linq [AAA_1(2, 2), AAA_2(3.9)].select(@λ begin
    AAA_1(a, b) -> a + b
    AAA_2(a)    -> 3
end)) == [4, 3]


text:
using Test
using FixedSizeArrays
CM = CompScienceMeshes

b = CM.sortneighbors([4,3,2,1], (x,y)-> y == x+1)
@test b == [1,2,3,4]

b = CM.sortneighbors([1,2,3,4], (x,y)->y == x+1)
@test b == [1,2,3,4]

b = CM.sortneighbors([4,5,6,1,2,3], (x,y)->y == x+1)
@test b == [1,2,3,4,5,6]

b = CM.sortneighbors([4,5,6,3,2,1], (x,y)->y == x+1)
@test b == [1,2,3,4,5,6]

b = CM.sortneighbors(randperm(10), (x,y)-> y == x+1)
@test b == collect(1:10)

a = [(1,2), (3,4), (2,3), (0,1)]
pred = (x,y) -> x[2] == y[1] || x[1] == y[2]
b = CM.sortneighbors(a, pred)
@test b == [(0,1), (1,2), (2,3), (3,4)]

mesh = CM.meshrectangle(1.0, 1.0, 0.5)
vtof, vton = CM.vertextocellmap(mesh)
supp = mesh.faces[vec(vtof[5,1:vton[5]])]
commonedge(x,y) = length(intersect(x,y)) == 2
supp = CM.sortneighbors(supp, commonedge)
@test length(supp) == 6
@test supp == Vec{3,Int}[
    Vec(2,5,4),
    Vec(2,3,5),
    Vec(3,6,5),
    Vec(5,6,8),
    Vec(5,8,7),
    Vec(4,5,7)
]
@test CM.isclosed(supp, commonedge) == true


text:
@testset "FeedbackChains" begin
    @testset "constructor" begin
        # 1. Can construct a FeedbackChain just like a normal chain.
        c = FeedbackChain(Dense(10, 5), relu)
        @test length(c.layers) == 2
        # 2. Can construct a FeedbackChain with forkers and mergers
        c = FeedbackChain(Merger("f1", identity, +), Dense(10, 10, relu), Splitter("f1"), Dense(10, 1))
        @test length(c.layers) == 4
    end # @testset "constructor"

    @testset "apply" begin
        h = Dict("f1"=>randn(10))
        x = randn(10)
        # apply empty chain
        c = FeedbackChain()
        @test c(h, x) == (h, x)
        # apply FeedbackChain without feedback
        c = FeedbackChain(x -> (2.0 .* x), x -> relu.(x), x -> x'x)
        state, y = c(Dict(), x)
        @test y ≈ relu.(2.0 .* x)'relu.(2.0 .* x)
        @test state isa Dict{String, Any}
        # apply feedback chain with feedback
        l1 = Dense(10, 10, relu)
        l2 = Dense(10, 1)
        c = FeedbackChain(Merger("f1", identity, +), l1, Splitter("f1"), l2)
        state, y = c(h, x)
        @test y ≈ l2(l1(h["f1"] + x))
        @test haskey(state, "f1")
        @test state["f1"] ≈ l1(h["f1"] + x)
    end # @testset "apply"

    @testset "params" begin
        l1 = Dense(10, 10)
        l2 = Dense(7, 10)
        c = FeedbackChain(l1, l2)
        @test hcat(params(c)...) == hcat(params(l1)..., params(l2)...)
    end # @testset "params"

    @testset "statefull" begin
        h = Dict("f1"=>randn(10))
        x = randn(10)
        l1 = Dense(10, 10, relu)
        l2 = Dense(10, 1)
        c = FeedbackChain(Merger("f1", identity, +), l1, Splitter("f1"), l2)
        c = StateFull(c, h)
        y = c(x)
        @test y ≈ l2(l1(h["f1"] + x))
        @test haskey(c.state, "f1")
        @test c.state["f1"] ≈ l1(h["f1"] + x)
        # check that truncation works
        # !!! Not clear how to do this with new Flux version !!!
        #@test Flux.Tracker.istracked(c.state["f1"])
        #Flux.truncate!(c)
        #@test !Flux.Tracker.istracked(c.state["f1"])
        #
        # check that reset works
        Flux.reset!(c)
        @test c.state["f1"] == c.init["f1"]
    end # @testset "statefull"

    @testset "indexing" begin
        c = FeedbackChain(Dense(1,2), Dense(2,3), Dense(3,1))
        @test c[1] == c.layers[1]
        @test c[2:3] == FeedbackChain(c.layers[2:3]...)
    end # @testset "indexing"

    @testset "iteration" begin
        c = FeedbackChain(Dense(1,2), Dense(2,3), Dense(3,1))
        layers = [layer for layer in c]
        @test layers == collect(c.layers)
    end # @testset "iteration"

    @testset "interface" begin
        # test that generic method `splitname` works
        c = FeedbackChain(
            Splitter("name1"),
            Merger("name2", nothing, +),
            Splitter("name3"),
            Dense(10,10)
        )
        @test splitnames(c) == ["name1", "name3"]
        # test that validation of naming works
        @test !namesvalid(c) # merger name does not fit splitter names
        c = FeedbackChain(Splitter("name1"), Merger("name1", nothing, +))
        @test namesvalid(c)
        c = FeedbackChain(Splitter("name1"), Splitter("name1"))
        @test !namesvalid(c)
    end # @testset "interface"
end # @testset "FeedbackChains"


text:
include("infdecomp.jl")
using InfDecomp

function my_cholesky!{FLOAT}(A::Matrix{FLOAT}) :: None
    @assert size(A)[1]==size(A)[2] "my_cholesky(): Need square matrix!"
    @assert size(A)[1] >=2         "my_cholesky(): Don't fuck with me!"

    const n::Int = size(A)[1]

    for ell in 1:n
        if A[ell,ell] <= 0
            @show ell A[ell,ell]
            @assert A[ell,ell] > 0
        end
        const d::FLOAT = sqrt(A[ell,ell])
        A[ell,ell] = d
        for k = (ell+1):n
            const c::FLOAT = A[k,ell] / d
            A[k,ell] = c
            for i = (ell+1):k
                A[k,i] -= c*A[i,ell]
            end # j
        end # k
    end # ell
    #
    for k in 1:n
        for ell in (k+1):n
            A[k,ell] = A[ell,k]
        end
    end
    return nothing;
end

import InfDecomp.My_Eval

# This function returns a pair (γ,H); the hessian is:  1/γ * H
function hessian_block{FLOAT,TFloat}(e::My_Eval, p::Vector{TFloat}, y::Int, z::Int, dummy::FLOAT) :: Tuple{FLOAT,Matrix{FLOAT}}
    @assert 1 <= y <= e.n_y
    @assert 1 <= z <= e.n_z

    P_yz  ::FLOAT = FLOAT(0.)
    for x = 1:e.n_x
        i = e.varidx[x,y,z]
        if i>0
            P_yz += p[i]
        end
    end

    H = zeros(FLOAT,e.n_x,e.n_x)

    # Start with the diagonal
    for x = 1:e.n_x
        i = e.varidx[x,y,z]
        if i>0
            P_xyz = FLOAT( p[i] )
            H[x,x] = P_yz/P_xyz - 1
        end
    end
    # Now off-diagonal.
    for x = 1:e.n_x
        for u = 1:(x-1)
            i_x = e.varidx[x,y,z]
            i_u = e.varidx[u,y,z]
            if i_x>0 && i_u>0
                counter += 1
                H[x,u] = FLOAT(-1)
                H[u,x] = FLOAT(-1)
            end
        end
    end

    return (P_yz,H)
    ;
end

The Quadratic Program
---------------------

[ a ]     [ H   -A^T ] [ u ]
[   ] ==  [          ] [   ]
[ b ]     [ A     0  ] [ v ]

For Schur-complement method, need to solve:

  A U (U^T H U)^{-1} U^T A^T v == c

where U is the isometry taking im(A^T) into the full space; H is invertible on that space.

The equation

   AU w = d

has a unique solution w in im(A^T).

Then solve A^T v = H w    for v.

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
### Process examples
# Always rerun examples
const EXAMPLES_OUT = joinpath(@__DIR__, "src", "examples")
ispath(EXAMPLES_OUT) && rm(EXAMPLES_OUT; recursive=true)
mkpath(EXAMPLES_OUT)

# Install and precompile all packages
# Workaround for https://github.com/JuliaLang/Pkg.jl/issues/2219
examples = filter!(isdir, readdir(joinpath(@__DIR__, "..", "examples"); join=true))
let script = "using Pkg; Pkg.activate(ARGS[1]); Pkg.instantiate()"
    for example in examples
        if !success(`$(Base.julia_cmd()) -e $script $example`)
            error(
                "project environment of example ",
                basename(example),
                " could not be instantiated",
            )
        end
    end
end
# Run examples asynchronously
processes = let literatejl = joinpath(@__DIR__, "literate.jl")
    map(examples) do example
        return run(
            pipeline(
                `$(Base.julia_cmd()) $literatejl $(basename(example)) $EXAMPLES_OUT`;
                stdin=devnull,
                stdout=devnull,
                stderr=stderr,
            );
            wait=false,
        )::Base.Process
    end
end

# Check that all examples were run successfully
isempty(processes) || success(processes) || error("some examples were not run successfully")

# Build documentation
using CalibrationErrors
using Documenter

makedocs(;
    modules=[CalibrationErrors],
    authors="David Widmann <david.widmann@it.uu.se>",
    repo="https://github.com/devmotion/CalibrationErrors.jl/blob/{commit}{path}#L{line}",
    sitename="CalibrationErrors.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://devmotion.github.io/CalibrationErrors.jl",
        assets=String[],
    ),
    pages=[
        "index.md",
        "introduction.md",
        "ece.md",
        "kce.md",
        "others.md",
        "Examples" =>
            map(filter!(filename -> endswith(filename, ".md"), readdir(EXAMPLES_OUT))) do x
                return joinpath("examples", x)
            end,
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/devmotion/CalibrationErrors.jl.git",
    push_preview=true,
    devbranch="main",
)


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__

d_str = "dag{k1 -> k2;k2 -> y;v -> x2;w -> k1;x1 -> v;x1 -> w;x2 -> y;x3 -> w;x3 -> y}"

dag = DAG("test_desc", d_str)

to_dagitty(dag)

#=
fname = joinpath(ProjDir, "test_descendents_03.dot")
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)
=#

bs = basis_set(dag)

@testset "basis_set_01" begin

  @test length(bs) == 48
  @test bs[18] == [:x2, :w, :v]
  @test bs[19] == [:x2, :w, :x1]
  @test bs[20] == [:x2, :w, :v, :x1]
  @test bs[44] == [:y, :k1, :w, :k2]
  @test bs[45] == [:y, :k1, :w, :x2, :k2]
  @test bs[46] == [:y, :k1, :w, :x3, :k2]
  @test bs[47] == [:y, :k1, :x2, :x3, :k2]
  @test bs[48] == [:y, :k1, :w, :x2, :x3, :k2]

end


text:
#####
##### `Symmetric`/`Hermitian`
#####

function frule((_, ΔA, _), T::Type{<:LinearAlgebra.HermOrSym}, A::AbstractMatrix, uplo)
    return T(A, uplo), T(ΔA, uplo)
end

function rrule(T::Type{<:LinearAlgebra.HermOrSym}, A::AbstractMatrix, uplo)
    Ω = T(A, uplo)
    function HermOrSym_pullback(ΔΩ)
        return (NO_FIELDS, _symherm_back(typeof(Ω), ΔΩ, Ω.uplo), DoesNotExist())
    end
    return Ω, HermOrSym_pullback
end

function frule((_, ΔA), TM::Type{<:Matrix}, A::LinearAlgebra.HermOrSym)
    return TM(A), TM(_symherm_forward(A, ΔA))
end
function frule((_, ΔA), ::Type{Array}, A::LinearAlgebra.HermOrSym)
    return Array(A), Array(_symherm_forward(A, ΔA))
end

function rrule(TM::Type{<:Matrix}, A::LinearAlgebra.HermOrSym)
    function Matrix_pullback(ΔΩ)
        TA = _symhermtype(A)
        T∂A = TA{eltype(ΔΩ),typeof(ΔΩ)}
        uplo = A.uplo
        ∂A = T∂A(_symherm_back(A, ΔΩ, uplo), uplo)
        return NO_FIELDS, ∂A
    end
    return TM(A), Matrix_pullback
end
rrule(::Type{Array}, A::LinearAlgebra.HermOrSym) = rrule(Matrix, A)

# Get type (Symmetric or Hermitian) from type or matrix
_symhermtype(::Type{<:Symmetric}) = Symmetric
_symhermtype(::Type{<:Hermitian}) = Hermitian
_symhermtype(A) = _symhermtype(typeof(A))

# for Ω = Matrix(A::HermOrSym), push forward ΔA to get ∂Ω
function _symherm_forward(A, ΔA)
    TA = _symhermtype(A)
    return if ΔA isa TA
        ΔA
    else
        TA{eltype(ΔA),typeof(ΔA)}(ΔA, A.uplo)
    end
end

# for Ω = HermOrSym(A, uplo), pull back ΔΩ to get ∂A
_symherm_back(::Type{<:Symmetric}, ΔΩ, uplo) = _symmetric_back(ΔΩ, uplo)
function _symherm_back(::Type{<:Hermitian}, ΔΩ::AbstractMatrix{<:Real}, uplo)
    return _symmetric_back(ΔΩ, uplo)
end
_symherm_back(::Type{<:Hermitian}, ΔΩ, uplo) = _hermitian_back(ΔΩ, uplo)
_symherm_back(Ω, ΔΩ, uplo) = _symherm_back(typeof(Ω), ΔΩ, uplo)

function _symmetric_back(ΔΩ, uplo)
    L, U, D = LowerTriangular(ΔΩ), UpperTriangular(ΔΩ), Diagonal(ΔΩ)
    return uplo == 'U' ? U .+ transpose(L) - D : L .+ transpose(U) - D
end
_symmetric_back(ΔΩ::Diagonal, uplo) = ΔΩ
_symmetric_back(ΔΩ::UpperTriangular, uplo) = Matrix(uplo == 'U' ? ΔΩ : transpose(ΔΩ))
_symmetric_back(ΔΩ::LowerTriangular, uplo) = Matrix(uplo == 'U' ? transpose(ΔΩ) : ΔΩ)

function _hermitian_back(ΔΩ, uplo)
    L, U, rD = LowerTriangular(ΔΩ), UpperTriangular(ΔΩ), real.(Diagonal(ΔΩ))
    return uplo == 'U' ? U .+ L' - rD : L .+ U' - rD
end
_hermitian_back(ΔΩ::Diagonal, uplo) = real.(ΔΩ)
function _hermitian_back(ΔΩ::LinearAlgebra.AbstractTriangular, uplo)
    ∂UL = ΔΩ .- Diagonal(_extract_imag.(diag(ΔΩ)))
    return if istriu(ΔΩ)
        return Matrix(uplo == 'U' ? ∂UL : ∂UL')
    else
        return Matrix(uplo == 'U' ? ∂UL' : ∂UL)
    end
end


text:
using HorizonSideRobots
include("BackPath.jl")
include("Structs.jl")

function count_borders!(robot::Robot)
    to_sud_west_corner = BackPath(r)
    robot = CountbordersRobot(robot)
    snake(robot)
    return get_num_borders(robot) - 1
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
module UnstructuredGeometryTests

using Test
using Gridap

c2v = [[1,2,4,5],[2,3,5,6]]
v2x = Point{2,Float64}[(0.,0.),(0.,1.),(0.,2.),(1.,0.),(1.,1.),(1.,2.)]

l = length(c2v)
order = 1
t = (HEX_AXIS,HEX_AXIS)
c2t = ConstantCellValue(t,l)
c2o = ConstantCellValue(order,l)

c2v_data, c2v_ptrs = generate_data_and_ptrs(c2v)

grid = UnstructuredGrid(v2x,c2v_data,c2v_ptrs,c2t,c2o)
@test isa(grid,UnstructuredGrid)
test_grid(grid,6,2)

grid = FlexibleUnstructuredGrid(v2x,c2v,c2t,c2o)
@test isa(grid,FlexibleUnstructuredGrid)
test_grid(grid,6,2)

end # module


text:
using Trello
using Test
using UUIDs

@testset "Trello.jl" begin

    cred = TrelloCred()
    @testset "End to End Usage" begin
        new_board_name = "ZZ_testing_$(uuid1())"

        # Should be able to create boards:
        new_board = create_board(cred, new_board_name; desc="testing")
        boards = get_boards(cred)
        @test haskey(boards, new_board_name)
        board = boards[new_board_name]
        @test board.desc == "testing"
        board_id = boards[new_board_name].id

        # the resturn from `create_board` at least has same key fields as can be read from `get_boards`
        @test new_board.id == board.id
        @test new_board.name == board.name

        # Test labels
        @test isempty(get_labels(cred, board_id))
        @test isempty(get_labels(cred, board))  # not passing ID but passing the board object

        # Should be able to create lists
        create_list(cred, board_id, "L1")
        create_list(cred, board_id, "L2")
        create_list(cred, board, "L3")  # not passing ID but passing the board object

        # should be able to get all the lists
        lists = get_lists(cred, board_id)
        @test lists == get_lists(cred, board)  # not passing ID but passing the board object

        @test collect(keys(lists)) == ["L1", "L2", "L3"]
        list = lists["L2"]
        list_id = list.id

        # Test cards
        create_card(cred, list_id, "C1")
        create_card(cred, list_id, "C2")
        create_card(cred, list, "C3")  # not passing ID but passing the list object

        cards = get_cards(cred, list_id)
        @test get_cards(cred, list) == cards   # same result wether using list_id or list object
        @test collect(keys(cards)) == ["C1", "C2", "C3"]

        # After deleting board should be gone
        delete_board(cred, board_id)
        @test !haskey(get_boards(cred), new_board_name)
    end

    @testset "delete_board by board object" begin
        new_board_name = "YY_testing_$(uuid1())"
        # Should be able to create boards:
        board = create_board(cred, new_board_name)
        @test haskey(get_boards(cred), new_board_name)
        delete_board(cred, board)
        @test !haskey(get_boards(cred), new_board_name)
    end
end


text:
#!/usr/bin/julia

function SelectionSort(arr)
    l = length(arr)
    sorted_list = []
    for i = 1:l
        m = minimum(arr)
        push!(sorted_list,m)
        deleteat!(arr, findfirst(x->x==m,arr))
    end
    return sorted_list
end

function error()
    println("Usage: please provide a list of at least two integers to sort in the format \"1, 2, 3, 4, 5\"")
end

function HandleInput(inp)

    a = split(inp,",")
    a = map(x->split(x," "),a)
    a = map(x->last(x),a)
    numbers = map(x->parse(Int,x),a)
    if length(numbers) == 1
        error()
        exit()
    end
    return numbers

end

function PrintOutput(out)
    for i = 1:length(out)
        print(out[i])
        if i != length(out)
            print(", ")
        end
    end
    println()
end

try

    n = HandleInput(ARGS[1])
    sorted = SelectionSort(n)
    PrintOutput(sorted)

catch e
    error()
end


text:
module LeetCode

const LC = LeetCode
export LC

include("common.jl")
include("problems/problems.jl")

end


text:
module TestReaderWriter

using TSML
using TSML.Utils
using TSML.TSMLTypes
using TSML.TSMLTransformers

using TSML.DataReaders
using TSML.DataWriters

using DataFrames
using Dates
using Test

function test_readerwriter()
    gcdims = (8761,2)
    ssum = 97564.0
    resdf=DataFrame()
    datapath=joinpath(dirname(pathof(TSML)),"../data")
    basefilename = "testdateval"
    fname = joinpath(datapath,basefilename*".csv")
    lcsv=DataReader(Dict(:filename=>fname))
    fit!(lcsv)
    dateval=transform!(lcsv)
    @test sum(size(dateval) .== gcdims ) == 2
    @test sum(dateval[:Value]) |> round == ssum
    csvname = replace(fname,"test"=>"out")
    wcsv = DataWriter(Dict(:filename=>csvname))
    fit!(wcsv)
    transform!(wcsv,dateval)
    pcsv = DataReader(Dict(:filename=>csvname))
    fit!(pcsv)
    resdf=transform!(pcsv)
    @test sum(size(resdf) .== gcdims) == 2
    @test sum(resdf[:Value]) |> round == ssum
    # check hdf5
    hdf5name = replace(fname,"csv"=>"h5")
    lhdf5 = DataWriter(Dict(:filename=>hdf5name))
    fit!(lhdf5)
    transform!(lhdf5,dateval)
    whdf5 = DataReader(Dict(:filename=>hdf5name))
    fit!(whdf5)
    resdf = transform!(whdf5)
    @test sum(size(resdf) .== gcdims) == 2
    @test sum(resdf[:Value]) |> round == ssum
    # check feather
    feathername = replace(fname,"csv"=>"feather")
    lfeather = DataWriter(Dict(:filename=>feathername))
    fit!(lfeather)
    transform!(lfeather,dateval)
    wfeather = DataReader(Dict(:filename=>feathername))
    fit!(wfeather)
    resdf = transform!(wfeather)
    @test sum(size(resdf) .== gcdims) == 2
    @test sum(resdf[:Value]) |> round == ssum
    # check jld
    jldname = replace(fname,"csv"=>"jld")
    ljld = DataWriter(Dict(:filename=>jldname))
    fit!(ljld)
    transform!(ljld,dateval)
    wjld = DataReader(Dict(:filename=>jldname))
    fit!(wjld)
    resdf = transform!(wjld)
    @test sum(size(resdf) .== gcdims) == 2
    @test sum(resdf[:Value]) |> round == ssum
    # cleanup
    rm(csvname,force=true)
    rm(hdf5name,force=true)
    rm(jldname,force=true)
    rm(feathername,force=true)
end
@testset "Data Readers/Writers: csv,hdf5,feather,jld" begin
    test_readerwriter()
end

end


text:
# """
# void ceres_init(){
# """
# function ceres_init()
#     @show "init ing ceres!"
#     ccall((:ceres_init, :libceres_c_api), Cvoid, ())
# end

function createproblem_c()::Ptr{Cvoid}
    return ccall((:ceres_create_problem, "libceres_c_api"), Ptr{Cvoid}, ())
end

function solve_c(problem, max_iter_num::Int64, solver_type::Int64)
    ccall( (:ceres_solve, "libceres_c_api"),
        Cvoid,
        (Ptr{Cvoid}, Cint, Cint),
        problem, convert(Cint, max_iter_num),
        convert(Cint, solver_type))
end

function freeproblem_c(problem)
    ccall((:ceres_free_problem, "libceres_c_api"), Cvoid, (Ptr{Cvoid},), problem)
end

"""
# void* ceres_problem_add_residual_block(
#    void* problem,
#    void* cost_function,
#    double* cost_function_data,
#    int cost_function_data_size,
#    void* loss_function,
#    double* loss_function_data,
#    int num_residuals,
#    int num_parameter_blocks,
#    int* parameter_block_sizes,
#    double** parameters);

lossfunction::Ptr{Cvoid},  # now is ceres c loss. It could be made to use
 cfunction in the future
"""
function addresidualblock_c(problem::Ptr{Cvoid},
                            costfunction::Ptr{Cvoid}, # cfunction
                            data::Ptr{Cdouble},
                            cost_function_data_size::Cint,
                            lossfunction::Ptr{Cvoid},  # ceres c loss
                            lossfunctionuserdata::Ptr{Cdouble},
                            numresiduals::Cint,
                            numparamblocks::Cint,
                            parametersizes::Ptr{Cint},
                            parameterpointers::Ptr{Ptr{Cdouble}})

    # void* ceres_problem_add_residual_block(
    #    void* problem,
    #    void* cost_function,
    #    double* cost_function_data,
    #    int cost_function_data_size,
    #    void* loss_function,
    #    double* loss_function_data,
    #    int num_residuals,
    #    int num_parameter_blocks,
    #    int* parameter_block_sizes,
    #    double** parameters);
    val = ccall( (:ceres_problem_add_residual_block, "libceres_c_api"),
              Ptr{Cvoid},
              (Ptr{Cvoid},            # ceres_problem_t* problem,
               Ptr{Cvoid},            # ceres_cost_function_t cost_function,
               Ptr{Cdouble},          # void* cost_function_data,
               Cint,                  # int cost_function_data_size,
               Ptr{Cvoid},            # ceres_loss_function_t loss_function,
               Ptr{Cdouble},          # void* loss_function_data,
               Cint,                  # int num_residuals,
               Cint,                  # int num_parameter_blocks,
               Ptr{Cint},             # int* parameter_block_sizes,
               Ptr{Ptr{Cdouble}}),    # double** parameters
              problem,
              costfunction,
              data,
              cost_function_data_size,
              lossfunction,
              lossfunctionuserdata,
              numresiduals,
              numparamblocks,
              parametersizes,
              parameterpointers )
    return val
end

function set_local_parameterization_c(problem::Ptr{Cvoid}, # why uint8? see line 6
                                      parameters::Ptr{Cdouble},
                                      local_parameterization::Ptr{Cvoid})
    ccall( (:ceres_problem_set_parameterization, "libceres_c_api"),
                Cvoid,
                (Ptr{UInt8},            # ceres_problem_t* problem,
                 Ptr{Cdouble},            # ceres_cost_function_t cost_function,
                 Ptr{Cvoid}),          # void* cost_function_data,
                problem,
                parameters,
                local_parameterization)

end

"""
void* ceres_create_quaternion_parameterization();
"""
function create_quaternion_parameterization_c()
    ret = ccall( (:ceres_create_quaternion_parameterization, "libceres_c_api"),
               Ptr{Cvoid}, () )
    return ret
end

# void ceres_set_parameter_block_constant(void* problem, double* data){
function ceres_set_parameter_block_constant_c(problem::Ptr{Cvoid}, # why uint8? see line 6
                                              parameters::Ptr{Cdouble})
    ccall((:ceres_set_parameter_block_constant, "libceres_c_api"),
          Cvoid, (Ptr{Cvoid},            # ceres_problem_t* problem,
                 Ptr{Cdouble}),
           problem,
           parameters)
end

# void ceres_set_parameter_block_variable(void* c_problem, double *values) {
function ceres_set_parameter_block_variable_c(problem::Ptr{Cvoid}, # why uint8? see line 6
                                              parameters::Ptr{Cdouble})
    ccall((:ceres_set_parameter_block_variable, "libceres_c_api"),
          Cvoid, (Ptr{Cvoid},            # ceres_problem_t* problem,
                 Ptr{Cdouble}),
           problem,
           parameters)
end

# void* ceres_create_huber_loss(double delta){
function ceres_create_huber_loss_c(δ::Float64)::Ptr{Cvoid}
    ccall((:ceres_create_huber_loss, "libceres_c_api"),
           Ptr{Cvoid}, (Cdouble,), δ)
end

# // bounds
#
# void ceres_SetParameterLowerBound(void* c_problem, double *values, int index, double lower_bound) {
#     ceres::Problem* problem = reinterpret_cast<ceres::Problem*>(c_problem);
#     problem->SetParameterLowerBound(values, index, lower_bound);
# }

function ceres_set_parameter_lower_bound_c(
      problem::Ptr{Cvoid}, # why uint8? see line 6
      parameters::Ptr{Cdouble},
      index_jl::Union{Int32,Int64},
      lower_bound::Float64
)
    ccall((:ceres_SetParameterLowerBound, "libceres_c_api"),
          Cvoid, (Ptr{Cvoid},            # ceres_problem_t* problem,
                  Ptr{Cdouble},
                  Cint,
                  Cdouble),
           problem,
           parameters,
           Cint(index_jl - 1),
           Cdouble(lower_bound))
end

# void ceres_SetParameterUpperBound(void* c_problem, double *values, int index, double upper_bound) {
#      ceres::Problem* problem = reinterpret_cast<ceres::Problem*>(c_problem);
#      problem->SetParameterUpperBound(values, index, upper_bound);
# }
function ceres_set_parameter_upper_bound_c(
      problem::Ptr{Cvoid},
      parameters::Ptr{Cdouble},
      index_jl::Union{Int32,Int64},
      upper_bound::Float64
)
    ccall((:ceres_SetParameterUpperBound, "libceres_c_api"),
          Cvoid, (Ptr{Cvoid},            # ceres_problem* problem,
                  Ptr{Cdouble},
                  Cint,
                  Cdouble),
           problem,
           parameters,
           Cint(index_jl - 1),
           Cdouble(upper_bound))
end

# // constant
# void ceres_SetParameterBlockConstant(void* c_problem, double *values) {
#      ceres::Problem* problem = reinterpret_cast<ceres::Problem*>(c_problem);
#      problem->SetParameterBlockConstant(values);
# }
# void ceres_SetParameterBlockVariable(void* c_problem, double *values) {
#      ceres::Problem* problem = reinterpret_cast<ceres::Problem*>(c_problem);
#      problem->SetParameterBlockVariable(values);
# }


text:
function displacement(alpha, x)
    u1 = alpha * x[2] * sin(pi * x[1])
    u2 = alpha * (x[1]^3 + cos(pi * x[2]))
    return [u1, u2]
end

function stress_field(lambda, mu, alpha, x)
    s11 =
        (lambda + 2mu) * alpha * pi * x[2] * cos(pi * x[1]) -
        lambda * alpha * pi * sin(pi * x[2])
    s22 =
        -(lambda + 2mu) * alpha * pi * sin(pi * x[2]) +
        lambda * alpha * pi * x[2] * cos(pi * x[1])
    s12 = alpha * mu * (3x[1]^2 + sin(pi * x[1]))
    return [s11, s22, s12]
end

function body_force(lambda, mu, alpha, x)
    b1 = alpha * (lambda + 2mu) * pi^2 * x[2] * sin(pi * x[1])
    b2 =
        -alpha * (6mu * x[1] + (lambda + mu) * pi * cos(pi * x[1])) +
        alpha * (lambda + 2mu) * pi^2 * cos(pi * x[2])
    return [b1, b2]
end

function onleftboundary(x, L, W)
    return x[1] ≈ 0.0
end

function onbottomboundary(x, L, W)
    return x[2] ≈ 0.0
end

function onrightboundary(x, L, W)
    return x[1] ≈ L
end

function ontopboundary(x, L, W)
    return x[2] ≈ W
end

function onboundary(x, L, W)
    return onbottomboundary(x, L, W) ||
           onrightboundary(x, L, W) ||
           ontopboundary(x, L, W) ||
           onleftboundary(x, L, W)
end

function construct_mesh_and_quadratures(
    meshwidth,
    nelmts,
    elasticitybasis,
    distancefunction,
    numqp,
)
    levelsetbasis = HermiteTensorProductBasis(2)

    dim, nf = size(interpolation_points(levelsetbasis))
    refpoints = interpolation_points(elasticitybasis)

    mesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], refpoints)
    cgmesh = CutCellDG.CGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], nf)

    levelset = CutCellDG.LevelSet(distancefunction, cgmesh, levelsetbasis)

    cutmesh = CutCellDG.CutMesh(mesh, levelset)
    cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
    interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
    facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)

    mergedmesh =
        CutCellDG.MergedMesh!(cutmesh, cellquads, facequads, interfacequads)

    return mergedmesh, cellquads, facequads, interfacequads, levelset
end


text:
using Revise
using Soss

m = @model μ begin
    x ~ Normal(μ,1)
end

m(μ=1.0)

rand(m(μ=1.0))

logdensityof(m(μ=1.0),(x=0.4,))
# rand(m1)

weightedSample(m(μ=1.0),(x=3,))

m = @model μ begin
    x ~ Normal(μ,1)
end

mμ = m(μ=1)

@btime weightedSample($mμ,(x=3,))
@btime weightedSample($mμ,NamedTuple())