text:
"""
    mutable struct MatrixOfConstraints{T,AT,BT,ST} <: MOI.ModelLike
        coefficients::AT
        constants::BT
        sets::ST
        caches::Vector{Any}
        are_indices_mapped::Vector{BitSet}
        final_touch::Bool
    end

Represent `ScalarAffineFunction` and `VectorAffinefunction` constraints in a
matrix form where the linear coefficients of the functions are stored in the
`coefficients` field, the constants of the functions or sets are stored in the
`constants` field. Additional information about the sets are stored in the
`sets` field.

This model can only be used as the `constraints` field of a
`MOI.Utilities.AbstractModel`.

When the constraints are added, they are stored in the `caches` field. They are
only loaded in the `coefficients` and `constants` fields once
`MOI.Utilities.final_touch` is called. For this reason, `MatrixOfConstraints`
should not be used by an incremental interface. Use `MOI.copy_to` instead.

The constraints can be added in two different ways:

1) With `add_constraint`, in which case a canonicalized copy of the function is
   stored in `caches`.
2) With `pass_nonvariable_constraints`, in which case the functions and sets are
   stored themselves in `caches` without mapping the variable indices. The
   corresponding index in `caches` is added in `are_indices_mapped`. This avoids
   doing a copy of the function in case the getter of
   `CanonicalConstraintFunction` does not make a copy for the source model,
   e.g., this is the case of `VectorOfConstraints`.

We illustrate this with an example. Suppose a model is copied from a
`src::MOI.Utilities.Model` to a bridged model with a `MatrixOfConstraints`. For
all the types that are not bridged, the constraints will be copied with
`pass_nonvariable_constraints`. Hence the functions stored in `caches` are
exactly the same as the ones stored in `src`. This is ok since this is only
during the `copy_to` operation during which `src` cannot be modified. On the
other hand, for the types that are bridged, the functions added may contain
duplicates even if the functions did not contain duplicates in `src` so
duplicates are removed with `MOI.Utilities.canonical`.

## Interface

The `.coefficients::AT` type must implement:
 * `AT()`
 * `MOI.empty(::AT)!`
 * [`MOI.Utilities.add_column`](@ref)
 * [`MOI.Utilities.set_number_of_rows`](@ref)
 * [`MOI.Utilities.allocate_terms`](@ref)
 * [`MOI.Utilities.load_terms`](@ref)
 * [`MOI.Utilities.final_touch`](@ref)

The `.constants::BT` type must implement:

 * `BT()`
 * `Base.empty!(::BT)`
 * `Base.resize(::BT)`
 * [`MOI.Utilities.load_constants`](@ref)
 * [`MOI.Utilities.function_constants`](@ref)
 * [`MOI.Utilities.set_from_constants`](@ref)

The `.sets::ST` type must implement:

 * `ST()`
 * `MOI.is_empty(::ST)`
 * `MOI.empty(::ST)`
 * `MOI.dimension(::ST)`
 * `MOI.is_valid(::ST, ::MOI.ConstraintIndex)`
 * `MOI.get(::ST, ::MOI.ListOfConstraintTypesPresent)`
 * `MOI.get(::ST, ::MOI.NumberOfConstraints)`
 * `MOI.get(::ST, ::MOI.ListOfConstraintIndices)`
 * [`MOI.Utilities.set_types`](@ref)
 * [`MOI.Utilities.set_index`](@ref)
 * [`MOI.Utilities.add_set`](@ref)
 * [`MOI.Utilities.rows`](@ref)
 * [`MOI.Utilities.final_touch`](@ref)
"""
mutable struct MatrixOfConstraints{T,AT,BT,ST} <: MOI.ModelLike
    coefficients::AT
    constants::BT
    sets::ST
    caches::Vector{Any}
    are_indices_mapped::Vector{BitSet}
    final_touch::Bool
    function MatrixOfConstraints{T,AT,BT,ST}() where {T,AT,BT,ST}
        model = new{T,AT,BT,ST}(AT(), BT(), ST(), Any[], BitSet[], false)
        MOI.empty!(model)
        return model
    end
end

###
### Interface for the .coefficients field
###

"""
    add_column(coefficients)::Nothing

Tell `coefficients` to pre-allocate datastructures as needed to store one
column.
"""
function add_column end

"""
    set_number_of_rows(coefficients, n)::Nothing

Tell `coefficients` to pre-allocate datastructures as needed to store `n` rows.
"""
function set_number_of_rows end

"""
    allocate_terms(coefficients, index_map, func)::Nothing

Tell `coefficients` that the terms of the function `func` where the variable
indices are mapped with `index_map` will be loaded with [`load_terms`](@ref).

The function `func` must be canonicalized before calling `allocate_terms`. See
[`is_canonical`](@ref).
"""
function allocate_terms end

"""
    load_terms(coefficients, index_map, func, offset)::Nothing

Loads the terms of `func` to `coefficients`, mapping the variable indices with
`index_map`.

The `i`th dimension of `func` is loaded at the `(offset + i)`th row of
`coefficients`.

The function must be allocated first with [`allocate_terms`](@ref).

The function `func` must be canonicalized, see [`is_canonical`](@ref).
"""
function load_terms end

"""
    final_touch(coefficients)::Nothing

Informs the `coefficients` that all functions have been added with `load_terms`.
No more modification is allowed unless `MOI.empty!` is called.

    final_touch(sets)::Nothing

Informs the `sets` that all functions have been added with `add_set`.
No more modification is allowed unless `MOI.empty!` is called.
"""
function final_touch end

"""
    extract_function(coefficients, row::Integer, constant::T) where {T}

Return the `MOI.ScalarAffineFunction{T}` function corresponding to row `row` in
`coefficients`.

    extract_function(
        coefficients,
        rows::UnitRange,
        constants::Vector{T},
    ) where{T}

Return the `MOI.VectorAffineFunction{T}` function corresponding to rows `rows`
in `coefficients`.
"""
function extract_function end

###
### Interface for the .constants field
###

"""
    load_constants(constants, offset, func_or_set)::Nothing

This function loads the constants of `func_or_set` in `constants` at an offset
of `offset`. Where `offset` is the sum of the dimensions of the constraints
already loaded. The storage should be preallocated with `resize!` before calling
this function.

This function should be implemented to be usable as storage of constants for
[`MatrixOfConstraints`](@ref).

The constants are loaded in three steps:
 1) `Base.empty!` is called.
 2) `Base.resize!` is called with the sum of the dimensions of all constraints.
 3) `MOI.Utilities.load_constants` is called for each function for vector
    constraint or set for scalar constraint.
"""
function load_constants end

"""
    function_constants(constants, rows)

This function returns the function constants that were loaded with
[`load_constants`](@ref) at the rows `rows`.

This function should be implemented to be usable as storage of constants for
[`MatrixOfConstraints`](@ref).
"""
function function_constants end

"""
    set_from_constants(constants, S::Type, rows)::S

This function returns an instance of the set `S` for which the constants where
loaded with [`load_constants`](@ref) at the rows `rows`.

This function should be implemented to be usable as storage of constants for
[`MatrixOfConstraints`](@ref).
"""
function set_from_constants end

###
### Interface for the .sets field
###

"""
    set_types(sets)::Vector{Type}

Return the list of the types of the sets allowed in `sets`.
"""
function set_types end

"""
    set_index(sets, ::Type{S})::Union{Int,Nothing} where {S<:MOI.AbstractSet}

Return an integer corresponding to the index of the set type in the list given
by [`set_types`](@ref).

If `S` is not part of the list, return `nothing`.
"""
function set_index end

"""
    add_set(sets, i)::Int64

Add a scalar set of type index `i`.

    add_set(sets, i, dim)::Int64

Add a vector set of type index `i` and dimension `dim`.

Both methods return a unique `Int64` of the set that can be used to reference
this set.
"""
function add_set end

"""
    rows(sets, ci::MOI.ConstraintIndex)::Union{Int,UnitRange{Int}}

Return the rows in `1:MOI.dimension(sets)` corresponding to the set of id
`ci.value`.

For scalar sets, this returns an `Int`. For vector sets, this returns an
`UnitRange{Int}`.
"""
function rows end

###
### MatrixOfConstraints
###

MOI.is_empty(v::MatrixOfConstraints) = MOI.is_empty(v.sets)

function MOI.empty!(v::MatrixOfConstraints{T}) where {T}
    MOI.empty!(v.coefficients)
    empty!(v.constants)
    MOI.empty!(v.sets)
    v.caches =
        [Tuple{_affine_function_type(T, S),S}[] for S in set_types(v.sets)]
    v.are_indices_mapped = [BitSet() for _ in eachindex(v.caches)]
    v.final_touch = false
    return
end

"""
    rows(model::MatrixOfConstraints, ci::MOI.ConstraintIndex)

Return the rows in `1:MOI.dimension(sets)` corresponding to the set of id
`ci.value`.

For scalar sets, this returns an `Int`. For vector sets, this returns an
`UnitRange{Int}`.
"""
rows(model::MatrixOfConstraints, ci::MOI.ConstraintIndex) = rows(model.sets, ci)

function _affine_function_type(
    ::Type{T},
    ::Type{<:MOI.AbstractScalarSet},
) where {T}
    return MOI.ScalarAffineFunction{T}
end

function _affine_function_type(
    ::Type{T},
    ::Type{<:MOI.AbstractVectorSet},
) where {T}
    return MOI.VectorAffineFunction{T}
end

function MOI.supports_constraint(
    v::MatrixOfConstraints{T},
    ::Type{F},
    ::Type{S},
) where {T,F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    return F == _affine_function_type(T, S) && set_index(v.sets, S) !== nothing
end

function MOI.is_valid(
    v::MatrixOfConstraints{T},
    ci::MOI.ConstraintIndex{F,S},
) where {T,F,S}
    return F == _affine_function_type(T, S) && MOI.is_valid(v.sets, ci)
end

function MOI.get(
    v::MatrixOfConstraints,
    attr::Union{
        MOI.ListOfConstraintTypesPresent,
        MOI.NumberOfConstraints,
        MOI.ListOfConstraintIndices,
    },
)
    return MOI.get(v.sets, attr)
end

_add_set(sets, i, ::MOI.AbstractScalarFunction) = add_set(sets, i)

function _add_set(sets, i, func::MOI.AbstractVectorFunction)
    return add_set(sets, i, MOI.output_dimension(func))
end

const _MATRIXOFCONSTRAINTS_MODIFY_NOT_ALLOWED_ERROR_MESSAGE = """
MatrixOfConstraints does not allow modifications to be made to the model once
`MOI.Utilities.final_touch` has been called. This is called at the end of
`MOI.copy_to` and in `MOI.Utilities.attach_optimizer` (which is called by
`MOI.optimize!` in a `MOI.Utilities.CachingOptimizer`). In order to be able to
apply modifications to this model, you should add a layer
`MOI.Utilities.CachingOptimizer(MOI.Utilities.Model{Float64}(), model)`
where `model` is the current model. This will automatically empty `model` when
modifications are done after `MOI.Utilities.final_touch` is called and copy the
model again in `MOI.Utilities.attach_optimizer`.
"""

function _add_constraint(
    model::MatrixOfConstraints,
    i::Int,
    index_map,
    func::F,
    set::S,
) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    allocate_terms(model.coefficients, index_map, func)
    # Without this type annotation, the compiler is unable to know the type
    # of `caches[i]` so this is slower and produce an allocation.
    push!(model.caches[i]::Vector{Tuple{F,S}}, (func, set))
    return MOI.ConstraintIndex{F,S}(_add_set(model.sets, i, func))
end

struct IdentityMap <: AbstractDict{MOI.VariableIndex,MOI.VariableIndex} end

Base.getindex(::IdentityMap, vi::MOI.VariableIndex) = vi

function MOI.add_constraint(
    model::MatrixOfConstraints{T},
    func::F,
    set::S,
) where {T,F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    i = set_index(model.sets, S)
    if i === nothing || F != _affine_function_type(T, S)
        throw(MOI.UnsupportedConstraint{F,S}())
    end
    if model.final_touch
        throw(
            MOI.AddConstraintNotAllowed{F,S}(
                _MATRIXOFCONSTRAINTS_MODIFY_NOT_ALLOWED_ERROR_MESSAGE,
            ),
        )
    end
    return _add_constraint(model, i, IdentityMap(), func, set)
end

function _allocate_constraints(
    model::MatrixOfConstraints{T},
    src,
    index_map,
    ::Type{F},
    ::Type{S},
) where {T,F,S}
    i = set_index(model.sets, S)
    if i === nothing || F != _affine_function_type(T, S)
        throw(MOI.UnsupportedConstraint{F,S}())
    end
    cis_src = MOI.get(
        src,
        MOI.ListOfConstraintIndices{_affine_function_type(T, S),S}(),
    )
    for ci_src in cis_src
        func = MOI.get(src, MOI.CanonicalConstraintFunction(), ci_src)
        set = MOI.get(src, MOI.ConstraintSet(), ci_src)
        push!(model.are_indices_mapped[i], length(model.caches[i]) + 1)
        index_map[ci_src] = _add_constraint(model, i, index_map, func, set)
    end
    return
end

function _load_constants(
    constants,
    offset,
    func::MOI.AbstractScalarFunction,
    set::MOI.AbstractScalarSet,
)
    MOI.throw_if_scalar_and_constant_not_zero(func, typeof(set))
    return load_constants(constants, offset, set)
end

function _load_constants(
    constants,
    offset,
    func::MOI.AbstractVectorFunction,
    ::MOI.AbstractVectorSet,
)
    return load_constants(constants, offset, func)
end

function _load_constraints(
    dest::MatrixOfConstraints,
    index_map,
    offset,
    func_sets,
    are_indices_mapped,
)
    for i in eachindex(func_sets)
        func, set = func_sets[i]
        if i in are_indices_mapped
            load_terms(dest.coefficients, index_map, func, offset)
        else
            load_terms(dest.coefficients, IdentityMap(), func, offset)
        end
        _load_constants(dest.constants, offset, func, set)
        offset += MOI.output_dimension(func)
    end
    return offset
end

_add_variable(model::MatrixOfConstraints) = add_column(model.coefficients)

function _add_variables(model::MatrixOfConstraints, n)
    return add_columns(model.coefficients, n)
end

function pass_nonvariable_constraints(
    dest::MatrixOfConstraints,
    src::MOI.ModelLike,
    index_map::IndexMap,
    constraint_types,
)
    for (F, S) in constraint_types
        _allocate_constraints(dest, src, index_map, F, S)
    end
    return
end

function final_touch(model::MatrixOfConstraints, index_map)
    if model.final_touch
        # If `default_copy_to` calls `final_touch`, then `index_map` is not
        # `nothing` and `final_touch` should be `false`.
        # When `CachingOptimizer` calls this, `index_map` is `nothing`
        # and `final_touch` might be `true`.
        # So we are always in a case where `index_map` is `nothing`.
        @assert index_map === nothing
        return
    end
    final_touch(model.sets)
    num_rows = MOI.dimension(model.sets)
    resize!(model.constants, num_rows)
    set_number_of_rows(model.coefficients, num_rows)
    offset = 0
    for (cache, mapped_indices) in zip(model.caches, model.are_indices_mapped)
        offset =
            _load_constraints(model, index_map, offset, cache, mapped_indices)
    end
    final_touch(model.coefficients)
    empty!(model.caches)
    empty!(model.are_indices_mapped)
    model.final_touch = true
    return
end

# Users of `MatrixOfConstraints` assume variable indices to be `1:n` where `n`
# is the number of columns so we don't support variable deletion.
function _throw_if_cannot_delete(
    ::MatrixOfConstraints,
    ::Vector{MOI.VariableIndex},
    vi::MOI.VariableIndex,
)
    return throw(MOI.DeleteNotAllowed(vi))
end

function _throw_if_cannot_delete(
    ::MatrixOfConstraints,
    vis::Vector{MOI.VariableIndex},
    ::Set{MOI.VariableIndex},
)
    return throw(MOI.DeleteNotAllowed(first(vis)))
end

###
### .constants::Vector
###

# Base.empty! is already implemented.
# Base.resize! is already implemented.

function load_constants(
    b::Vector{T},
    offset,
    func::MOI.VectorAffineFunction{T},
) where {T}
    copyto!(b, offset + 1, func.constants)
    return
end

function_constants(b::Vector, rows) = b[rows]

"""
    set_with_dimension(::Type{S}, dim) where {S<:MOI.AbstractVectorSet}

Returns the instance of `S` of [`VecMathOptInterface.dimension`](@ref) `dim`.
This needs to be implemented for sets of type `S` to be useable with
[`MatrixOfConstraints`](@ref).
"""
function set_with_dimension(::Type{S}, dim) where {S<:MOI.AbstractVectorSet}
    return S(dim)
end

function set_with_dimension(
    ::Type{S},
    dim,
) where {S<:MOI.AbstractSymmetricMatrixSetTriangle}
    side_dimension = side_dimension_for_vectorized_dimension(dim)
    return S(side_dimension)
end

function set_with_dimension(
    ::Type{S},
    dim,
) where {S<:MOI.AbstractSymmetricMatrixSetSquare}
    return S(isqrt(dim))
end

function set_with_dimension(::Type{MOI.LogDetConeTriangle}, dim)
    side_dimension = side_dimension_for_vectorized_dimension(dim - 2)
    return MOI.LogDetConeTriangle(side_dimension)
end

function set_with_dimension(::Type{MOI.LogDetConeSquare}, dim)
    return MOI.LogDetConeSquare(isqrt(dim - 2))
end

function set_with_dimension(::Type{MOI.RootDetConeTriangle}, dim)
    side_dimension = side_dimension_for_vectorized_dimension(dim - 1)
    return MOI.RootDetConeTriangle(side_dimension)
end

function set_with_dimension(::Type{MOI.RootDetConeSquare}, dim)
    return MOI.RootDetConeSquare(isqrt(dim - 1))
end

function set_from_constants(::Vector, ::Type{S}, rows) where {S}
    return set_with_dimension(S, length(rows))
end

function MOI.get(
    model::MatrixOfConstraints,
    ::Union{MOI.CanonicalConstraintFunction,MOI.ConstraintFunction},
    ci::MOI.ConstraintIndex,
)
    @assert model.final_touch
    MOI.throw_if_not_valid(model, ci)
    r = rows(model, ci)
    return extract_function(
        model.coefficients,
        r,
        function_constants(model.constants, r),
    )
end

function MOI.get(
    model::MatrixOfConstraints,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    @assert model.final_touch
    MOI.throw_if_not_valid(model, ci)
    return set_from_constants(model.constants, S, rows(model, ci))
end


text:
"""
Full positive definite matrix together with a Cholesky factorization object.
"""
struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
    dim::Int
    mat::S
    chol::CholType{T,S}

    PDMat{T,S}(d::Int,m::AbstractMatrix{T},c::CholType{T,S}) where {T,S} = new{T,S}(d,m,c)
end

function PDMat(mat::AbstractMatrix,chol::CholType)
    d = size(mat, 1)
    size(chol, 1) == d ||
        throw(DimensionMismatch("Dimensions of mat and chol are inconsistent."))
    PDMat{eltype(mat),typeof(mat)}(d, mat, chol)
end

PDMat(mat::Matrix) = PDMat(mat, cholesky(mat))
PDMat(mat::Symmetric) = PDMat(Matrix(mat))
PDMat(fac::CholType) = PDMat(Matrix(fac), fac)

### Conversion
Base.convert(::Type{PDMat{T}},         a::PDMat) where {T<:Real} = PDMat(convert(AbstractArray{T}, a.mat))
Base.convert(::Type{AbstractArray{T}}, a::PDMat) where {T<:Real} = convert(PDMat{T}, a)
Base.convert(::Type{AbstractArray{T}}, a::PDMat{T}) where {T<:Real} = a

### Basics

dim(a::PDMat) = a.dim
Base.Matrix(a::PDMat) = copy(a.mat)
LinearAlgebra.diag(a::PDMat) = diag(a.mat)
LinearAlgebra.cholesky(a::PDMat) = a.chol

### Arithmetics

function pdadd!(r::Matrix, a::Matrix, b::PDMat, c)
    @check_argdims size(r) == size(a) == size(b)
    _addscal!(r, a, b.mat, c)
end

*(a::PDMat{S}, c::T) where {S<:Real, T<:Real} = PDMat(a.mat * c)
*(a::PDMat, x::AbstractVecOrMat) = a.mat * x
\(a::PDMat, x::AbstractVecOrMat) = a.chol \ x

### Algebra

Base.inv(a::PDMat) = PDMat(inv(a.chol))
LinearAlgebra.logdet(a::PDMat) = logdet(a.chol)
LinearAlgebra.eigmax(a::PDMat) = eigmax(a.mat)
LinearAlgebra.eigmin(a::PDMat) = eigmin(a.mat)
Base.kron(A::PDMat, B::PDMat) = PDMat(kron(A.mat, B.mat), Cholesky(kron(A.chol.U, B.chol.U), 'U', A.chol.info))

### whiten and unwhiten

function whiten!(r::StridedVecOrMat, a::PDMat, x::StridedVecOrMat)
    cf = a.chol.U
    v = _rcopy!(r, x)
    ldiv!(transpose(cf), v)
end

function unwhiten!(r::StridedVecOrMat, a::PDMat, x::StridedVecOrMat)
    cf = a.chol.U
    v = _rcopy!(r, x)
    lmul!(transpose(cf), v)
end

### quadratic forms

quad(a::PDMat, x::AbstractVector) = sum(abs2, a.chol.U * x)
invquad(a::PDMat, x::AbstractVector) = sum(abs2, a.chol.L \ x)

"""
    quad!(r::AbstractArray, a::AbstractPDMat, x::StridedMatrix)

Overwrite `r` with the value of the quadratic form defined by `a` applied columnwise to `x`
"""
quad!(r::AbstractArray, a::PDMat, x::StridedMatrix) = colwise_dot!(r, x, a.mat * x)

"""
    invquad!(r::AbstractArray, a::AbstractPDMat, x::StridedMatrix)

Overwrite `r` with the value of the quadratic form defined by `inv(a)` applied columnwise to `x`
"""
invquad!(r::AbstractArray, a::PDMat, x::StridedMatrix) = colwise_dot!(r, x, a.mat \ x)

### tri products

function X_A_Xt(a::PDMat, x::StridedMatrix)
    cf = a.chol.U
    z = rmul!(copy(x), transpose(cf))
    return z * transpose(z)
end

function Xt_A_X(a::PDMat, x::StridedMatrix)
    cf = a.chol.U
    z = lmul!(cf, copy(x))
    return transpose(z) * z
end

function X_invA_Xt(a::PDMat, x::StridedMatrix)
    cf = a.chol.U
    z = rdiv!(copy(x), cf)
    return z * transpose(z)
end

function Xt_invA_X(a::PDMat, x::StridedMatrix)
    cf = a.chol.U
    z = ldiv!(transpose(cf), copy(x))
    return transpose(z) * z
end

ext: jl


text:
"""
    compensator(ts::Number, events::Array{<:Number}, bg::Number, kappa::Number, kern::Distributions.Distribution)

The compensator or integrated intensity function.

# Arguments

* `ts::Number` time to which evaluate the integrated intensity function.
* `events::Array{<:Number}` events from the intensity function.
* `bg::Number` Background rate.
* `kappa::Number` Kappa parameter
* `kern::Distributions.Distribution` Distribution of the kernel.

# Examples

"""
function compensator(ts::Number, events::Array{<:Number}, bg::Number, kappa::Number, kern::Distributions.Distribution)
    bg * ts + kappa * sum(Distributions.cdf.(kern, ts .- events[events .< ts]))
end


text:
"""
    Webscale(m::Int)

Online k-means algorithm based on [1].

[1] Sculley, D. Web-scale k-means clustering. in Proceedings of the 19th international conference on World wide web - WWW ’10 1177 (ACM Press, 2010). doi:10.1145/1772690.1772862.
"""
struct Webscale <: OnIPSA
    m::Int
    v::Vector{Int}
end

function Webscale(m::Int)
    return Webscale(m, zeros(Int, m))
end

Base.show(io::IO, alg::Webscale) = print(io, "Webscale (m = ", alg.m, ")")

function initZ(
    rng::AbstractRNG, alg::Webscale, X::AbstractVector; arraytype=Vector{Float64}, kwargs...
)
    length(X) >= alg.m || error(
        "Input data not big enough given desired number of inducing points : $(alg.m)"
    )
    Z = to_vec_of_vecs(X[sample(rng, 1:length(X), alg.m; replace=false)], arraytype)
    return Z
end

function add_point!(
    ::AbstractRNG, Z::AbstractVector{T}, alg::Webscale, X::AbstractVector; kwargs...
) where {T}
    d = zeros(Int, length(X))
    for i in 1:length(X)
        d[i] = find_nearest_center(X[i], Z)[1] # Save the closest IP index for each point
    end
    for i in 1:length(X)
        alg.v[d[i]] += 1
        η = 1 / alg.v[d[i]]
        Z[d[i]] = (1 - η) * Z[d[i]] + η * X[i] # Update the IP position
    end
    return Z
end

function add_point(
    ::AbstractRNG, Z::AbstractVector{T}, alg::Webscale, X::AbstractVector; kwargs...
) where {T}
    d = zeros(Int, length(X))
    Z = copy(Z)
    for i in 1:length(X)
        d[i] = find_nearest_center(X[i], Z)[1] # Save the closest IP index for each point
    end
    for i in 1:length(X)
        alg.v[d[i]] += 1
        η = 1 / alg.v[d[i]]
        Z[d[i]] = (1 - η) * Z[d[i]] + η * X[i] # Update the IP position
    end
    return Z
end


text:
##===----------------------------------------------------------------------===##
##                                     BUILD.JL                               ##
##                                                                            ##
##===----------------------------------------------------------------------===##
##                                                                            ##
##   This file builds Vectorize.jl during package installation and runs the   ##
##           benchmarking suite to determine the best implementations         ##
##                                                                            ##
##===----------------------------------------------------------------------===##

# Clean up directory and build status
include("clean.jl")

# compute package directory location
currdir = @__FILE__
pkgdir = currdir[1:end-13]
@static if is_windows()
    function_location = pkgdir*"src\Functions.jl"
else
    function_location = pkgdir*"src/Functions.jl"
end

# detect architecture
if Sys.ARCH != :x86_64
    error("Vectorize.jl currently only supports x86_64; your detected architecture is $(Sys.ARCH)")
    exit(1)
end

@static if is_windows()
    include("windows.jl")
else
    include("unix.jl")
end

# Have to import vectorize after Yeppp is downloaded
@static if is_windows()
    push!(LOAD_PATH, parsedir("$(pkgdir)src\\"))
else
    push!(LOAD_PATH, parsedir("$(pkgdir)src/"))
end

# import vectorize in order to benchmark
using Vectorize

# available functions - Yeppp, VML, Accelerate register against this dictionary
functions = Vectorize.functions

# Run the benchmarking process
N = 1_000 # length of each vector
@static if is_windows()
    file = open("$(pkgdir)src\Functions.jl", "a")
else
    file = open("$(pkgdir)src/Functions.jl", "a")
end

for ((f, T), options) in functions
    if length(T) == 1
        benchmarkSingleArgFunction(f, options, T[1], file, 1_000)
    elseif length(T) == 2
        benchmarkTwoArgFunction(f, options, T, file, 1_000)
    elseif length(T) == 3
        benchmarkThreeArgFunction(f, options, T, file, 1_000)
    end
end
close(file)

Base.compilecache("Vectorize")


text:
const out_dir   = "../../docs/src/integrators/runge_kutta/lotka_volterra_2d/"
const run_id    = "lotka_volterra_2d_firk"

# create list of figures to put in list
figures_custom = (,)

include("lotka_volterra_2d_settings.jl")

# solver settings
# set_config(:nls_atol, 8eps())
# set_config(:nls_rtol, 2eps())
# set_config(:nls_stol, 2eps())
set_config(:nls_nmax, 20)


text:
using DashBootstrapComponents, DashHtmlComponents

popover = html_div([
    dbc_button(
        "Toggle",
        id="toggle",
        color="success",
        className="mr-4",
        n_clicks=0,
    ),
    dbc_button("Target", id="target", color="danger", n_clicks=0),
    dbc_popover(
        [
            dbc_popoverheader("Popover header"),
            dbc_popoverbody("And here's some amazing content. Cool!"),
        ],
        id="popover",
        is_open=false,
        target="target",
    ),
]);

callback!(
    app,
    Output("popover", "is_open"),
    Input("toggle", "n_clicks"),
    State("popover", "is_open"),
) do n, is_open
    return n > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using Gadfly
using Colors
using DataFrames
using Query

Gadfly.push_theme(Theme(major_label_font="Droid Sans",
                        minor_label_font="Droid Sans",
                        major_label_font_size=18pt, minor_label_font_size=16pt,
                        line_width=0.8mm, key_label_font="Droid Sans",
                        lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2),
                        key_label_font_size=14pt, key_position=:right))

colors = [colorant"#1f78b4", colorant"#33a02c", colorant"#e31a1c",
          colorant"#ff7f00", colorant"#6a3d9a", colorant"#a6cee3",
          colorant"#b2df8a", colorant"#fb9a99", colorant"#fdbf6f",
          colorant"#cab2d6"]

function std_or_z(x)
    s = std(x)
    if isnan(s)
        s = 0.0
    end
    s
end

function get_res_slgym(filename::String)
    names = [:ltype, :seed, :timestamp, :Problem, :Generation, :Evaluations,
             :Steps, :Reward]
    res = readtable(filename, separator=',', header=false, names=names)
    res[:Reward][(res[:Problem] .== "Acrobot-v1").&(res[:Steps] .== 0)] .-= 500
    res
end

function get_gens(res::DataFrame)
    maxs = by(res, [:Problem, :Steps, :Generation, :seed],
              df->DataFrame(Reward=maximum(df[:Reward])))
    gens = by(maxs, [:Problem, :Steps, :Generation],
              df->DataFrame(Rmean=mean(df[:Reward]), Rstd=std_or_z(df[:Reward])))
    gens[:Rlow] = gens[:Rmean] .- 0.5*gens[:Rstd]
    gens[:Rhigh] = gens[:Rmean] .+ 0.5*gens[:Rstd]
    gens
end

function plot_gym_gens(gens::DataFrame, problem::String="Acrobot-v1",
                       outfile::String="acrobot.pdf";
                       ylabel="Reward", xlabel="Generation",
                       key_position=:right)

    pgens = @from i in gens begin
        @where i.Problem == problem
        @select i; @collect DataFrame
    end

    xmin=minimum(pgens[:Generation])
    xmax=maximum(pgens[:Generation])
    ymin=minimum(pgens[:Rlow])
    ymax=maximum(pgens[:Rhigh])

    plt = plot(pgens, x=:Generation, y=:Rmean, ymin=:Rlow, ymax=:Rhigh,
               color=:Steps, Geom.line, Geom.ribbon,
               Scale.color_discrete_manual(colors...),
               Guide.ylabel(ylabel), Guide.xlabel(xlabel),
               Coord.cartesian(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),
               Guide.title(problem))
    draw(PDF(outfile, 8inch, 6inch), plt);
    plt
end

function get_sizes(filename::String)
    names = [:id, :epochs, :gen, :size]
    res = readtable(filename, separator=',', header=false, names=names)

    stats = by(res, [:epochs, :gen], df->DataFrame(
        msize=mean(df[:size]), ssize=std_or_z(df[:size])))
    stats[:low_size] = stats[:msize] - 0.5 * stats[:ssize]
    stats[:high_size] = stats[:msize] + 0.5 * stats[:ssize]

    res, stats
end

function plot_sizes(stats::DataFrame, title::String="Boston Housing - GRN sizes",
                    logfile::String="boston_sizes.pdf";
                    xmin=minimum(stats[:gen]),
                    xmax=maximum(stats[:gen]),
                    ymin=minimum(stats[:low_size]),
                    ymax=maximum(stats[:high_size]),
                    ylabel="N",
                    xlabel="Generation",
                    key_position=:right)
    plt = plot(stats, x=:gen, y=:msize, ymin=:low_size, ymax=:high_size,
               color=:epochs, Geom.line, Geom.ribbon,
               Scale.color_discrete_manual(colors...),
               Guide.ylabel(ylabel), Guide.xlabel(xlabel),
               Coord.cartesian(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),
               Guide.title(title))
    draw(PDF(logfile, 8inch, 6inch), plt);
    plt
end

function get_error_res(filename::String, maxgen=0)
    # names = [:m, :time, :gen, :error, :accuracy, :mean, :std, :test]
    names = [:id, :epochs, :time, :gen, :fit, :test]
    error_res = readtable(filename, separator=',', header=false, names=names)
    # error_res[:fit] = log10.(error_res[:fit])
    error_res[:epochs] = string.(error_res[:epochs])
    bests = by(error_res, [:id, :epochs, :gen], df->DataFrame(
        error=minimum(df[:fit]), test_fit=minimum(df[:test])))

    error_stats = by(bests, [:epochs, :gen], df->DataFrame(
        merror=mean(df[:error]), serror=std_or_z(df[:error])))
    error_stats[:low_error] = error_stats[:merror] - 0.5 * error_stats[:serror]
    error_stats[:high_error] = error_stats[:merror] + 0.5 * error_stats[:serror]

    error_res, error_stats
end

function get_serror_res(filename::String, maxgen=0)
    # names = [:m, :time, :gen, :error, :accuracy, :mean, :std, :test]
    names = [:id, :epochs, :time, :gen, :epoch, :fit, :test]
    error_res = readtable(filename, separator=',', header=false, names=names)
    error_res[:epochs] = string.(error_res[:epochs])
    # error_res[:fit] = log10.(error_res[:fit])
    bests = by(error_res, [:id, :epochs, :gen], df->DataFrame(
        error=minimum(df[:fit]), test_fit=minimum(df[:test])))

    error_stats = by(bests, [:epochs, :gen], df->DataFrame(
        merror=mean(df[:error]), serror=std_or_z(df[:error])))
    error_stats[:low_error] = error_stats[:merror] - 0.5 * error_stats[:serror]
    error_stats[:high_error] = error_stats[:merror] + 0.5 * error_stats[:serror]

    error_res, error_stats
end

function get_learn_res(filename::String)
    all_res = readtable(filename, separator=',', header=false,
                        names = [:type, :id, :evo_epoch, :tgen, :time, :gen,
                                 :epoch, :loss, :test])

    not_grns = @from i in all_res begin
        @where i.type != "grn"
        @select i
        @collect DataFrame
    end;

    l_grns = @from i in all_res begin
        @where (i.type == "grn" && i.evo_epoch == "10" && i.tgen == 50)
        @select i
        @collect DataFrame
    end;

    res = vcat(not_grns, l_grns)

    stats = by(res, [:type, :evo_epoch, :tgen, :epoch], df->DataFrame(
        mloss=mean(df[:loss]), sloss=std(df[:loss]),
        mtest=mean(df[:test]), stest=std(df[:test])))
    stats[:low_loss] = stats[:mloss] - 0.5 * stats[:sloss]
    stats[:high_loss] = stats[:mloss] + 0.5 * stats[:sloss]
    stats[:low_test] = stats[:mtest] - 0.5 * stats[:stest]
    stats[:high_test] = stats[:mtest] + 0.5 * stats[:stest]

    res, stats
end

function get_final(filename::String)
    names = [:type, :id, :epoch, :tgen, :time, :gen, :mse, :test]
    final = readtable(filename, separator=',', header=false, names=names)
    final[:epoch] = string.(final[:epoch])
    final[:tgen][final[:tgen] .== 49] = 50
    final
end

function plot_final(final::DataFrame, title::String="Boston Housing - best GRNs",
                    filename::String="boston_final.pdf")
    grns = @from i in final begin
        @where i.type == "grn"
        @select i
        @collect DataFrame
    end;
    grns[:Generation] = string.(grns[:tgen])
    grns[:E] = "Train"
    grns[:yax] = grns[:mse]
    cgrns = deepcopy(grns)
    cgrns[:yax] = cgrns[:test]
    cgrns[:E] = "Test"

    full = vcat(grns, cgrns)

    plt = plot(full, xgroup=:epoch, ygroup=:E, x=:Generation, y=:yax, color=:epoch,
                     Scale.y_log10,
                     Guide.ylabel("Mean Squared Error"),
                     Guide.xlabel("Generation by number of epochs during evolution"),
                     Guide.title(title),
                     Guide.xlabel(nothing),
                     style(key_position=:none),
                     Geom.subplot_grid(Geom.boxplot),
                     Scale.color_discrete_manual(colors...));
    draw(PDF(filename, 18inch, 8inch), plt);
    plt
end

function plot_comparison(final::DataFrame, title::String="Boston Housing",
                    filename::String="boston_comparison.pdf", tgen=50)
    not_grns = @from i in final begin
        @where i.type != "grn"
        @select i
        @collect DataFrame
    end;

    l_grns = @from i in final begin
        @where (i.type == "grn" && i.epoch == "10" && i.tgen == tgen)
        @select i
        @collect DataFrame
    end;

    f = vcat(not_grns, l_grns);
    f[:const] = ""

    plt = plot(f, xgroup=:type, x=:const, y=:test,
               Scale.y_log10,
               Guide.xlabel("Model"),
               Guide.ylabel("Mean Squared Error"),
               # Geom.boxplot,
               # Scale.x_discrete(order=[4, 1, 2, 3]),
               Geom.subplot_grid(Geom.boxplot),
               Guide.title(title),
               style(key_position=:none, default_color=colorant"#000000"),
               Scale.color_discrete_manual(colors...));
    draw(PDF(filename, 8inch, 6inch), plt);
    plt
end

function plot_error(error_stats::DataFrame, title::String="Air Quality Evolution",
                    logfile::String="error.pdf";
                    xmin=minimum(error_stats[:gen]),
                    xmax=maximum(error_stats[:gen]),
                    ymin=log10(minimum(error_stats[:low_error])),
                    ymax=log10(maximum(error_stats[:high_error])),
                    ylabel="Mean Squared Error",
                    xlabel="Generation",
                    key_position=:right)
    zerr = @from i in error_stats begin
        @where ((i.epochs == "0") && (i.gen == 0))
        @select [i.merror, i.low_error, i.high_error]
        @collect
    end

    error_stats[:merror][error_stats[:gen] .== 0] = zerr[1][1].value
    error_stats[:low_error][error_stats[:gen] .== 0] = zerr[1][2].value
    error_stats[:high_error][error_stats[:gen] .== 0] = zerr[1][3].value

    plt = plot(error_stats, x=:gen, y=:merror, ymin=:low_error, ymax=:high_error,
               color=:epochs, Geom.line, Geom.ribbon,
               Scale.color_discrete_manual(colors...),
               Scale.y_log10,
               Guide.ylabel(ylabel), Guide.xlabel(xlabel),
               Coord.cartesian(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),
               Guide.title(title))
    draw(PDF(logfile, 8inch, 6inch), plt);
    plt
end

function plot_learn(learn_stats::DataFrame, title::String="Energy Training",
                    logfile::String="learn.pdf")

    plt = plot(learn_stats, x=:epoch, y=:mloss, ymin=:low_loss, ymax=:high_loss,
               color=:type, Geom.line, Geom.ribbon,
               Scale.color_discrete_manual(colors...),
               Guide.ylabel("Loss"),
               Guide.title(title))
    draw(PDF(logfile, 8inch, 6inch), plt);
    plt
end

function plot_learn_distributions(learn_res::DataFrame,
                                  title::String="Boston learning distributions",
                                  logfile::String="boston_learn_dist.pdf")

    downsampled = @from i in res begin
        @where ((mod(i.epoch, 10) == 1) && (i.epoch <= 51))
        @select i
        @collect DataFrame
    end;

    plt = plot(downsampled, xgroup=:epoch, x=:type, y=:loss, color=:type,
               Scale.y_log10,
               Guide.ylabel("Loss [log10]"),
               Geom.subplot_grid(Geom.boxplot), Guide.xlabel("Epoch"),
               Guide.title(title), Scale.color_discrete_manual(colors...));

    draw(PDF(logfile, 12inch, 6inch), plt);
    plt
end


text:
module ParallelProgressMeter

using ProgressMeter

export @parallelprogress

globalProgressMeters = Dict()
globalProgressValues = Dict()
globalPrintLock      = Dict()

"""
    @parallelprogress p for i = 1:50
        # parallel computation goes here
    end

displays progress in parallel `for` loop using the `Progress` object `p`.
"""
macro parallelprogress(args...)
    if length(args) < 1
        throw(ArgumentError("@parallelprogress requires at least one argument."))
    end

    progressargs = args[1:end-1]
    loop = args[end]
    metersym = gensym("meter")

    if !(isa(loop, Expr) && loop.head === :for)
        throw(ArgumentError("Final argument to @parallelprogress must be a for loop."))
    end

    id = randstring(20)

    setup = quote
        iterable = $(esc(loop.args[1].args[2]))
        $(esc(metersym)) = Progress(length(iterable), $([esc(arg) for arg in progressargs]...))

        ParallelProgressMeter.globalProgressMeters[$(esc(id))] = $(esc(metersym))
        ParallelProgressMeter.globalProgressValues[$(esc(id))] = 0
        ParallelProgressMeter.globalPrintLock[$(esc(id))]      = ReentrantLock()
    end

    loop.args[2] = ProgressMeter.showprogress_process_expr(loop.args[2], metersym)
    push!(loop.args[2].args, :(remotecall(ParallelProgressMeter.updateProgressMeter, 1, $(esc(id)), 1)))

    return quote
        $setup
        @sync @parallel $loop
    end
end

function updateProgressMeter(id, n)
    global globalProgressMeters
    global globalProgressValues
    global globalPrintLock

    lock(globalPrintLock[id])
    globalProgressValues[id] += n
    update!(globalProgressMeters[id], globalProgressValues[id])
    unlock(globalPrintLock[id])
end

end # module


text:
#!/usr/bin/env julia --project
# This experiment file establishes the initial conditions, boundary conditions,
# source terms and simulation parameters (domain size + resolution) for the

# Convective Boundary Layer LES case (Kitamura et al, 2016).

## ### Convective Boundary Layer LES
## [Nishizawa2018](@cite)
#
# To simulate the experiment, type in
#
# julia --project experiments/AtmosLES/convective_bl_les.jl

using ArgParse
using Distributions
using DocStringExtensions
using LinearAlgebra
using Printf
using Random
using StaticArrays
using UnPack
using Test

using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.ODESolvers
using ClimateMachine.Orientations
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.TurbulenceConvection
using ClimateMachine.VariableTemplates

using ClimateMachine.BalanceLaws
import ClimateMachine.BalanceLaws: source, prognostic_vars

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine.Atmos: altitude, recover_thermo_state

"""
    ConvectiveBL Geostrophic Forcing (Source)
"""
struct ConvectiveBLGeostrophic{FT} <: TendencyDef{Source}
    "Coriolis parameter [s⁻¹]"
    f_coriolis::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end
prognostic_vars(::ConvectiveBLGeostrophic) = (Momentum(),)

function source(::Momentum, s::ConvectiveBLGeostrophic, m, args)
    @unpack state, aux = args
    @unpack f_coriolis, u_geostrophic, u_slope, v_geostrophic = s

    z = altitude(m, aux)
    # Note z dependence of eastward geostrophic velocity
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(m, aux)
    fkvector = f_coriolis * ẑ
    # Accumulate sources
    return -fkvector × (state.ρu .- state.ρ * u_geo)
end

"""
  ConvectiveBL Sponge (Source)
"""
struct ConvectiveBLSponge{FT} <: TendencyDef{Source}
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Sponge exponent"
    γ::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end
prognostic_vars(::ConvectiveBLSponge) = (Momentum(),)

function source(::Momentum, s::ConvectiveBLSponge, m, args)
    @unpack state, aux = args

    @unpack z_max, z_sponge, α_max, γ = s
    @unpack u_geostrophic, u_slope, v_geostrophic = s

    z = altitude(m, aux)
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(m, aux)
    # Accumulate sources
    if z_sponge <= z
        r = (z - z_sponge) / (z_max - z_sponge)
        β_sponge = α_max * sinpi(r / 2)^s.γ
        return -β_sponge * (state.ρu .- state.ρ * u_geo)
    else
        FT = eltype(state)
        return SVector{3, FT}(0, 0, 0)
    end
end

"""
  Initial Condition for ConvectiveBoundaryLayer LES
"""
function init_problem!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord

    # Problem floating point precision
    param_set = parameter_set(bl)
    FT = eltype(state)
    R_gas::FT = R_d(param_set)
    c_p::FT = cp_d(param_set)
    c_v::FT = cv_d(param_set)
    p0::FT = MSLP(param_set)
    _grav::FT = grav(param_set)
    γ::FT = c_p / c_v
    # Initialise speeds [u = Eastward, v = Northward, w = Vertical]
    u::FT = 4
    v::FT = 0
    w::FT = 0
    # Assign piecewise quantities to θ_liq and q_tot
    θ_liq::FT = 0
    q_tot::FT = 0
    # functions for potential temperature

    θ_liq = FT(288) + FT(4 / 1000) * z
    θ = θ_liq
    π_exner = FT(1) - _grav / (c_p * θ) * z # exner pressure
    ρ = p0 / (R_gas * θ) * (π_exner)^(c_v / R_gas) # density
    # Establish thermodynamic state and moist phase partitioning
    TS = PhaseEquil_ρθq(param_set, ρ, θ_liq, q_tot)

    # Compute momentum contributions
    ρu = ρ * u
    ρv = ρ * v
    ρw = ρ * w

    # Compute energy contributions
    e_kin = FT(1 // 2) * (u^2 + v^2 + w^2)
    e_pot = _grav * z
    ρe_tot = ρ * total_energy(e_kin, e_pot, TS)

    # Assign initial conditions for prognostic state variables
    state.ρ = ρ
    state.ρu = SVector(ρu, ρv, ρw)
    state.energy.ρe = ρe_tot
    if !(moisture_model(bl) isa DryModel)
        state.moisture.ρq_tot = ρ * q_tot
    end

    if z <= FT(400) # Add random perturbations to bottom 400m of model
        state.energy.ρe += rand() * ρe_tot / 100
    end
    init_state_prognostic!(turbconv_model(bl), bl, state, aux, localgeo, t)
end

function surface_temperature_variation(bl, state, t)
    FT = eltype(state)
    ρ = state.ρ
    θ_liq_sfc = FT(291.15) + FT(20) * sinpi(FT(t / 12 / 3600))
    param_set = parameter_set(bl)
    if moisture_model(bl) isa DryModel
        TS = PhaseDry_ρθ(param_set, ρ, θ_liq_sfc)
    else
        q_tot = state.moisture.ρq_tot / ρ
        TS = PhaseEquil_ρθq(param_set, ρ, θ_liq_sfc, q_tot)
    end
    return air_temperature(TS)
end

function convective_bl_model(
    ::Type{FT},
    config_type,
    zmax,
    surface_flux;
    turbconv = NoTurbConv(),
    moisture_model = "dry",
) where {FT}

    ics = init_problem!     # Initial conditions

    C_smag = FT(0.23)     # Smagorinsky coefficient
    C_drag = FT(0.001)    # Momentum exchange coefficient
    z_sponge = FT(2560)     # Start of sponge layer

    α_max = FT(0.75)       # Strength of sponge layer (timescale)

    γ = 2                  # Strength of sponge layer (exponent)
    u_geostrophic = FT(4)        # Eastward relaxation speed
    u_slope = FT(0)              # Slope of altitude-dependent relaxation speed
    v_geostrophic = FT(0)        # Northward relaxation speed
    f_coriolis = FT(1.031e-4) # Coriolis parameter
    u_star = FT(0.3)
    q_sfc = FT(0)
    moisture_flux = FT(0)

    # Assemble source components
    source_default = (
        Gravity(),
        ConvectiveBLSponge{FT}(
            zmax,
            z_sponge,
            α_max,
            γ,
            u_geostrophic,
            u_slope,
            v_geostrophic,
        ),
        ConvectiveBLGeostrophic{FT}(
            f_coriolis,
            u_geostrophic,
            u_slope,
            v_geostrophic,
        ),
        turbconv_sources(turbconv)...,
    )

    if moisture_model == "dry"
        source = source_default
        moisture = DryModel()
    elseif moisture_model == "equilibrium"
        source = source_default
        moisture = EquilMoist(; maxiter = 5, tolerance = FT(0.1))
    elseif moisture_model == "nonequilibrium"
        source = (source_default..., CreateClouds())
        moisture = NonEquilMoist()
    else
        @warn @sprintf(
            """
%s: unrecognized moisture_model in source terms, using the defaults""",
            moisture_model,
        )
        source = source_default
    end
    # Set up problem initial and boundary conditions
    if surface_flux == "prescribed"
        energy_bc = PrescribedEnergyFlux((state, aux, t) -> LHF + SHF)
        moisture_bc = PrescribedMoistureFlux((state, aux, t) -> moisture_flux)
    elseif surface_flux == "bulk"
        energy_bc = BulkFormulaEnergy(
            (bl, state, aux, t, normPu_int) -> C_drag,
            (bl, state, aux, t) ->
                (surface_temperature_variation(bl, state, t), q_sfc),
        )
        moisture_bc = BulkFormulaMoisture(
            (state, aux, t, normPu_int) -> C_drag,
            (state, aux, t) -> q_sfc,
        )
    else
        @warn @sprintf(
            """
%s: unrecognized surface flux; using 'prescribed'""",
            surface_flux,
        )
    end

    # Define the physics
    physics = AtmosPhysics{FT}(
        param_set;
        turbulence = SmagorinskyLilly{FT}(C_smag),
        moisture = moisture,
        turbconv = turbconv,
    )

    moisture_bcs = moisture_model == "dry" ? () : (; moisture = moisture_bc)
    boundary_conditions = (
        AtmosBC(
            physics;
            momentum = Impenetrable(DragLaw(
                # normPu_int is the internal horizontal speed
                # P represents the projection onto the horizontal
                (state, aux, t, normPu_int) -> (u_star / normPu_int)^2,
            )),
            energy = energy_bc,
            moisture_bcs...,
            turbconv = turbconv_bcs(turbconv)[1],
        ),
        AtmosBC(physics; turbconv = turbconv_bcs(turbconv)[2]),
    )

    problem = AtmosProblem(
        init_state_prognostic = ics,
        boundaryconditions = boundary_conditions,
    )

    # Assemble model components
    model =
        AtmosModel{FT}(config_type, physics; problem = problem, source = source)

    return model
end

function config_diagnostics(driver_config)
    default_dgngrp = setup_atmos_default_diagnostics(
        AtmosLESConfigType(),
        "2500steps",
        driver_config.name,
    )
    core_dgngrp = setup_atmos_core_diagnostics(
        AtmosLESConfigType(),
        "2500steps",
        driver_config.name,
    )
    return ClimateMachine.DiagnosticsConfiguration([
        default_dgngrp,
        core_dgngrp,
    ])
end


text:
"""Module containing functions for creating and operating on a `RTree`, to quickly store and locate georeferenced polygons."""
module GeoTrees

using Revise
using ArchGDAL
using SpatialIndexing
using JLD2

const agd = ArchGDAL
const si = SpatialIndexing

#=
POTREBBE ESSERE UTILE AVERE DEI VETTORI DI APPOGGIO CHE CONTENGANO TUTTI I POSSIBILI VALORI DEI VARI ATTRIBUTI DEI POLIGONI,
IN QUESTO VERIFICARE SE UN POLIGONO CON UN CERTO VALORE DI UN CERTO ATTRIBUTO ESISTE DIVENTA TRIVIALE, INVECE DI RICHIEDERE
DI RICERCARE UN POLIGONO CON QUELL'ATTRIBUTO IN TUTTO L'ALBERO.
=#
mutable struct CCSTree
    terrain_types::Dict{Float64, String}
    tree::SpatialIndexing.RTree{Float64, 2}

    CCSTree(tree) = new(Dict{Float64, String}(), tree)
end

"""
    Base.show(io::IO, rtree::SpatialIndexing.RTree{T, N})

Write a text reppresentation of `SpatialIndexing.RTree{T, N}` object `rtree` to the output stream `io`.
"""
function Base.show( io::IO, rtree::SpatialIndexing.RTree{T, N} ) where {T, N}
    println(io, typeof(rtree), ":")
    println(io, "\tPolygons: ", rtree.nelems)
    println(io, "\tLevels: ", rtree.root.level)
    println(io, "\tNodes per level:")
    println(io, "\t", reverse(rtree.nnodes_perlevel) )
    println(io, "\tRoot MBR:", rtree.root.mbr)
end

"""
    Base.show(io::IO, ccstree::CCSTree)

Write a text reppresentation of `CCSTree` object `ccstree` to the output stream `io`.
"""
function Base.show( io::IO, ccstree::CCSTree )
    println(io, typeof(ccstree), ":")
    println(io, "\tTerrain typologies: ", length(ccstree.terrain_types))
    println(io, "\tPolygons: ", ccstree.tree.nelems)
    println(io, "\tLevels: ", ccstree.tree.root.level)
    println(io, "\tNodes per level:")
    println(io, "\t", reverse(ccstree.tree.nnodes_perlevel) )
    println(io, "\tRoot MBR:", ccstree.tree.root.mbr)
end

"""
    SpatialIndexing.check( ccstree::CCSTree )

Check wether a `CCSTree` is well formed.
"""
function check( ccstree::CCSTree )
    return si.check(ccstree.tree)
end

"""
    mbrtype( polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} )

Return the type of the mbr returned by function `mbr()` applied to `polygon`
"""
function mbrtype( polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} )
    return SpatialIndexing.HasMBR{ SpatialIndexing.Rect{Float64, 2} }
end

"""
    mbr( polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} )

Return the minimum bounding rectangle of `polygon` as a `SpatialIndexing.Rect{T, N}`
"""
function mbr( polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} )
    boundingbox = agd.envelope(polygon)
    return SpatialIndexing.Rect{Float64, 2}( (boundingbox.MinX, boundingbox.MinY), (boundingbox.MaxX, boundingbox.MaxY) )
end

"""
    printTree( node, n::Int64=0 )

Given a "SpatialIndexing.RTree" or "Node" either("Spatialindexing.Branch" or "Spatialindexing.Leaf") print the content of the tree with `node` as root.
"""
function printTree( node::Union{SpatialIndexing.RTree, SpatialIndexing.Node}, n::Int64=0 )
    if node isa SpatialIndexing.Leaf
        println("LEVEL: 1")
        println("MBR: $(node.mbr)")
        println("ELEMENTS:")
        for (i, el) in enumerate(node.children)
            println("$n.$i) $(el.id)  -  $(el.val[1])  -  $(el.mbr)")
        end
        println("\n")
    elseif node isa SpatialIndexing.Branch
        println("LEVEL: $(node.level)")
        println("MBR: $(node.mbr)")
        println("CHILDREN:")
        for (i, el) in enumerate(node.children)
            println("$n.$i")
            printTree(el, i)
        end
        println("\n")
    else
        printTree(node.root, 0)
    end
end

"""
    distance( element::Union{SpatialIndexing.SpatialElem, SpatialIndexing.Node}, region::SpatialIndexing.Region )

Compute the distance between `element` and `region`.
"""
distance( element::Union{SpatialIndexing.SpatialElem, SpatialIndexing.Node}, region::SpatialIndexing.Region ) = agd.distance(
    region isa SpatialIndexing.Point ? agd.createpoint(region.coord...) :
        let (xl, yl) = region.low, (xh, yh) = region.high
            agd.createpolygon([(xl, yh), (xh, yh), (xh, yl), (xl, yl), (xl, yh)])
        end,
    element isa SpatialIndexing.SpatialElem ? element.val[2] :
        let (xl, yl) = element.mbr.low, (xh, yh) = element.mbr.high
            agd.createpolygon([(xl, yh), (xh, yh), (xh, yl), (xl, yl), (xl, yh)])
        end
)

"""
    createTree( file_path::String, branch_capacity::Int64=7, leaf_capacity::Int64=7 )

Create a `SpatialIndexing.RTree{Float64, 2}` filled with the polygons contained in a shapefile at `file_path`.
"""
function createRTree( file_path::String, branch_capacity::Int64=7, leaf_capacity::Int64=7 )
 # Obtain all the "features" ("polygons" + additional information) of the file
    features = collect(agd.getlayer(agd.read(file_path), 0))
 # Create the RTree (It contains an empty vector as the "terrain_types" field and the empty tree that we just defined as the field "tree")
    tree = RTree{Float64, 2}(Float64, NamedTuple, branch_capacity=branch_capacity, leaf_capacity=leaf_capacity)
 # Insert all of the relevant data of the features in the tree as nodes
    for feature in features
     # Obtain the geometry (polygon) of the feature
        geom = agd.getgeom(feature, 0)
     # Find the MBR of the polygon
        feature_mbr = mbr(geom)
     # Check if the MBR has a valid value
        if si.isvalid(feature_mbr)
         # If so insert a new node in the tree
            si.insert!( tree, feature_mbr, agd.getfield(feature, :objectid), ( code=agd.getfield(feature, :codice_num), geometry=geom ) )
        else
            println("Omitted feature $(agd.getfield(feature, :objectid)), the feature mbr is not valid.")
        end
    end
 # Check that the tree is well-build
    return si.check(tree) ? tree : throw(DomainError("Error in the creation of the tree, check input parameters"))
end

"""
    createCCSTree( file_path::String, branch_capacity::Int64=7, leaf_capacity::Int64=7 )

Create a CCSTree based on a vectorial file at `file_path` containing a series of polygons decribing the terrain composition of an area.
The  created object will contain a field `terrain_types` containing a `Vector{Tuple{Int64, String}}` holding the codes and descriptions of
all the possible kinds of terrain found in the vectorial file.
The "tree" field of the object will instead contain a `SpatialIndexing.RTree{Float64, 2}` with all the polygons of the file.
"""
function createCCSTree( file_path::String, branch_capacity::Int64=7, leaf_capacity::Int64=7 )
 # Obtain all the "features" ("polygons" + additional information) of the file
    features = collect(agd.getlayer(agd.read(file_path), 0))
 # Create the CCSTree from an empty RTree
    ccstree = CCSTree( RTree{Float64, 2}(Float64, NamedTuple, branch_capacity=branch_capacity, leaf_capacity=leaf_capacity) )
  # Insert all of the relevant data of all the features in the tree as nodes
    for feature in features
     # Obtain the geometry (polygon) of the feature
        geom = agd.getgeom(feature, 0)
     # Find the MBR of the polygon
        feature_mbr = mbr(geom)
     # Check if the MBR has a valid value
        if si.isvalid(feature_mbr)
         # Get the code that identifies the kind of landcover of the polygon
            num_code = agd.getfield(feature, :codice_num)
         # Insert the new node corresponding to the feature in the "tree" field of the CCSTree
            si.insert!( ccstree.tree, feature_mbr, agd.getfield(feature, :objectid), ( code=num_code, geometry=geom ) )
         # Insert the landcover type in the dictionary in the field "terrain_types" of the CCSTree along with its description
            push!( ccstree.terrain_types, num_code => agd.getfield(feature, :legenda) )
        else
            println("Omitted feature $(agd.getfield(feature, :objectid)), the feature mbr is not valid.")
        end
    end
  # Check that the tree is well-build
    return check(ccstree) ? ccstree : throw(DomainError("Error in the creation of the tree, check input parameters"))
end

"""
    exists_terrain_type( tree::CCSTree, type_code::Int64 )

Check if the type of terrain corresponding to `type_code` exists in at least one polygon of `tree`.
"""
exists_terrain_type( tree::CCSTree, type_code::Int64 ) = type_code in keys(tree.terrain_types)

"""
    exists_terrain_type( tree::CCSTree, type_description::String )

Check if the type of terrain described by `type_description` exists in at least one polygon of `tree`.
Note: `type_description` must match exactly the description of the terrain type, if the complete description is not
known, consider using function `find_terrain_type_similar` instead.
"""
exists_terrain_type( tree::CCSTree, type_description::String ) = type_description in values(tree.terrain_types)

"""
    find_terrain_type_similar( tree::CCSTree, type::String )

Return a `Dict{Int64, String}` containing all the codes and the respective descriptions in which the latter contains `type` as a substring.
"""
find_terrain_type_similar( tree::CCSTree, type::String ) = filter( (k, v) -> contains(v, type), tree.terrain_types )

# --------------------------------------------------------------- POLYGON SEARCH ---------------------------------------------------------------------------------------
"""
    findPolygons( node::SpatialIndexing.Leaf{T,N}, polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} )::Vector{SpatialIndexing.SpatialElem} where {T, N}

Find and return a `Vector{SpatialIndexing.SpatialElem}` that contains all the `SpatialIndexing.SpatialElem`s children of `node` that intersect, contain or are containded in `polygon`.
"""
function findPolygons( node::SpatialIndexing.Leaf{T,N}, polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} )::Vector{SpatialIndexing.SpatialElem} where {T, N}
    # Return a Vector of all the children of the node that contain, are contained in, or intersect with `polygon`
    return [
        child
        for child in node.children
        if agd.contains(polygon, child.val.geometry) || agd.contains(child.val.geometry, polygon) || agd.intersects(polygon, child.val.geometry)
    ]
end

"""
    findPolygons( node::SpatialIndexing.Branch{T,N,V}, polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} ) where {T, N, V}

Find and return a `Vector{SpatialIndexing.SpatialElem}` that contains all the elements of the subtree rooted in `node` that intersect, contain or are contained in `polygon`.
"""
function findPolygons( node::SpatialIndexing.Branch{T,N,V}, polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} ) where {T, N, V}
    res = Vector{SpatialIndexing.SpatialElem}()
    for child in node.children
        poly_mbr = mbr(polygon)
        if si.intersects(child.mbr, poly_mbr) || si.contains(child.mbr, poly_mbr) || si.contains(poly_mbr, child.mbr)
            results = findPolygon(child, polygon)
            if !isempty(results)
                append!(res, results)
            end
        end
    end
    return res
end

"""
    findPolygons( tree::SpatialIndexing.RTree{T,N}, polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} ) where {T, N}

Return all the polygons of `tree` intersected by, contained in or containing `polygon`, or the polygon that conains it, or an empty `Vector{SpatialIndexing.SpatialElem}`, if there is no such polygon.
"""
function findPolygons( tree::SpatialIndexing.RTree{T,N}, polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} ) where {T, N}
    return findPolygon(tree.root, polygon)
end

"""
    findPolygons( tree::CCSTree, polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} ) where {T, N}

Return all the polygons of `tree` intersected by, contained in or containing `polygon`, or the polygon that conains it, or an empty `Vector{SpatialIndexing.SpatialElem}`, if there is no such polygon.
"""
function findPolygons( tree::CCSTree, polygon::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon} ) where {T, N}
    return findPolygon(tree.tree.root, polygon)
end

# ----------------------------------------------------------------- KNN SEARCH -----------------------------------------------------------------------------------------
"""
    findKNN( node::SpatialIndexing.Leaf{T,N}, point::SpatialIndexing.Point{T,N}, k::Int64 ) where {T, N}

Return the `k` `SpatialIndexing.SpatialElem`s children of `node` that are closest to `point` (at most k elements will be returned).
"""
function findKNN( node::SpatialIndexing.Leaf{T,N}, point::SpatialIndexing.Point{T,N}, k::Int64 ) where {T, N}
    agd_point = agd.createpoint(point.coord...)
    results = sort!( # 2) Sort the children in increasing order of distance from `point`
        map( # 1) Map every children (they will all be `SpatialElements`) of the node to the distance between its centroid and `point`
            element -> (
                element,
                agd.distance( # Compute distance between the centroid of the polygon and `point`
                    agd.centroid(element.val.geometry),
                    agd_point
                )
            ),
            node.children
        ),
        lt=(x, y) -> x[2] < y[2]
    )
    # 3) Return the, at most, `k` elements that are closer to `point`
    return results[1:min(k, length(results))]
end

"""
    findKNN( node::SpatialIndexing.Branch{T,N,V}, point::SpatialIndexing.Point{T,N}, k::Int64 ) where {T, N, V}

Return the `k` `SpatialIndexing.SpatialElem`s of the subtree rooted in `node` that are closest to `point` (at most k elements will be returned).
"""
function findKNN( node::SpatialIndexing.Branch{T,N,V}, point::SpatialIndexing.Point{T,N}, k::Int64 ) where {T, N, V}
    # 1) Order the children of the current node by increasing distance from point
    candidates = sort( node.children, lt=(x, y) -> distance(x, point) < distance(y, point) )
    results = sort!( # 5) Sort the resulting vector by increasing distance from `point`
        reduce( # 4) For each child of the current node a Vector of the results will be returned, join the vectors in a single one
            vcat,
            findKNN.( # 3) For each of the children (they could be `Leaf`s or other `Branch`s) find its `k` children closest to `point`
                candidates[1:min(k, length(candidates))], # 2) Take, at most, the first `k` children (if there are less than `k` children, take all of them)
                Ref(point),
                k
            )
        ),
        lt=(x, y) -> x[2] < y[2]
    )
    # 6) Take, at most, the first `k` results (if there are less than `k` children, take all of them)
    return results[1:min(k, length(results))]
end

"""
    findKNN( tree::SpatialIndexing.RTree{T,N}, point::SpatialIndexing.Point{T,N}, k::Int64 ) where {T, N}

Return the `k` `SpatialIndexing.SpatialElem`s contained in `tree` that are closest to `point` (at most k elements will be returned).
"""
function findKNN( tree::SpatialIndexing.RTree{T,N}, point::SpatialIndexing.Point{T,N}, k::Int64 ) where {T, N}
    return findKNN(tree.root, point, k)
end

"""
    findKNN( tree::CCSTree, point::SpatialIndexing.Point{T,N}, k::Int64 ) where {T, N}

Return the `k` `SpatialIndexing.SpatialElem`s contained in `tree` that are closest to `point` (at most k elements will be returned).
"""
function findKNN( tree::CCSTree, point::SpatialIndexing.Point{T,N}, k::Int64 ) where {T, N}
    return findKNN(tree.tree.root, point, k)
end

# ------------------------------------------------------------- DISTANCE BASED SEARCH ----------------------------------------------------------------------------------
"""
    findPolygons( node::SpatialIndexing.Leaf{T,N}, point::SpatialIndexing.Point{T,N} distance::Float64 ) where {T, N}

Return all the `SpatialIndexing.SpatialElem`s children of `node` that are within `distance` from point.
"""
function findPolygons( node::SpatialIndexing.Leaf{T,N}, point::SpatialIndexing.Point{T,N}, distance::Float64 ) where {T, N}
    agd_point = agd.createpoint(point.coord...)
    results = sort!( # 2) Sort the children in increasing order of distance from `point`
        map( # 1) Map every children (they will all be `SpatialElements`) of the node to the distance between its centroid and `point`
            element -> (
                element,
                agd.distance( # Compute distance between the centroid of the polygon and `point`
                    agd.centroid(element.val.geometry),
                    agd_point
                )
            ),
            node.children
        ),
        lt=(x, y) -> x[2] < y[2]
    )
    # 3) Return all the children within `distance` from `point`
    return results[1:something( findfirst(res -> res[2] > distance, result), length(result) )]
end

"""
    findPolygons( node::SpatialIndexing.Branch{T,N,V}, point::SpatialIndexing.Point{T,N} distance::Float64 ) where {T, N, V}

Return all the `SpatialIndexing.SpatialElem`s of the subtree rooted in `node` that are within `distance` from `point`.
"""
function findPolygons( node::SpatialIndexing.Branch{T,N,V}, point::SpatialIndexing.Point{T,N}, distance::Float64 ) where {T, N, V}
    results = sort!( # 5) Sort the results by distance from `point`
        reduce( # 4) For each child of the current node a Vector of the results will be returned, join the vectors in a single one
            vcat,
            findPolygons.( # 3) For each of the children (they could be `Leaf`s or other `Branch`s) find all its children whithin `distance` from `point`
                sort( # 1) Sort the children by increasing distance from `point`
                    node.children,
                    lt=(x, y) -> distance(x, point) < distance(y, point)
                )[1:something( findfirst(res -> res[2] > distance, result), length(result) )], # 2) All children up to the first farther from `point` than `distance`
                Ref(point),
                distance
            )
        ),
        lt=(x, y) -> x[2] < y[2]
    )
    # 6) Return all the results up to the first farther than `distance` from `point`
    return results[1:something( findfirst(res -> res[2] > distance, result), length(result) )]
end

"""
    findPolygons( tree::SpatialIndexing.RTree{T,N}, point::SpatialIndexing.Point{T,N}, distance::Float64 ) where {T, N}

Return all the `SpatialIndexing.SpatialElem`s contained in `tree` that are within `distance` from `point`.
"""
function findPolygons( tree::SpatialIndexing.RTree{T,N}, point::SpatialIndexing.Point{T,N}, distance::Float64 ) where {T, N}
    return findPolygons(tree.root, point, distance)
end

"""
    findPolygons( tree::CCSTree, point::SpatialIndexing.Point{T,N}, distance::Float64 ) where {T, N}

Return all the `SpatialIndexing.SpatialElem`s contained in `tree` that are within `distance` from `point`.
"""
function findPolygons( tree::CCSTree, point::SpatialIndexing.Point{T,N}, distance::Float64 ) where {T, N}
    return findPolygons(tree.tree.root, point, distance)
end

# ----------------------------------------------------------- LOADING AND SAVING -----------------------------------------------------------------------------
"""
    saveTree( tree::Union{SpatialIndexing.RTree{T, N}, CCSTree}, output_file::AbstractString=".\\tree.jld2" ) where {T, N}

Save `tree` as a file in the directory and with the name specified by `output_file` in Julia Data format 2 (JLD2).
"""
function saveTree( tree::Union{SpatialIndexing.RTree{T, N}, CCSTree}, output_file::AbstractString=".\\tree.jld2" ) where {T, N}
    if !occursin(output_file, ".jld2")
        output_file *= ".jld2"
    end
    save_object(tree, output_file)
end

"""
    loadTree( tree_file::AbstractString )

Load a `SpatialIndexing.RTree{T, N}` or a `CCSTree` saved as `tree_file`.
"""
function loadTree( tree_file::AbstractString )
    return load_object(tree_file)
end

#=
# Serve un modo per simulare la possibilità di aggiungere campi in modo da poter creare strutture che si adattino al numero di
# attributi rilevanti di cui si vuole ottenere la lista degli unici.
mutable struct Fields
    fields::Dict{Symbol, Any}

    Fields(fields) = new(fields)
end

mutable struct Cesto
    frutti::Fields

    function Cesto( frutti..., numeri... )
        if length(frutti) != length(numeri)
            throw(DomainError("`frutti` e `numeri` per ogni frutto deve esserci un numero"))
        end
        for (f, n) in zip(frutti, numeri)
            dict = Dict( Symbol(f) => n )
        end
    end
end
=#

using ArchGDAL

const agd = ArchGDAL

ccs_file = "C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\ccs\\ccs.shp"
ccs = agd.read(ccs_file)
features =  collect(agd.getlayer(ccs, 0))
features[1]

ccstree = createCCSTree("C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\ccs\\ccs.shp")

rtree = createRTree("C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\ccs\\ccs.shp")

#= ---------------------------- TREE TESTING [V] --------------------------------

#   ccs_shp = agd.read("C:\\Users\\Lenovo\\Documents\\GitHub\\Tirocinio\\resources\\Analysis data\\ccs\\ccs.shp")
ccs_shp = agd.read("D:\\Z_Tirocinio_Dati\\ccs WGS84\\ccs.shp")
features =  collect(agd.getlayer(ccs_shp, 0))

tree = nothing
GC.gc()
tree = RTree{Float64, 2}(Float64, NamedTuple, branch_capacity=7, leaf_capacity=7)
for feature in features
    geom = agd.getgeom(feature, 0)
    feature_mbr = mbr(geom)
    if si.isvalid(feature_mbr)
        si.insert!( tree, feature_mbr, agd.getfield(feature, :objectid), ( code=agd.getfield(feature, :codice_num), geometry=geom ) )
    end
end
si.check(tree)
tree.nnodes_perlevel

=#
#= ---------------------------- QUERY TESTING [V] --------------------------------

# Multilinea di contorno al poligono
geom = agd.getgeom( agd.getgeom(features[311113]), 0 )
# NUmber of vertexes of the poligon
num_points = agd.ngeom(geom)
# Coordinates of points inside the polygon:
#   Average of each coordinate of the vertexes
xt1 = sum( j -> agd.getx(geom, j), 0:num_points-1 ) / Float64(num_points)
yt1 = sum( j -> agd.gety(geom, j), 0:num_points-1 ) / Float64(num_points)
#   Average of 3 vertexes
xt2 = sum( j -> agd.getx(geom, j), 0:2 ) / 3.0
yt2 = sum( j -> agd.gety(geom, j), 0:2 ) / 3.0
#   In the middle of an edge (mean of 2 vertexes)
xt3 = sum( j -> agd.getx(geom, j), 0:1 ) / 2.0
yt3 = sum( j -> agd.gety(geom, j), 0:1 ) / 2.0
# In the first half of the polygon
xt4 = sum( j -> agd.getx(geom, j), 0:(num_points÷2)-1 ) / (num_points÷2)
yt4 = sum( j -> agd.gety(geom, j), 0:(num_points÷2)-1 ) / (num_points÷2)

#= Results for each ponint
    res = findPolygon(tree, si.Point((xt1, yt1)))
    res = findPolygon(tree, si.Point((xt2, yt2)))
    res = findPolygon(tree, si.Point((xt3, yt3)))
    res = findPolygon(tree, si.Point((xt4, yt4)))
=#

res1 = knn( tree, si.Point((xt1, yt1)), 3 )
res2 = knn( tree, si.Point((xt2, yt2)), 3 )
res3 = knn( tree, si.Point((xt3, yt3)), 3 )
res4 = knn( tree, si.Point((xt4, yt4)), 3 )

sat_shp = agd.read(split( @__DIR__ , "\\Porting\\")[1] * "\\Mappe\\sat WGS84\\sette_sorelle.shp")
sat = agd.getgeom(collect(agd.getlayer(sat_shp, 0))[1])
res = findPolygon(tree, sat)
# Get the result in decreasing order based on the surface of the polygon
sort!( res, lt=(x, y) -> agd.geomarea(x.val.geometry) < agd.geomarea(y.val.geometry), rev=true )

# Id di tutti i poligoni intersecati dal poligono "sette_sorelle" (ottenuti con QGIS)
objectids = [
    192271, 192307, 197591, 214336, 214356, 214422,
    214518, 214606, 214777, 214840, 214967, 214999,
    215264, 215314, 215391, 215555, 215918, 215944,
    221275, 221300, 221335, 221415, 221641, 222633,
    223271, 231330, 232459, 232478, 232492, 232541,
    240004, 240325, 240533, 240803, 240903, 241012,
    241266, 241369, 241738, 249519, 249534, 249542,
    249588, 249619, 249622, 252658, 252687, 252737,
    252754, 252931, 258981, 261396, 261404, 261417,
    262501, 263548, 263567, 263578, 263584, 263602,
    263625, 265696, 265697, 265699, 265700, 265701,
    265702, 265703, 265704, 265705, 265706, 265708,
    265709, 265711, 265712, 265714, 265715, 267016,
    267021, 267029, 267043, 267064, 395922, 396242,
    396281
]

# Check that all the returned polygons are correct
all(r -> r.id in objectids, res)
resids = [ Int64(r.id) for r in res ]
# Check whether there are any missing polygon
all(id -> id in resids, objectids)

#=  PER SALVARE GLI ID TROVATI
    path = "C:\\Users\\DAVIDE-FAVARO\\Desktop\\Dati\\res_ids.txt"
    path = "D:\\Documents and Settings\\DAVIDE-FAVARO\\My Documents\\GitHub\\Tirocinio\\Mappe\\polygon_ids.txt"
    open( path, "w" ) do io
        for id in objectids
            write(io, "$id\n")
        end
    end
=#

=#
#= ---------------------------- SAVE TESTING [V] --------------------------------

path = "D:\\Z_Tirocinio_Dati\\tree_test.jld2"
save_object(path, tree)
loaded_tree = load_object(path);
si.check(loaded_tree)
l_res = findPolygon(loaded_tree, si.Point((xt1, yt1)))
l_res = findPolygon(loaded_tree, si.Point((xt2, yt2)))
l_res = findPolygon(loaded_tree, si.Point((xt3, yt3)))

=#
# ---------------------------- PLOT TESTING [V] --------------------------------
# PRESO DA "plot_utils.jl" CONTENUTO IN "https://github.com/alyst/SpatialIndexing.jl/tree/master/examples"

using PlotlyJS
using Printf: @sprintf

const SI = SpatialIndexing

# webcolor palette for R-tree levels
const LevelsPalette = Dict(
    1 => "#228B22", # Forest Green
    2 => "#DAA520", # Goldenrod
    3 => "#FF6347", # Tomato
    4 => "#B22222", # Firebrick
    5 => "#4B0082", # Indigo
    6 => "#800080", # Purple
    7 => "#4169E1", # Royal blue
    8 => "#008080", # Teal
    9 => "#00FFFF", # Cyan
)

# create plotly trace for a single R-tree `node` (rectangle edges)
function node_trace(node::SI.Node{<:Any,2}, ix::Union{Integer, Nothing};
                    showlegend::Bool)
    nbr = SI.mbr(node)
    res = scatter(x = [nbr.low[1], nbr.high[1], nbr.high[1], nbr.low[1], nbr.low[1]],
                  y = [nbr.low[2], nbr.low[2], nbr.high[2], nbr.high[2], nbr.low[2]],
                  mode=:lines, line_color=get(LevelsPalette, SI.level(node), "#708090"),
                  line_width=SI.level(node),
                  hoverinfo=:text, hoveron=:points,
                  text="lev=$(SI.level(node)) ix=$(ix !== nothing ? ix : "none")<br>nchildren=$(length(node)) nelems=$(SI.nelements(node)) area=$(@sprintf("%.2f", SI.area(SI.mbr(node))))",
                  name="level $(SI.level(node))",
                  legendgroup="level $(SI.level(node))", showlegend=showlegend)
    return res
end

# create plotly trace for a single R-tree `node` (cube edges)
function node_trace(node::SI.Node, ix::Union{Integer, Nothing};
                    showlegend::Bool)
    nbr = SI.mbr(node)
    res = scatter3d(x = [nbr.low[1],  nbr.high[1], nbr.high[1], nbr.low[1], nbr.low[1],
                         nbr.low[1],  nbr.high[1], nbr.high[1], nbr.low[1], nbr.low[1],
                         NaN, nbr.low[1], nbr.low[1],
                         NaN, nbr.high[1], nbr.high[1],
                         NaN, nbr.high[1], nbr.high[1]],
                    y = [nbr.low[2],  nbr.low[2],  nbr.high[2], nbr.high[2], nbr.low[2],
                         nbr.low[2],  nbr.low[2],  nbr.high[2], nbr.high[2], nbr.low[2],
                         NaN, nbr.high[2], nbr.high[2],
                         NaN, nbr.high[2], nbr.high[2],
                         NaN, nbr.low[2], nbr.low[2]],
                    z = [nbr.low[3],  nbr.low[3],  nbr.low[3],  nbr.low[3], nbr.low[3],
                         nbr.high[3], nbr.high[3], nbr.high[3], nbr.high[3], nbr.high[3],
                         NaN, nbr.low[3], nbr.high[3],
                         NaN, nbr.low[3], nbr.high[3],
                         NaN, nbr.low[3], nbr.high[3]],
                  mode=:lines, line_color=get(LevelsPalette, SI.level(node), "#708090"),
                  line_width=SI.level(node),
                  hoverinfo=:text, hoveron=:points,
                  text="lev=$(SI.level(node)) ix=$(ix !== nothing ? ix : "none")<br>nchildren=$(length(node)) nelems=$(SI.nelements(node)) area=$(@sprintf("%.2f", SI.area(SI.mbr(node))))",
                  name="level $(SI.level(node))",
                  legendgroup="level $(SI.level(node))", showlegend=showlegend)
    return res
end

# create plotly traces for the nodes in a subtree with the `node` root
# and append them to `node_traces`
function append_subtree_traces!(node_traces::Vector{PlotlyBase.AbstractTrace},
                                node::SI.Node, ix::Union{Integer, Nothing},
                                levels::Set{Int})
    push!(node_traces, node_trace(node, ix, showlegend = SI.level(node) ∉ levels))
    push!(levels, SI.level(node)) # show each level once in legend
    if node isa SI.Branch
        for (i, child) in enumerate(SI.children(node))
            append_subtree_traces!(node_traces, child, i, levels)
        end
    end
    return nothing
end

data_trace(tree::RTree{<:Any, 2}) =
    scatter(mode=:markers, name=:data, marker_color = "#333333", marker_size = 2,
            x=[SI.center(SI.mbr(x)).coord[1] for x in tree],
            y=[SI.center(SI.mbr(x)).coord[2] for x in tree],
            text=["id=$(SI.id(x))" for x in tree], hoverinfo=:text)

data_trace(tree::RTree) =
    scatter3d(mode=:markers, name=:data, marker_color = "#333333", marker_size = 2,
              x=[SI.center(SI.mbr(x)).coord[1] for x in tree],
              y=[SI.center(SI.mbr(x)).coord[2] for x in tree],
              z=[SI.center(SI.mbr(x)).coord[3] for x in tree],
              text=["id=$(SI.id(x))" for x in tree], hoverinfo=:text)

# create Plotly plot of the given tree
function PlotlyJS.plot(tree::RTree)
    ndims(tree) == 1 && throw(ArgumentError("1-D R-trees not supported"))
    ndims(tree) > 3 && @warn("Only 1st-3rd dimensions would be shown for $(ndims(tree))-D trees")
    node_traces = Vector{PlotlyBase.AbstractTrace}()
    append_subtree_traces!(node_traces, tree.root, nothing, Set{Int}())
    PlotlyJS.plot([node_traces; [data_trace(tree)]], Layout(hovermode=:closest))
end

# MOLTO PESANTE A LIVELLO DI PRESTAZIONI
tree_plot = plot(tree);
open( "C:\\Users\\DAVIDE-FAVARO\\Desktop\\Dati\\tree_plot.png", "w" ) do io
    PlotlyJS.savefig(io, tree_plot, width=900, height=1000)
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
module FSIDriversTest

using GridapFSI.Drivers
using Test

problem = FSIProblem{:not_implemented_problem}()
@test_throws ErrorException execute(problem)
kwargs = Dict(:a=>1.0)
@test Drivers._get_kwarg(:a,kwargs) == 1.0
@test Drivers._get_kwarg(:b,kwargs,2.0) == 2.0
@test_throws UndefVarError Drivers._get_kwarg(:b,kwargs)

end


text:
""" A singleton type that indicates `ScratchSpace` is deallocated. """
struct Cleared end

"""
    ScratchSpace(f::F, value::T)
    ScratchSpace{F,T}(f::F)

An object for appropriately caching the `value` returned from `f`. The
`value` is thrown away when `serialize`d.

Construct this by `x = ScratchSpace(f, f())`. At the use site, update the
variable with `x = allocated(x)` and then fetch the value by `x.value`.
"""
struct ScratchSpace{F,T} <: AbstractScratchSpace
    f::F
    value::Union{T,Cleared}
    ScratchSpace(f::F, value::T) where {F,T} = new{F,T}(f, value)
    # ScratchSpace{F,T}(f::F, value::T) where {F,T} = new{F,T}(f, value)
    ScratchSpace{F,T}(f::F) where {F,T} = new{F,T}(f, Cleared())  # deallocate
end
# Note: Subtyping `AbstractScratchSpace` so that FoldsCUDA can remove
# `ScratchSpace` (that can contain mutables) before sharing the accumulator
# between threads without importing FLoops.jl.
# Note: Using `::Union{T,Cleared}` so that serialize-deserialize does not
# change the type (after throwing away the `.value`). It seems that the
# compiler can optimize away the type-instability.

@inline function allocate(x::ScratchSpace)
    if x.value isa Cleared
        return ScratchSpace(x.f, x.f())
    else
        return x
    end
end

# Idea: use this instead of `x.value` to help(?) the compiler
#=
@inline function getvalue(x::ScratchSpace)
    value = x.value
    if value isa Cleared
        return x.f()
    else
        return value
    end
end
=#

# Throw away temporary object `.value` when crossing the wire (i.e., called via
# Distributed).
Serialization.serialize(ser::S, x::ScratchSpace{F,T}) where {S<:AbstractSerializer,F,T} =
    invoke(Serialization.serialize, Tuple{S,Any}, ser, ScratchSpace{F,T}(x.f))


text:
const LIGHT_SPEED = 299792458 # m/s^2

"""
    lerp(t, a, b)

Linear interpolation on the interval `[a, b]` as `t ∈ [0, 1].

"""
lerp(t, a, b) = (one(t)-t)*a + t*b

"""
    quadratic(A, B, C)

Solve quadratic equation `A x^2 + B x + C = 0` for real coefficients.

If a real solution exists, returns `(true, t0, t1)` where `t0` is the smaller
and `t1` the larger solution. If no real solution exists, returns `(false, Inf, Inf)`.

"""
function quadratic(A::Real, B::Real, C::Real)
    dd = B^2 - 4*A*C
    if dd < 0.0
        return false, Inf, Inf
    end
    q = B<0 ? -0.5 * (B - sqrt(dd)) : -0.5 * (B + sqrt(dd))
    t0 = q/A
    t1 = C/q
    return true, min(t0,t1), max(t0,t1)
end

"""
    round_pow2(n::Integer)

Round an integer up to nearest power of two.
"""
round_pow2(n::Integer) = 2^convert(typeof(n), ceil(log(2, n)))

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# Practical Machine learning
# Decision Trees based learning - Random Forests example
# Chapter 5

import .Trees: TabularData, Tree, Node, Leaf, Criterion, Gini, CrossEntropy, MSE, getroot, getleft, getright, isnode, n_samples, impurity

type RandomForest{T}
    # parameters
    n_estimators::Int
    max_features::Any
    max_depth::Int
    min_samples_split::Int
    criterion::Criterion

    # learner
    learner::Union{T,Void}

    function RandomForest(n_estimators, max_features, max_depth, min_samples_split, criterion)
        if n_estimators < 1
            error("n_estimators is too small (got: $n_estimators)")
        end

        if isa(max_features, Integer) && max_features < 1
            error("max_features is too small (got: $max_features)")
        elseif isa(max_features, AbstractFloat) && !(0. < max_features <= 1.)
            error("max_features should be in (0, 1] (got: $max_features)")
        elseif isa(max_features, Symbol) && !(is(max_features, :sqrt) || is(max_features, :third))
            error("max_features should be :sqrt or :third (got: $max_features)")
        end

        if is(max_depth, nothing)
            max_depth = typemax(Int)
        elseif isa(max_depth, Integer) && max_depth <= 1
            error("max_depth is too small (got: $max_depth)")
        end

        if min_samples_split <= 1
            error("min_sample_split is too small (got: $min_samples_split)")
        end

        if is(criterion, :gini)
            criterion = Gini
        elseif is(criterion, :entropy)
            criterion = CrossEntropy
        elseif is(criterion, :mse)
            criterion = MSE
        else
            error("criterion is invalid (got: $criterion)")
        end

        new(n_estimators, max_features, max_depth, min_samples_split, criterion, nothing)
    end
end

function resolve_max_features(max_features::Any, n_features::Int)
    if is(max_features, :sqrt)
        floor(Int, sqrt(n_features))
    elseif is(max_features, :third)
        div(n_features, 3)
    elseif isa(max_features, Integer)
        max(max_features, n_features)
    elseif isa(max_features, AbstractFloat)
        floor(Int, n_features * max_features)
    elseif is(max_features, nothing)
        n_features
    else
        error("max_features is invalid: $max_features")
    end
end

function feature_importances(rf::RandomForest)
    if is(rf.learner, nothing)
        error("not yet trained")
    end
    rf.learner.improvements
end

function set_improvements!(learner)
    improvements = learner.improvements

    for tree in learner.trees
        root = getroot(tree)
        add_improvements!(tree, root, improvements)
    end
    normalize!(improvements)
end

function add_improvements!(tree::Tree, node::Node, improvements::Vector{Float64})
    left = getleft(tree, node)
    right = getright(tree, node)
    n_left_samples = n_samples(left)
    n_right_samples = n_samples(right)
    averaged_impurity = (impurity(left) * n_left_samples + impurity(right) * n_right_samples) / (n_left_samples + n_right_samples)
    improvement = impurity(node) - averaged_impurity
    improvements[node.feature] += improvement

    add_improvements!(tree, left, improvements)
    add_improvements!(tree, right, improvements)
    return
end

function add_improvements!(::Tree, ::Leaf, ::Vector{Float64})
    # do nothing!
    return
end


text:
# Two main data structures used here
#   (1) RatingsList  -- stores an array of players, and a dictionary of their ratings
#   (2) RatingsTable -- for storing an array of ratings (with players as column headings)

#####################################
#  RatingsList
#####################################

"""
    RatingsList

 Self-contained set of players and their ratings

## Parameters
* `players::Array{String,1}`: ordered list of player names
* `ratings::Dict{String, Float64}`: the ratings of the players (where these have been set)
```
"""
struct RatingsList
    players::Array{String,1} # names of players
    ratings::Dict{String, Float64}
#    K::Dict{String, Float64}
#    β::Dict{String, Float64}
end
copy(r::RatingsList) = RatingsList( r.players, r.ratings )
length(r::RatingsList) = length(r.players)
==(r1::RatingsList, r2::RatingsList) = (r1.players == r2.players) && (r1.ratings == r2.ratings)

# get inverse mapping from an array of players to their indexes
function player_indexes( players::Array{String,1} )
    d = Dict{String, Int64}()
    for i in 1:length(players)
        d[players[i]] = i
    end
    return d
end
player_indexes( R::RatingsList ) = player_indexes( R.players )

# create a blank ratings list from a list of players
function RatingsList( player_list::Dict{String,Int64} )
    tmp = sort( collect( keys(player_list) ) )
    return RatingsList( tmp, Dict{String, Int64}() )
end

#####################################
#  RatingsTable
#####################################
"""
    RatingsTable

 Table of ratings (with players as column headings)

## Parameters
* `players::Array{String,1}`: ordered list of player names
* `ratings::DataFrame`: containing a series of ratings for each player as a column
```
"""
struct RatingsTable
    players::Array{String,1}
    ratings::DataFrame
end
copy(R::RatingsTable) = RatingsTable( R.players, R.ratings )
size(R::RatingsTable) = size(R.ratings)
size(R::RatingsTable, d) = size(R.ratings, d)
lastindex(R::RatingsTable) = size(R.ratings, 1)

# create empty RatingsTable with players as names of columns
function RatingsTable(players::Array{String,1})
    R = DataFrame( Matrix{Union{Missing,Float64}}(undef, 0, length(players)) )
    names!(R, Symbol.(players))
    return RatingsTable(players, R)
end

# create empty RatingsTable with players as names of columns, and n empty rows
function RatingsTable(players::Array{String,1}, n::Int)
    R = DataFrame( Matrix{Union{Missing,Float64}}(undef, n, length(players)) )
    names!(R, Symbol.(players))
    return RatingsTable(players, R)
end

# create a RatingsTable from a RatingsList (with same players, and one row giving the input ratings)
function RatingsTable(r::RatingsList)
    A = DataFrame( Matrix{Union{Missing,Float64}}(undef, 1, length(r)) )
    names!(A, Symbol.(r.players))
    for k in keys(r.ratings)
        A[ Symbol(k) ] = r.ratings[k]
    end
    return RatingsTable(r.players, A)
end

# append one RatingsTable onto another (they must have matching players lists)
function append!(R1::RatingsTable, R2::RatingsTable)
    append!(R1.ratings, R2.ratings)
    return R1
end

# add a row to a RatingsTable from a RatingsList (they must have matching players lists)
function push!(R::RatingsTable, r::RatingsList)
    if R.players != r.players
        error("must have matching player lists")
    end
    x = RatingsTable(r)
    append!(R, x)
    return R
end

# insert a row into a RatingsTable from a ratings list -- not sure if this is needed?
# function insert!(R::RatingsTable, i::Int, r::RatingsList) # insert a row into a existing RatingsTable
#    not implemented as I haven't needed it yet
# end

# address rows of a RatingsTable
getindex(R::RatingsTable, i::Int) = R.ratings[i,:]
getindex(R::RatingsTable, range::UnitRange{Int}) = R.ratings[range,:]
function setindex!(R::RatingsTable, r::RatingsList, i::Int)
    R.ratings[i,:] = RatingsTable(r).ratings
end

# compare RatingsTables
==(R1::RatingsTable, R2::RatingsTable) = (R1.players == R2.players) && (R1.ratings == R2.ratings)


text:
#############################################################################
#############################################################################
#
# This file implements polynomial GCD
#
#############################################################################
#############################################################################

"""
The extended euclid algorithm for polynomials modulo prime.
"""

function extended_euclid_alg(a::PolynomialModP, b::PolynomialModP)
    prime = a.mod
    old_r, r = a, b
    old_s, s = one(PolynomialModP,prime), zero(PolynomialModP,prime)
    old_t, t = zero(PolynomialModP,prime), one(PolynomialModP,prime)

    while !iszero(r)

        q = divide(old_r, r) |> first

        old_r, r = r, old_r - q*r

        old_s, s = s, old_s - q*s

        old_t, t = t, old_t - q*t

    end
    g, s, t = old_r, old_s, old_t
    @assert s*a + t*b - g == 0
    return g, s, t
end

"""
The GCD of two polynomials modulo prime.
"""

gcd(a::PolynomialModP, b::PolynomialModP) = extended_euclid_alg(a,b) |> first


text:
module MolenaarRemake2
using JuMP
using KNITRO

sub_out_con = 5.0

function molenaar_remake2(sub_out_con)

    molenaar = Model(
        optimizer_with_attributes(
            KNITRO.Optimizer,
            "ms_enable" => 1,
            "opttol" => 1E-16,
            "opttolabs" => 1e-16,
            "honorbnds" => 1,
            "ms_maxsolves" => 10,
        ),
    )
    JuMP.set_silent(molenaar)

    # Constraints (used the same for now)
    UB = 100
    LB = 1e-8

    @variables molenaar begin
        LB <= trS <= UB # substrate transporter
        LB <= rib <= UB # ribosome
        LB <= catef <= UB # catalytic enzyme
        LB <= metef <= UB # metabolic enzyme
        LB <= prc <= UB # precursor enzyme (activating enzyme)
        LB <= lpb <= UB # lipid synthesis enzyme

        LB <= Si <= UB # intracellular substrate
        LB <= M <= UB # intermediate
        LB <= P <= UB # precursor
        LB <= lip <= UB # Lipid

        LB <= atp <= UB
        LB <= adp <= UB

        LB <= S <= UB # substrate concentration outside

        LB <= mu <= UB

        0 <= alpha[1:6] <= 1 # fraction of ribosome engaged in synthesis

        LB <= beta <= UB

        LB <= v_trS <= UB # transporter rate
        LB <= v_lpb <= UB # lipid biosynthesis rate
        LB <= v_rib <= UB # ribosomale rate
        LB <= v_catef <= UB # catabolic enzyme rate
        LB <= v_metef <= UB # metabolic enzyme rate
        LB <= v_prc <= UB
    end

    @NLparameters molenaar begin
        kcat_trS == 7
        kcat_prc == 5
        kcat_rib == 3
        kcat_catef == 10
        kcat_metef == 5
        kcat_lpb == 5

        Km_trS == 1
        Km_prc == 1
        Km_rib == 1
        Km_catef == 1
        Km_metef == 1
        Km_lpb == 1

        Km_atp == 1   # atp affinity of the activating enzyme prc
        Km_adp == 0.5 # adp affinity of catef and metef
        gamm == 0.6 # efficiency of catef

        sA_trS == 1 # surface area of the transporter
        sA_lip == 1 # surface area of the lipids

        S == sub_out_con
    end

    @NLconstraints molenaar begin
        (v_trS - v_metef - v_catef) - mu * Si == 0
        (v_prc - v_rib - v_lpb) - mu * P == 0
        (v_catef + v_metef - v_prc) - mu * M == 0
        (v_lpb) - mu * lip == 0

        beta * (lip + trS) == 1

        sum(alpha[i] for i = 1:6) == 1
        alpha[1] * v_rib - mu * trS == 0
        alpha[2] * v_rib - mu * prc == 0
        alpha[3] * v_rib - mu * catef == 0
        alpha[4] * v_rib - mu * metef == 0
        alpha[5] * v_rib - mu * rib == 0
        alpha[6] * v_rib - mu * lpb == 0

        (v_metef + gamm * v_catef - v_prc) - mu * atp == 0
        (v_prc - v_metef - gamm * v_catef) - mu * adp + mu * (atp + adp) == 0 # strange

        v_rib == (kcat_rib * rib * P) / (Km_rib + P)
        v_trS == (kcat_trS * trS * S) / (Km_trS + S)
        v_prc == (kcat_prc * prc * M * atp) / (M * atp + Km_prc * atp + Km_atp * M)
        v_catef ==
        (kcat_catef * catef * Si * adp) / (Si * adp + Si * Km_adp + adp * Km_catef)
        v_metef ==
        (kcat_metef * metef * Si * adp) / (Si * adp + Si * Km_adp + adp * Km_metef)
        v_lpb == (kcat_lpb * lpb * P) / (Km_lpb + P)

        rib + catef + metef + prc + lpb == 1 # maximal intracellular protein concentration
        atp + adp == 1 # total concentration of energy intermediates

        trS <= lip  # membrane integrity ((lip*trS) is the PLmax from the model)

    end

    @objective(molenaar, Max, mu)
    optimize!(molenaar)
    value(mu)

    res = Dict(
        "S" => value(S),
        "mu" => value(mu),
        "beta" => value(beta),
        "alpha" => value.(alpha),
        "v_rib" => value(v_rib),
        "v_trS" => value(v_trS),
        "v_catef" => value(v_catef),
        "v_metef" => value(v_metef),
        "v_lpb" => value(v_lpb),
        "atp" => value(atp),
        "adp" => value(adp),
    )
    return res
end

end #module


text:
@testset "Graph utils" begin
  # TODO
end