text:
"""
    ClusterAggregation

Functor object for using cluster aggregation in an L-shaped algorithm. Create by supplying a [`ClusterAggregate`](@ref) object through `aggregate ` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

The following cluster rules are available
- [`StaticCluster`](@ref)
- [`ClusterByReference`](@ref)
- [`Kmedoids`](@ref
- [`Hierarchical`](@ref)

...
# Parameters
- `rule::ClusterRule`: Rule that determines how cuts should be sorted into clusters
- `lock_after::Function = (τ,n)->false`: Function that determines if the current aggregation scheme should be fixed, based on the current optimality gap `τ` and the number of iterations `n`
...
"""
struct ClusterAggregation{T <: AbstractFloat, C <: AbstractClusterRule} <: AbstractAggregation
    buffer::Vector{SparseOptimalityCut{T}}
    rule::C
    partitioning::Dict{Int,Int}
    lock::Function
    aggregates::Vector{AggregatedOptimalityCut{T}}

    function ClusterAggregation(rule::AbstractClusterRule, lock::Function, ::Type{T}) where T <: AbstractFloat
        C = typeof(rule)
        aggregates = [zero(AggregatedOptimalityCut{T}) for _ = 1:5]
        return new{T,C}(Vector{SparseOptimalityCut{T}}(), rule, Dict{Int,Int}(), lock, aggregates)
    end
end

function aggregate_cut!(lshaped::AbstractLShaped, ::ClusterAggregation, cut::HyperPlane)
    return add_cut!(lshaped, cut)
end

function aggregate_cut!(lshaped::AbstractLShaped, aggregation::ClusterAggregation{T}, cut::HyperPlane{OptimalityCut}) where T <: AbstractFloat
    if aggregation.lock(gap(lshaped), num_iterations(lshaped)) && haskey(aggregation.partitioning, cut.id)
        aggregation.aggregates[aggregation.partitioning[cut.id]] += cut
        return false
    end
    push!(aggregation.buffer, cut)
    return false
end

function aggregate_cut!(cutqueue::CutQueue, aggregation::ClusterAggregation, ::MetaData, t::Integer, cut::HyperPlane, x::AbstractArray)
    put!(cutqueue, (t, cut))
    return nothing
end

function aggregate_cut!(cutqueue::CutQueue, aggregation::ClusterAggregation{T}, metadata::MetaData, t::Integer, cut::HyperPlane{OptimalityCut}, x::AbstractArray) where T <: AbstractFloat
    gap = fetch(metadata, t, :gap)
    if aggregation.lock(gap, t) && haskey(aggregation.partitioning, cut.id)
        aggregation.aggregates[aggregation.partitioning[cut.id]] += cut
        return nothing
    end
    push!(aggregation.buffer, cut)
    return nothing
end

function num_thetas(num_subproblems::Integer, ::ClusterAggregation)
    return num_subproblems
end

function num_thetas(num_subproblems::Integer, ::ClusterAggregation, ::AbstractScenarioProblems)
    return num_subproblems
end

function flush!(lshaped::AbstractLShaped, aggregation::ClusterAggregation{T}) where T <: AbstractFloat
    added = false
    if aggregation.lock(gap(lshaped), num_iterations(lshaped)) && !isempty(aggregation.partitioning)
        for (idx,aggregate) in enumerate(aggregation.aggregates)
            if !iszero(aggregate)
                added |= add_cut!(lshaped, aggregate)
                aggregation.aggregates[idx] = zero(AggregatedOptimalityCut{T})
            end
        end
    else
        if isempty(aggregation.buffer)
            return false
        end
        aggregates = cluster(aggregation.rule, aggregation.buffer)
        for (idx, aggregate) in enumerate(aggregates)
            for id in aggregate.ids
                aggregation.partitioning[id] = idx
            end
            added |= add_cut!(lshaped, aggregate)
        end
        n = maximum(values(aggregation.partitioning))
        resize!(aggregation.aggregates, n)
        for i in 1:n
            aggregation.aggregates[i] = zero(AggregatedOptimalityCut{T})
        end
    end
    empty!(aggregation.buffer)
    return added
end

function flush!(cutqueue::CutQueue, aggregation::ClusterAggregation{T}, metadata::MetaData, t::Integer, x::AbstractArray) where T <: AbstractFloat
    gap = fetch(metadata, t, :gap)
    if aggregation.lock(gap, t) && !isempty(aggregation.partitioning)
        for (idx,aggregate) in enumerate(aggregation.aggregates)
            if !iszero(aggregate)
                put!(cutqueue, (t, aggregate))
                aggregation.aggregates[idx] = zero(AggregatedOptimalityCut{T})
            end
        end
    else
        if isempty(aggregation.buffer)
            return nothing
        end
        aggregates = cluster(aggregation.rule, aggregation.buffer)
        for (idx, aggregate) in enumerate(aggregates)
            for id in aggregate.ids
                aggregation.partitioning[id] = idx
            end
            put!(cutqueue, (t, aggregate))
        end
        n = maximum(values(aggregation.partitioning))
        resize!(aggregation.aggregates, n)
        for i in 1:n
            aggregation.aggregates[i] = zero(AggregatedOptimalityCut{T})
        end
    end
    empty!(aggregation.buffer)
    return nothing
end

# API
# ------------------------------------------------------------
"""
    ClusterAggregate(rule::AbstractClusterRule; lock_after::Function = (τ,n)->false)

Factory object for [`ClusterAggregation`](@ref). Pass to `aggregate ` in the `LShapedSolver` factory function. See ?ClusterAggregation for parameter descriptions.

"""
mutable struct ClusterAggregate <: AbstractAggregator
    rule::AbstractClusterRule
    lock::Function

    function ClusterAggregate(rule::AbstractClusterRule; lock_after = (τ,n) -> false)
        return new(rule, lock_after)
    end
end

struct ClusterRule <: AggregationParameter end

function MOI.get(aggregator::ClusterAggregate, ::ClusterRule)
    return aggregator.rule
end

function MOI.set(aggregator::ClusterAggregate, ::ClusterRule, rule::AbstractClusterRule)
    return aggregator.rule = rule
end

function (aggregator::ClusterAggregate)(num_subproblems::Integer, T::Type{<:AbstractFloat})
    return ClusterAggregation(aggregator.rule, aggregator.lock, T)
end

function remote_aggregator(aggregation::ClusterAggregation, ::AbstractScenarioProblems, ::Integer)
    return ClusterAggregate(aggregation.rule; lock_after = aggregation.lock)
end

function str(aggregator::ClusterAggregate)
    return "cluster aggregation ruled by $(str(aggregator.rule))"
end


text:
# This file is included from within RadialBasisFunctionModels.jl #src

# ## Getting the Coefficients
const VecOfVecs{T} = AbstractVector{<:AbstractVector}

# ###  Polynomial Basis

# Any constructor of an `RBFModel` must solve for the coefficients in ``\eqref{eqn:coeff_basic}``.
# To build the equation system, we need a basis ``\{p_j\}_{1 \le j \le Q}`` of ``Π_d(ℝ^n)``.
# For the interpolation system to be solvable we have to choose the
# right polynomial space for ``p``.
# Basically, if the RBF Kernel (or the radial function) is
# *conditionally positive definite* of order ``D`` we have to
# find a polynomial ``p`` with ``\deg p \ge D-1``.[^wendland]
# If the kernel is CPD of order ``D=0`` we do not have to add an polynomial
# and can interpolate arbitrary (distinct) data points. \

# The canonical basis is ``x_1^{α_1} x_2^{α_2} … x_n^{α_n}`` with
# ``α_i ≥ 0`` and ``Σ_i α_i ≤ d``.
# For ``\bar{d} \le d`` we can recursively get the non-negative integer solutions for
# ``Σ_i α_i = \bar{d}`` with the following function:

@doc """
    non_negative_solutions( d :: Int, n :: Int)

Return a matrix with columns that correspond to solution vectors
``[x_1, …, x_n]`` to the equation ``x_1 + … + x_n = d``,
where the variables are non-negative integers.
"""
function non_negative_solutions( d :: Int, n :: Int ) :: Matrix{Int}
    if n == 1
        return fill(d,1,1)
    else
        num_sols = binomial( d + n - 1, n - 1)
        sol_matrix = Matrix{Int}(undef, n, num_sols)
        j = 1
        for i = 0 : d
            ## find all solutions of length `n-1` that sum to `i`
            ## if we add `d-i` to each column, then each column
            ## has `n` elements and sums to `d`
            padded_shorter_solutions = vcat( d-i, non_negative_solutions(i, n-1) )
            num_shorter_sols = size( padded_shorter_solutions, 2 )
            sol_matrix[:, j : j + num_shorter_sols - 1] .= padded_shorter_solutions
            j += num_shorter_sols
        end
        return sol_matrix
    end
end

# The polyonmial basis exponents are then given by all possible
# ``\bar{d}\le d``:
@doc """
    non_negative_solutions_ineq( d :: Int, n :: Int)

Return a matrix with columns that correspond to solution vectors
``[x_1, …, x_n]`` to the equation ``x_1 + … + x_n <= d``,
where the variables are non-negative integers.
"""
@memoize ThreadSafeDict function non_negative_solutions_ineq( d :: Int, n :: Int ) :: Matrix{Int}
    no_cols = binomial( n + d, n)
    ret_mat = Matrix{Int}(undef, n, no_cols )

    i = 1
    for d̄ = 0 : d
        sub_mat = non_negative_solutions( d̄, n)
        ret_mat[:, i:i + size(sub_mat,2) - 1] .= sub_mat[:,:]
        i += size(sub_mat, 2)
    end

    return ret_mat
end

# !!! note
#     I did an unnecessary rewrite of `non_negative_solutions` to be
#     Zygote-compatible. Therefore the matrices etc.
#     `Combinatorics` has `multiexponents` which should do the same...

# We **don't** use `DynamicPolynomials.jl` to generate the Polyomials **anymore**.
# Zygote did overflow when there were calculations with those polynomials.
# Not a problem for calculating the basis (because of we are `ignore()`ing
# the basis calculation now, assuming we never want to differentiate
# with respect to `n,d`),
# but when constructing the outputs from them.
# Instead we directly construct `StaticPolynomial`s and define a
# `PolynomialSystem` that evaluates all basis polynomials.

@doc """
    canonical_basis( n:: Int, d :: Int ) :: Union{PolynomialSystem, EmptyPolySystem}

Return the canonical basis of the space of `n`-variate
polynomials of degree at most `d`.
"""
function canonical_basis(n :: Int, d::Int, OneType :: Type = Float64)
    if d < 0
        return EmptyPolySystem{n}()
    else
        exponent_matrix = non_negative_solutions_ineq( d, n )
        one_float = OneType(1)  # `one_float` is used as coefficient(s) to guarantee floating point output
        return PolynomialSystem(
             ( Polynomial( [one_float,], e[:,:] ) for e ∈ eachcol(exponent_matrix) )...
        )
    end
end

# ### Solving the Equation System

# Now let ``\{p_j\}_{1\le j\le Q}`` be a basis of the polynomial space.
# Set ``P = [ p_j(x^i) ] ∈ ℝ^{N × Q}`` and ``Φ = φ(\| x^i - x^j \|)``.
# In case of interpolation, the linear equation system for the
# coefficients of $r$ is
# ```math
# S c := \begin{equation}
#     \begin{bmatrix}
#     Φ & P \\
#     P^T & 0_{Q × Q}
#     \end{bmatrix}
#     \begin{bmatrix}
#         w \\
#         λ
#     \end{bmatrix}
#     \stackrel{!}=
#     \begin{bmatrix}
#     Y
#     \\
#     0_Q
#     \end{bmatrix}.
#     \tag{I}
#     \label{eqn:coeff}
# \end{equation}
# ```

# We can also use differing feature vectors and centers.
# ``Φ`` then becomes
# ``Φ = [k_j(x^i)]_{1\le i \le N_d, 1\le j \le N_c} = [φ(‖ x^i - ξ^j ‖)]``,
# where we denote the number of kernel centers by ``N_c`` and the number
# of feauters (``d``ata) by ``N_d``.
# In the overdetermined least squares case (with pair-wise different centers and
# pair-wise different features), we do away with the second row of equations in \eqref{eqn:coeff}.
# The solution ``c = [w, λ]^T`` is then given by the
# Moore-Penrose Pseudo-Inverse:
# ```math
#     c = \underbrace{ ( S^T S )^{-1} S^T}_{=S^\dagger} \begin{bmatrix}
#     Y
#     \\
#     0_Q
#     \end{bmatrix}.
# ```
# Julia automatically computes the LS solution with `S\RHS`.

# !!! note
#     When we have vector data ``Y ⊂ ℝ^k``, e.g. from modelling MIMO functions, then
#     Julia easily allows for multiple columns in the righthand side of the interpolation
#     equation system and we get weight vectors for multiple models, that can
#     be thought of as one vector model ``r\colon ℝ^n \to ℝ^k``.

# Some helpers to build the matrices:

"Provided an array of sites with `n` variables, return a polynomial basis system of degree `poly_deg`."
function make_polys(sites, poly_deg = 1)
    poly_precision = promote_type(Float16, inner_type(sites))
    poly_basis_sys = Zyg.ignore() do
        canonical_basis( length(sites[1]), poly_deg, poly_precision )
    end
    return poly_basis_sys
end

function make_kernel( φ :: RadialFunction, center :: AbstractVector{<:Real})
    return ShiftedKernel(φ,center)
end

"Return array of `ShiftedKernel`s based functions in `φ_arr` with centers from `centers`."
function make_kernels( φ_arr :: AbstractVector{<:RadialFunction}, centers :: VecOfVecs )
    @assert length(φ_arr) == length(centers)
    return [ make_kernel(φ,c) for (φ,c) ∈ zip( φ_arr, centers) ]
end

"Return array of `ShiftedKernel`s based function `φ` with centers from `centers`."
function make_kernels( φ :: RadialFunction, centers :: VecOfVecs )
    return [ make_kernel(φ,c) for c ∈ centers ]
end

"Return RBF basis matrix by applying each kernel to each site. Kernels vary with the columns."
function _rbf_matrix( kernels, sites )
    return transpose( hcat( map(kernels, sites)... ) )
end

"Return polynomial basis matrix by applying each basis polynomial to each site. Polynomials vary with the columns."
function _poly_matrix( polys, sites )
    return transpose( hcat( map(polys, sites)... ) )
end

"Return RBF and Polynomial basis matrices as well as the vector of kernels and the polynomial basis system."
function get_matrices( φ, sites, centers = []; poly_deg = 1 )

    rbf_centers = isempty(centers) ? sites : centers

    kernels = make_kernels( φ, rbf_centers )
    polys = make_polys( sites, poly_deg )

    return _rbf_matrix( kernels, sites ), _poly_matrix( polys, sites ), kernels, polys
end

# Later, the construtor will call the above helper functions and then use
# these methods to retrieve the coefficients:

@doc """
    coefficients(sites, values, kernels, rad_funcs, polys )

Return the coefficient matrices `w` and `λ` for an rbf model
``r(x) = Σ_{i=1}^N wᵢ φ(\\|x - x^i\\|) + Σ_{j=1}^M λᵢ pᵢ(x)``,
where ``N`` is the number of centers and ``M``
is the number of n-variate basis polynomials.
"""
function coefficients( Φ, P, values;  mode :: Symbol = :ls)

    N_d, N_c = size( Φ )
    N_dp, Q = size( P )

    @assert N_d == N_dp "Row count of RBF and Poly basis matrices does not match."

    if N_d < N_c
        error("Underdetermined models not supported yet.")
    end
    if N_d < Q
        error("Too few data sites for selectod polynomial degree. (Need at least $(Q).)")
    end

    ## system matrix S and right hand side
    S = [Φ P]
    RHS = transpose( hcat(values... ) )

    return _coefficients( Φ, P, S, RHS, Val(:ls) )
end

# The actual work is delegated to `_coefficients`.

function _coefficients( Φ, P, S, RHS, ::Val{:ls} )
    N_c = size(Φ,2); Q = size(P,2);
    coeff = S \ RHS
    return _coeff_matrices(coeff, S, RHS, N_c, Q )
end

## interpolation requires zero padding
function _coefficients( Φ, P, S, RHS, ::Val{:interpolation} )
    N_d, N_c = size(Φ); Q = size(P,2);
    @assert N_d == N_c "Interpolation requires same number of features and centers." # TODO remove assertion
    S̃ = [ S ;                               # N_d × (N_c + Q)
          P' zeros(eltype(S), Q, Q )]       # Q × N_d and Q × Q
    RHS_padded = [ RHS;
        zeros( eltype(RHS), Q, size(RHS,2) ) ];
    coeff = S̃ \ RHS_padded
    return _coeff_matrices( coeff, S̃, RHS_padded, N_c, Q )
end

## treat sized arrays in a special way, so that zero padding preserves the types
## can potentially be removed once https://github.com/JuliaArrays/StaticArrays.jl/issues/856 is resolved
 function _coefficients( Φ, P, S :: StaticMatrix, RHS :: StaticMatrix, ::Val{:interpolation} )
    N_d, N_c = size(Φ); Q = size(P,2);
    @assert N_d == N_c "Interpolation requires same number of features and centers." # TODO remove assertion

    S̃ = [S ;
         P' @SMatrix(zeros(eltype(S),Q,Q)) ];
    RHS_padded = [ RHS;
        @SMatrix(zeros(eltype(RHS), Q ,size(RHS,2)))];
    coeff = S̃ \ RHS_padded
    return _coeff_matrices( coeff, S̃, RHS_padded, N_c, Q )
end

## this method returns random coefficients
function _coefficients( Φ, P, S, RHS, ::Val{:rand} )
    F = Base.promote_eltype(S, RHS)
    N_c = size(Φ,2); Q = size(P,2)
    m = N_c + Q
    k = size(RSH, 2)
    return _coeff_matrices(rand( F, m, k ), S, RHS, N_c, Q )
end

# We want statically sized matrices when appropriate, that is why we
# call `_coeff_matrices`:
"Return rbf weights, polynomial regressor weights, system matrix and right hand side of the RBF equations."
function _coeff_matrices(coeff :: AbstractMatrix, S, RHS, N_c, Q )
    return view(coeff, 1 : N_c, :), view(coeff, N_c + 1 : N_c + Q, :), S, RHS
end

## index into `coeff` using static array so that the results are sized, too
function _coeff_matrices(coeff :: StaticMatrix, S, RHS, N_c, Q )
    return coeff[ SVector{N_c}(1 : N_c), :], coeff[ SVector{Q}( N_c + 1 : N_c + Q ), :], S, RHS
end

# We can easily impose linear equality constraints,
# for example requiring interpolation only on a subset of features.
# In matrix form, $I$ linear equality constraints (for ``k`` outputs) can be written
# as
# ```math
# E c = b, \quad E ∈ ℝ^{I×(N_c + Q)}, b ∈ ℝ^{I\times k},\, I,k ∈ ℕ_0.
# ```
# Now, let $ĉ$ be the least squares solution from above.
# The constrained solution is
# ```math
#  c = ĉ - Ψ E^T ( E Ψ E^T)^{-1} ( E ĉ - b ), \; Ψ := (S^T S)^{-1}
# \tag{CLS1}
# \label{eqn:cls1}
# ```
# This results from forming the Lagrangian of an equivalent minimization problem.
# Let ``δ = ĉ - c ∈ ℝ^{q\times k}, q = N_c + Q,`` and define the constraint residuals
# as ``γ = Eĉ - b ∈ ℝ^{I\times k}``.
# The Lagrangian for minimizing ``δ^TS^TSδ`` under ``Eδ=γ`` is
# ```math
# \begin{aligned}
#     L &= δ^T S^T S δ + 2 λ^T( E δ - γ )\\
#     D_δL &= 2 δ^T S^T S + 2λ^T E \\
#     D_λL &= 2 δ^T E^T - 2 γ^T
# \end{aligned}
# ```
# Setting the derivatives to zero leads to \eqref{eqn:cls1} via
# ```math
#     \begin{bmatrix}
#         S^T S & E^T \\
#         E & 0_{I\times I}
#     \end{bmatrix}
#     \begin{bmatrix}
#     δ \\ λ
#     \end{bmatrix}
#     = \begin{bmatrix}
#     0_{q\times k} \\ γ
#     \end{bmatrix}
# \tag{L}
# \label{eqn:cls2}
# ```
# See [^adv_eco] for details.

# TODO make this respect StaticMatrices too!! #src
function constrained_coefficients(
        w :: AbstractMatrix{<:Real},
        λ :: AbstractMatrix{<:Real},
        S :: AbstractMatrix{<:Real},
        E :: AbstractMatrix{<:Real},
        b :: AbstractMatrix{<:Real}
    )
    ## Using Lagrangian approach:

    ĉ = [w; λ]  # least squares solution
    γ = E*ĉ - b # constraint residuals

    I, q = size(E)
    k = size(w,2)

    A = vcat(
        [S'S E'],
        [E zeros(Int,I,I)]
    )

    RHS = [
        zeros(Int, q, k);
        γ
    ]

    δλ = A \ RHS
    δ = δλ[1 : q, :]

    c = ĉ - δ  # coefficients for constrained problem

    N_c = size(w,1)

    return c[1 : N_c, :], c[N_c+1:end, :]
end

# For the case that mentioned above, that is, interpolation at a
# subset of sites, we can easily build the ``E`` matrix from the ``S``
# matrix by taking the corresponding rows.
function constrained_coefficients(
        w :: AbstractMatrix{<:Real},
        λ :: AbstractMatrix{<:Real},
        S :: AbstractMatrix{<:Real},
        RHS_ls :: AbstractMatrix{<:Real},
        interpolation_indices :: AbstractVector{Int}
    )

    E = S[interpolation_indices, :]
    b = RHS_ls[interpolation_indices, :]
    return constrained_coefficients( w, λ, S, E, b )
end

# ### The Actual, Usable Constructor

# We want the user to be able to pass 1D data as scalars and use the following helpers:

ensure_vec_of_vecs( before :: AbstractVector{<:AbstractVector{<:Real}} ) = before
ensure_vec_of_vecs( before :: AbstractVector{ <:Real }) = [[x,] for x in before ]

function inner_type( vec_of_vecs :: AbstractVector{<:AbstractVector{T}}) where T
    if Base.isabstracttype(T)   # like Any if data is of mixed precision
        return Float64
    else
        return T
    end
end

function _check_data(features,labels,centers; type_checks = true)
    ## Basic Data integrity checks

    @assert !isempty(features) "Provide at least 1 feature vector."
    @assert !isempty(labels) "Provide at least 1 label vector."
    num_vars = length(features[1])
    num_outputs = length(labels[1])
    @assert all( length(s) == num_vars for s ∈ features ) "All features must have same dimension."
    @assert all( length(v) == num_outputs for v ∈ labels ) "All labels must have same dimension."

    num_sites = length(features)
    num_vals = length(labels)
    @assert num_sites == num_vals "Provide as many features as labels."

    if type_checks
        sites = ensure_vec_of_vecs(features)
        values = ensure_vec_of_vecs(labels)
        if !isempty(centers)
            @assert all( length(c) == num_vars for c ∈ centers ) "All centers must have dimension $(num_vars)."
            C = ensure_vec_of_vecs(centers)
        else
            C = sites
        end
    else
        sites = features
        values = labels
        C = centers
    end

    num_centers = length(C)

    return num_vars, num_outputs, num_sites, num_vals, num_centers, sites, values, C
end

# We now have all ingredients for the basic outer constructor:

@doc """
    RBFModel( features, labels, φ = Multiquadric(), poly_deg = 1; kwargs ... )

Construct a `RBFModel` from the feature vectors in `features` and
the corresponding labels in `lables`, where `φ` is a `RadialFunction` or a vector of
`RadialFunction`s.\n
Scalar data can be used, it is transformed internally. \n
StaticArrays can be used, e.g., `features :: Vector{<:SVector}`.
Providing `SVector`s only might speed up the construction.\n
If the degree of the polynomial tail, `poly_deg`, is too small it will be set to `cpd_order(φ)-1`.

If the RBF centers do not equal the the `features`, you can use the keyword argument `centers` to
pass a list of centers. If `φ` is a vector, then the length of `centers` and `φ` must be equal and
`centers[i]` will be used in conjunction with `φ[i]` to build a `ShiftedKernel`. \n
If `features` has 1D data, the output of the model will be a 1D-vector.
If it should be a scalar instead, set the keyword argument `vector_output` to `false`.
"""
function RBFModel(
        features :: AbstractVector{ <:NumberOrVector },
        labels :: AbstractVector{ <:NumberOrVector },
        φ :: Union{RadialFunction, AbstractVector{<:RadialFunction}} = Multiquadric(),
        poly_deg :: Int = 1;
        centers :: AbstractVector{ <:NumberOrVector } = Vector{Float16}[],
        interpolation_indices :: AbstractVector{ <: Int } = Int[],
        vector_output :: Bool = true,
        coeff_mode :: Symbol = :auto,
        save_matrices :: Bool = true
    )

    num_vars, num_outputs, num_sites, num_vals, num_centers, sites, values, C = _check_data( features, labels, centers )

    if coeff_mode == :auto
        can_interpolate_uniquely = φ isa RadialFunction ? poly_deg >= cpd_order(φ) - 1 : all( poly_deg >= cpd_order(phi) - 1 for phi in φ )
        coeff_mode = num_sites == num_centers && can_interpolate_uniquely ? :interpolation : :ls
    end

    Φ, P, kernels, poly_basis_sys = get_matrices( φ, sites, C; poly_deg )

    w, λ, S, RHS = coefficients( Φ, P, values; mode = coeff_mode )

    meta = save_matrices ? (rbf_mat = Φ, poly_mat = P) : nothing

    if !isempty(interpolation_indices)
        w, λ = constrained_coefficients( w, λ, S, RHS, interpolation_indices)
    end

    ## build output polynomials
    poly_sum = PolySum( poly_basis_sys, transpose(λ) )

    ## build RBF system
    rbf_sys = RBFSum(kernels, transpose(w), num_outputs)

    ## vector output? (dismiss user choice if labels are vectors)
    vec_output = num_outputs == 1 ? vector_output : true

    return RBFModel(rbf_sys, poly_sum, num_vars, num_outputs, num_centers, meta; vec_output)
end

# ### Special Constructors

# We offer some specialized models (that simply wrap the main type).
struct RBFInterpolationModel
    model :: RBFModel
end
(mod :: RBFInterpolationModel)(args...) = mod.model(args...)
@forward RBFInterpolationModel.model grad, jac, jacT, auto_grad, auto_jac

# The constructor is a tiny bit simpler and additional checks take place:
"""
    RBFInterpolationModel(features, labels, φ, poly_deg; kwargs… )

Build a model interpolating the feature-label pairs.
Does not accept `center` keyword argument.
"""
function RBFInterpolationModel(
        features :: AbstractVector{ <:NumberOrVector },
        labels :: AbstractVector{ <:NumberOrVector },
        φ :: Union{RadialFunction,AbstractVector{<:RadialFunction}} = Multiquadric(),
        poly_deg :: Int = 1;
        vector_output :: Bool = true,
        save_matrices :: Bool = true
    )
    @assert length(features) == length(labels) "Provide as many features as labels!"

    if poly_deg < cpd_order(φ) - 1
        @warn "Polyonmial degree too small for interpolation. Using $(cpd_order(φ)-1)."
        poly_deg = max( poly_deg,  cpd_order(φ) - 1 )
    end

    mod = RBFModel(features, labels, φ, poly_deg; vector_output, save_matrices, coeff_mode = :interpolation)
    return RBFInterpolationModel( mod )
end

# We want to provide a convenient alternative constructor for interpolation models
# so that the radial function can be defined by passing a `Symbol` or `String`.

const SymbolToRadialConstructor = NamedTuple((
    :gaussian => Gaussian,
    :multiquadric => Multiquadric,
    :inv_multiquadric => InverseMultiquadric,
    :cubic => Cubic,
    :thin_plate_spline => ThinPlateSpline
))

"Obtain a `RadialFunction` from its name and constructor arguments."
function _get_rad_func( φ_symb :: Union{Symbol, String}, φ_args )

    ## which radial function to use?
    radial_symb = Symbol( lowercase( string( φ_symb ) ) )
    if !(radial_symb ∈ keys(SymbolToRadialConstructor))
        @warn "Radial Funtion $(radial_symb) not known, using Gaussian."
        radial_symb = :gaussian
        φ_args = nothing
    end

    constructor = SymbolToRadialConstructor[radial_symb]
    if isnothing(φ_args)
        φ = constructor()
    else
        φ = constructor( φ_args... )
    end

    return φ
end

# The alternative constructors are build programmatically:
for op ∈ [ :RBFInterpolationModel, :RBFModel ]
    @eval begin
        function $op(
                features :: AbstractVector{ <:NumberOrVector },
                labels :: AbstractVector{ <:NumberOrVector },
                φ_symb :: Union{Symbol, String},
                φ_args = nothing,
                poly_deg :: Int = 1; kwargs...
            )

            φ = _get_rad_func( φ_symb, φ_args )
            return $op(features, labels, φ, poly_deg; kwargs... )
        end
    end
end

# ### Container with Training Data

# The RBF Machine is similar in design to what an MLJ machine does:
# Training data (feature and label **vectors**) are stored and can be added.
# The inner model is trained with `fit!`.

# **TODO** In the future, we can customize the `fit!` method when updating a model
# to only consider *new* training data.
# This also makes type conversion of the whole data arrays unnecessary.

"""
    RBFMachineWithKernel(; φ, features, labels, poly_deg)

A container holding an inner `model :: RBFModel` (or `model == nothing`).
The inner model interpolates the data after calling `fit!`.
An array of arrays of features is stored in the `features` field.
Likewise for `labels`.
After `add_data!`, the model has to be fitted again, indicated by its field `is_valid`.

In constrast to the more general RBFModel, only one single RBF function (i.e. one set of shape parameters)
should be used.
"""
@with_kw mutable struct RBFMachineWithKernel{
        FT <: AbstractVector{<:AbstractVector{<:AbstractFloat}},
        LT <: AbstractVector{<:AbstractVector{<:AbstractFloat}},
        R <: RadialFunction,
    }

    φ :: R = Gaussian()

    features :: FT = Vector{Float64}[]
    labels :: LT = Vector{Float64}[]

    poly_deg :: Int = 1

    model :: Union{Nothing, RBFModel} = nothing
    valid :: Bool = false   # is model trained on all data sites?

    function RBFMachineWithKernel{FT,LT,R}(φ::R, features::FT, labels::LT, poly_deg, model, valid ) where{R,FT,LT}
        @assert length(features) == length(labels) "Need same number of `features` and `labels`."
        features_copied = copy(features)
        labels_copied = copy(labels)
        return new{FT,LT,R}(φ,features_copied, labels_copied, poly_deg, model, valid )
    end

end

function RBFMachine(; kernel_name :: Symbol = :gaussian, kernel_args = nothing, poly_deg = 1,
    labels :: FT = Vector{Float64}[], features :: LT = Vector{Float64}[], kwargs... ) where {FT, LT}
    T = promote_type(eltype(eltype(FT)), eltype(eltype(LT)))
    KARGS = isnothing(kernel_args) ? nothing : T.(kernel_args)
    φ = _get_rad_func( kernel_name, KARGS )
    @assert poly_deg >= cpd_order(φ) - 1 "Polynomial degree too low for interpolation."
    return RBFMachineWithKernel(; φ, poly_deg, labels, features, kwargs...)
end

"Return floating point type of training data elements."
_precision( :: RBFMachineWithKernel{FT,LT} ) where {FT,LT} = eltype( Base.promote_eltype(FT, LT) )

"Fit `mach :: RBFMachineWithKernel` to the training data."
function fit!( mach :: RBFMachineWithKernel )::Nothing
    @assert length(mach.features) > 0 "Cannot `fit!` without data."
    @assert length(mach.features) == length(mach.labels) "Interpolation requires same number of features and labels."
    num_needed =  binomial( mach.poly_deg + length(mach.features[1]), mach.poly_deg)
    @assert length(mach.features) >= num_needed "Too few data sites for selected polynomial degree (need $(num_needed))."

    mach.model = RBFModel( mach.features, mach.labels, mach.φ, mach.poly_deg; save_matrices = true )

    mach.valid = true
    return nothing
end

# Forward evaluation methods of inner model:
( mach :: RBFMachineWithKernel )(args...) = mach.model(args...)
@forward RBFMachineWithKernel.model grad, jac, jacT, auto_grad, auto_jac

# Methods to add features and labels:
"Add a feature vector(s) and a label(s) to the `machine` container."
function add_data!(
        m :: RBFMachineWithKernel, features :: AbstractVector{<:AbstractVector}, labels :: AbstractVector{<:AbstractVector}
    ) :: Nothing
    @assert length(features) == length(labels) "Provide same number of features and labels."
    @assert all( length(f) == length(features[1]) for f in features ) "Features must have same length."
    @assert all( length(l) == length(labels[1]) for l in labels ) "Labels must have same length"
    @assert isempty(m.features) || length(m.features[1]) == length(features[1]) && length(m.labels[1]) == length(labels[1]) "Length doesnt match previous data."
    append!(m.features, features)
    append!(m.labels, labels)
    m.valid = false
    return nothing
end

function add_data!(
        m :: RBFMachineWithKernel, feature :: AbstractVector{<:AbstractFloat}, label:: AbstractVector{<:AbstractFloat}
    ) :: Nothing
    return add_data!(m, [ feature, ], [label, ])
end

# Convenience methods to "reset" a machine:

function Base.empty!( m :: RBFMachineWithKernel ) :: Nothing
    empty!(m.features)
    empty!(m.labels)
    m.model = nothing
    m.valid = false
    return nothing
end

function Base.isempty(m :: RBFMachineWithKernel ) :: Bool
    isempty( m.features ) && isempty( m.labels ) && isnothing(m.model)
end

ext: jl


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

module ResolveTest
import ..Pkg # ensure we are using the correct Pkg

using Test
using Pkg.Types
using Pkg.Types: VersionBound
using UUIDs
using Pkg.Resolve
import Pkg.Resolve: VersionWeight, add_reqs!, simplify_graph!, ResolverError, Fixed, Requires

# print info, stats etc.
const VERBOSE = false

# Check that VersionWeight keeps the same ordering as VersionNumber

vlst = [
    v"0.0.0",
    v"0.0.1",
    v"0.1.0",
    v"0.1.1",
    v"1.0.0",
    v"1.0.1",
    v"1.1.0",
    v"1.1.1"
    ]

for v1 in vlst, v2 in vlst
    vw1 = VersionWeight(v1)
    vw2 = VersionWeight(v2)
    clt = v1 < v2
    @test clt == (vw1 < vw2)
    ceq = v1 == v2
    @test ceq == (vw1 == vw2)
end

# auxiliary functions
const uuid_package = UUID("cfb74b52-ec16-5bb7-a574-95d9e393895e")
pkguuid(p::String) = uuid5(uuid_package, p)
function storeuuid(p::String, uuid_to_name::Dict{UUID,String})
    uuid = p == "julia" ? Resolve.uuid_julia : pkguuid(p)
    if haskey(uuid_to_name, uuid)
        @assert uuid_to_name[uuid] == p
    else
        uuid_to_name[uuid] = p
    end
    return uuid
end
wantuuids(want_data) = Dict{UUID,VersionNumber}(pkguuid(p) => v for (p,v) in want_data)

function graph_from_data(deps_data)
    uuid_to_name = Dict{UUID,String}()
    uuid(p) = storeuuid(p, uuid_to_name)
    fixed = Dict{UUID,Fixed}()
    all_versions = Dict{UUID,Set{VersionNumber}}()
    all_deps = Dict{UUID,Dict{VersionNumber,Dict{String,UUID}}}()
    all_compat = Dict{UUID,Dict{VersionNumber,Dict{String,VersionSpec}}}()

    deps = Dict{String,Dict{VersionNumber,Dict{String,VersionSpec}}}()
    for d in deps_data
        p, vn, r = d[1], d[2], d[3:end]
        if !haskey(deps, p)
            deps[p] = Dict{VersionNumber,Dict{String,VersionSpec}}()
        end
        if !haskey(deps[p], vn)
            deps[p][vn] = Dict{String,VersionSpec}()
        end
        isempty(r) && continue
        rp = r[1]
        rvs = VersionSpec(r[2:end])
        deps[p][vn][rp] = rvs
    end
    for (p,preq) in deps
        u = uuid(p)
        all_versions[u] = Set(keys(deps[p]))

        deps_pkgs = Dict{String,Set{VersionNumber}}()
        for (vn,vreq) in deps[p], rp in keys(vreq)
            push!(get!(Set{VersionNumber}, deps_pkgs, rp), vn)
        end
        all_deps[u] = Dict{VersionNumber,Dict{String,UUID}}()
        all_compat[u] = Dict{VersionNumber,Dict{String,VersionSpec}}()
        for (vn,vreq) in preq
            all_deps[u][vn] = Dict{String,UUID}()
            all_compat[u][vn] = Dict{String,VersionSpec}()
            for (rp,rvs) in vreq
                all_deps[u][vn][rp] = uuid(rp)
                all_compat[u][vn][rp] = rvs
            end
        end
    end
    return Graph(all_versions, all_deps, all_compat, uuid_to_name, Requires(), fixed, VERBOSE)
end
function reqs_from_data(reqs_data, graph::Graph)
    reqs = Dict{UUID,VersionSpec}()
    function uuid_check(p)
        uuid = pkguuid(p)
        @assert graph.data.uuid_to_name[uuid] == p
        return uuid
    end
    for r in reqs_data
        p = uuid_check(r[1])
        reqs[p] = VersionSpec(r[2:end])
    end
    reqs
end
function sanity_tst(deps_data, expected_result; pkgs=[])
    if VERBOSE
        println()
        @info("sanity check")
        # @show deps_data
        # @show pkgs
    end
    graph = graph_from_data(deps_data)
    id(p) = pkgID(pkguuid(p), graph)
    result = sanity_check(graph, Set(pkguuid(p) for p in pkgs), VERBOSE)

    length(result) == length(expected_result) || return false
    expected_result_uuid = [(id(p), vn) for (p,vn) in expected_result]
    for r in result
        r ∈ expected_result_uuid || return  false
    end
    return true
end
sanity_tst(deps_data; kw...) = sanity_tst(deps_data, []; kw...)

function resolve_tst(deps_data, reqs_data, want_data = nothing; clean_graph = false)
    if VERBOSE
        println()
        @info("resolving")
        # @show deps_data
        # @show reqs_data
    end
    graph = graph_from_data(deps_data)
    reqs = reqs_from_data(reqs_data, graph)
    add_reqs!(graph, reqs)
    simplify_graph!(graph, clean_graph = clean_graph)
    want = resolve(graph)
    return want == wantuuids(want_data)
end

@testset "schemes" begin
    VERBOSE && @info("SCHEME 1")
    ## DEPENDENCY SCHEME 1: TWO PACKAGES, DAG
    deps_data = Any[
        ["A", v"1", "B", "1-*"],
        ["A", v"2", "B", "2-*"],
        ["B", v"1"],
        ["B", v"2"]
    ]

    @test sanity_tst(deps_data)
    @test sanity_tst(deps_data, pkgs=["A", "B"])
    @test sanity_tst(deps_data, pkgs=["B"])
    @test sanity_tst(deps_data, pkgs=["A"])

    # require just B
    reqs_data = Any[
        ["B", "*"]
    ]

    want_data = Dict("B"=>v"2")
    resolve_tst(deps_data, reqs_data, want_data)
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require just A: must bring in B
    reqs_data = Any[
        ["A", "*"]
    ]
    want_data = Dict("A"=>v"2", "B"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    VERBOSE && @info("SCHEME 2")
    ## DEPENDENCY SCHEME 2: TWO PACKAGES, CYCLIC
    deps_data = Any[
        ["A", v"1", "B", "2-*"],
        ["A", v"2", "B", "1-*"],
        ["B", v"1", "A", "2-*"],
        ["B", v"2", "A", "1-*"]
    ]

    @test sanity_tst(deps_data)

    # require just A
    reqs_data = Any[
        ["A", "*"]
    ]
    want_data = Dict("A"=>v"2", "B"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require just B, force lower version
    reqs_data = Any[
        ["B", "1"]
    ]
    want_data = Dict("A"=>v"2", "B"=>v"1")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require just A, force lower version
    reqs_data = Any[
        ["A", "1"]
    ]
    want_data = Dict("A"=>v"1", "B"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    VERBOSE && @info("SCHEME 3")
    ## DEPENDENCY SCHEME 3: THREE PACKAGES, CYCLIC, TWO MUTUALLY EXCLUSIVE SOLUTIONS
    deps_data = Any[
        ["A", v"1", "B", "2-*"],
        ["A", v"2", "B", "1"],
        ["B", v"1", "C", "2-*"],
        ["B", v"2", "C", "1"],
        ["C", v"1", "A", "1"],
        ["C", v"2", "A", "2-*"]
    ]

    @test sanity_tst(deps_data)

    # require just A (must choose solution which has the highest version for A)
    reqs_data = Any[
        ["A", "*"]
    ]
    want_data = Dict("A"=>v"2", "B"=>v"1", "C"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require just B (must choose solution which has the highest version for B)
    reqs_data = Any[
        ["B", "*"]
    ]
    want_data = Dict("A"=>v"1", "B"=>v"2", "C"=>v"1")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require just A, force lower version
    reqs_data = Any[
        ["A", "1"]
    ]
    want_data = Dict("A"=>v"1", "B"=>v"2", "C"=>v"1")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require A and C, incompatible versions
    reqs_data = Any[
        ["A", "1"],
        ["C", "2-*"]
    ]
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)

    VERBOSE && @info("SCHEME 4")
    ## DEPENDENCY SCHEME 4: TWO PACKAGES, DAG, WITH TRIVIAL INCONSISTENCY
    deps_data = Any[
        ["A", v"1", "B", "2-*"],
        ["B", v"1"]
    ]

    @test sanity_tst(deps_data, [("A", v"1")])
    @test sanity_tst(deps_data, pkgs=["B"])

    # require B (must not give errors)
    reqs_data = Any[
        ["B", "*"]
    ]
    want_data = Dict("B"=>v"1")
    @test resolve_tst(deps_data, reqs_data, want_data)

    VERBOSE && @info("SCHEME 5")
    ## DEPENDENCY SCHEME 5: THREE PACKAGES, DAG, WITH IMPLICIT INCONSISTENCY
    deps_data = Any[
        ["A", v"1", "B", "2-*"],
        ["A", v"1", "C", "2-*"],
        # ["A", v"1", "julia", "10"],
        ["A", v"2", "B", "1"],
        ["A", v"2", "C", "1"],
        ["B", v"1", "C", "2-*"],
        ["B", v"2", "C", "2-*"],
        ["C", v"1"],
        ["C", v"2"]
    ]

    @test sanity_tst(deps_data, [("A", v"2")])
    @test sanity_tst(deps_data, pkgs=["B"])
    @test sanity_tst(deps_data, pkgs=["C"])

    # require A, any version (must use the highest non-inconsistent)
    reqs_data = Any[
        ["A", "*"]
    ]
    want_data = Dict("A"=>v"1", "B"=>v"2", "C"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require A, force highest version (impossible)
    reqs_data = Any[
        ["A", "2-*"]
    ]
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)

    VERBOSE && @info("SCHEME 6")
    ## DEPENDENCY SCHEME 6: TWO PACKAGES, CYCLIC, TOTALLY INCONSISTENT
    deps_data = Any[
        ["A", v"1", "B", "2-*"],
        ["A", v"2", "B", "1"],
        ["B", v"1", "A", "1"],
        ["B", v"2", "A", "2-*"]
    ]

    @test sanity_tst(deps_data, [("A", v"1"), ("A", v"2"),
                                ("B", v"1"), ("B", v"2")])

    # require A (impossible)
    reqs_data = Any[
        ["A", "*"]
    ]
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)

    # require B (impossible)
    reqs_data = Any[
        ["B", "*"]
    ]
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)

    VERBOSE && @info("SCHEME 7")
    ## DEPENDENCY SCHEME 7: THREE PACKAGES, CYCLIC, WITH INCONSISTENCY
    deps_data = Any[
        ["A", v"1", "B", "1"],
        ["A", v"2", "B", "2-*"],
        ["B", v"1", "C", "1"],
        ["B", v"2", "C", "2-*"],
        ["C", v"1", "A", "2-*"],
        ["C", v"2", "A", "2-*"],
    ]

    @test sanity_tst(deps_data, [("A", v"1"), ("B", v"1"),
                                ("C", v"1")])

    # require A
    reqs_data = Any[
        ["A", "*"]
    ]
    want_data = Dict("A"=>v"2", "B"=>v"2", "C"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require C
    reqs_data = Any[
        ["C", "*"]
    ]
    want_data = Dict("A"=>v"2", "B"=>v"2", "C"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require C, lowest version (impossible)
    reqs_data = Any[
        ["C", "1"]
    ]
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)

    VERBOSE && @info("SCHEME 8")
    ## DEPENDENCY SCHEME 8: THREE PACKAGES, CYCLIC, TOTALLY INCONSISTENT
    deps_data = Any[
        ["A", v"1", "B", "1"],
        ["A", v"2", "B", "2-*"],
        ["B", v"1", "C", "1"],
        ["B", v"2", "C", "2-*"],
        ["C", v"1", "A", "2-*"],
        ["C", v"2", "A", "1"],
    ]

    @test sanity_tst(deps_data, [("A", v"1"), ("A", v"2"),
                                ("B", v"1"), ("B", v"2"),
                                ("C", v"1"), ("C", v"2")])

    # require A (impossible)
    reqs_data = Any[
        ["A", "*"]
    ]
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)

    # require B (impossible)
    reqs_data = Any[
        ["B", "*"]
    ]
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)

    # require C (impossible)
    reqs_data = Any[
        ["C", "*"]
    ]
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)

    VERBOSE && @info("SCHEME 9")
    ## DEPENDENCY SCHEME 9: SIX PACKAGES, DAG
    deps_data = Any[
        ["A", v"1"],
        ["A", v"2"],
        ["A", v"3"],
        ["B", v"1", "A", "1"],
        ["B", v"2", "A", "*"],
        ["C", v"1", "A", "2"],
        ["C", v"2", "A", "2-*"],
        ["D", v"1", "B", "1-*"],
        ["D", v"2", "B", "2-*"],
        ["E", v"1", "D", "*"],
        ["F", v"1", "A", "1-2"],
        ["F", v"1", "E", "*"],
        ["F", v"2", "C", "2-*"],
        ["F", v"2", "E", "*"],
    ]

    @test sanity_tst(deps_data)

    # require just F
    reqs_data = Any[
        ["F", "*"]
    ]
    want_data = Dict("A"=>v"3", "B"=>v"2", "C"=>v"2",
                    "D"=>v"2", "E"=>v"1", "F"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require just F, lower version
    reqs_data = Any[
        ["F", "1"]
    ]
    want_data = Dict("A"=>v"2", "B"=>v"2", "D"=>v"2",
                    "E"=>v"1", "F"=>v"1")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require F and B; force lower B version -> must bring down F, A, and D versions too
    reqs_data = Any[
        ["F", "*"],
        ["B", "1"]
    ]
    want_data = Dict("A"=>v"1", "B"=>v"1", "D"=>v"1",
                    "E"=>v"1", "F"=>v"1")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require F and D; force lower D version -> must not bring down F version
    reqs_data = Any[
        ["F", "*"],
        ["D", "1"]
    ]
    want_data = Dict("A"=>v"3", "B"=>v"2", "C"=>v"2",
                    "D"=>v"1", "E"=>v"1", "F"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require F and C; force lower C version -> must bring down F and A versions
    reqs_data = Any[
        ["F", "*"],
        ["C", "1"]
    ]
    want_data = Dict("A"=>v"2", "B"=>v"2", "C"=>v"1",
                    "D"=>v"2", "E"=>v"1", "F"=>v"1")
    @test resolve_tst(deps_data, reqs_data, want_data)

    VERBOSE && @info("SCHEME 10")
    ## DEPENDENCY SCHEME 10: FIVE PACKAGES, SAME AS SCHEMES 5 + 1, UNCONNECTED
    deps_data = Any[
        ["A", v"1", "B", "2-*"],
        ["A", v"1", "C", "2-*"],
        ["A", v"2", "B", "1"],
        ["A", v"2", "C", "1"],
        ["B", v"1", "C", "2-*"],
        ["B", v"2", "C", "2-*"],
        ["C", v"1"],
        ["C", v"2"],
        ["D", v"1", "E", "1-*"],
        ["D", v"2", "E", "2-*"],
        ["E", v"1"],
        ["E", v"2"]
    ]

    @test sanity_tst(deps_data, [("A", v"2")])
    @test sanity_tst(deps_data, pkgs=["B"])
    @test sanity_tst(deps_data, pkgs=["D"])
    @test sanity_tst(deps_data, pkgs=["E"])
    @test sanity_tst(deps_data, pkgs=["B", "D"])

    # require A, any version (must use the highest non-inconsistent)
    reqs_data = Any[
        ["A", "*"]
    ]
    want_data = Dict("A"=>v"1", "B"=>v"2", "C"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require just D: must bring in E
    reqs_data = Any[
        ["D", "*"]
    ]
    want_data = Dict("D"=>v"2", "E"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)

    # require A and D, must be the merge of the previous two cases
    reqs_data = Any[
        ["A", "*"],
        ["D", "*"]
    ]
    want_data = Dict("A"=>v"1", "B"=>v"2", "C"=>v"2", "D"=>v"2", "E"=>v"2")
    @test resolve_tst(deps_data, reqs_data, want_data)
end

@testset "realistic" begin
    VERBOSE && @info("SCHEME REALISTIC")
    ## DEPENDENCY SCHEME 11: A REALISTIC EXAMPLE
    ## ref Julia issue #21485

    include("resolvedata1.jl")

    @test sanity_tst(ResolveData.deps_data, ResolveData.problematic_data)
    @test resolve_tst(ResolveData.deps_data, ResolveData.reqs_data, ResolveData.want_data)
end

@testset "nasty" begin
    VERBOSE && @info("SCHEME NASTY")
    ## DEPENDENCY SCHEME 12: A NASTY CASE

    include("NastyGenerator.jl")
    deps_data, reqs_data, want_data, problematic_data = NastyGenerator.generate_nasty(5, 20, q=20, d=4, sat = true)

    @test sanity_tst(deps_data, problematic_data)
    @test resolve_tst(deps_data, reqs_data, want_data)

    deps_data, reqs_data, want_data, problematic_data = NastyGenerator.generate_nasty(5, 20, q=20, d=4, sat = false)

    @test sanity_tst(deps_data, problematic_data)
    @test_throws ResolverError resolve_tst(deps_data, reqs_data)
end

end # module


text:
# convenience functions for simulation studies

#####################################################################################################
# change observed (data) without reconstructing the whole model
#####################################################################################################

# use the same observed type as before
swap_observed(model::AbstractSemSingle; kwargs...) = swap_observed(model, typeof(observed(model)).name.wrapper; kwargs...)

# construct a new observed type
swap_observed(model::AbstractSemSingle, observed_type; kwargs...) = swap_observed(model, observed_type(;kwargs...); kwargs...)

swap_observed(model::AbstractSemSingle, new_observed::SemObs; kwargs...) =
    swap_observed(model, observed(model), imply(model), loss(model), diff(model), new_observed; kwargs...)

function swap_observed(
        model::AbstractSemSingle,
        old_observed,
        imply,
        loss,
        diff,
        new_observed::SemObs;
        kwargs...)

    kwargs = Dict{Symbol, Any}(kwargs...)

    # get field types
    kwargs[:observed_type] = typeof(new_observed)
    kwargs[:old_observed_type] = typeof(old_observed)
    kwargs[:imply_type] = typeof(imply)
    kwargs[:loss_types] = [typeof(lossfun) for lossfun in loss.functions]
    kwargs[:diff_type] = typeof(diff)

    # update imply
    imply = update_observed(imply, new_observed; kwargs...)
    kwargs[:imply] = imply
    kwargs[:n_par] = n_par(imply)

    # update loss
    loss = update_observed(loss, new_observed; kwargs...)
    kwargs[:loss] = loss

    # update diff
    diff = update_observed(diff, new_observed; kwargs...)

    #new_imply = update_observed(model.imply, new_observed; kwargs...)

    return Sem(
        new_observed,
        update_observed(model.imply, new_observed; kwargs...),
        update_observed(model.loss, new_observed; kwargs...),
        update_observed(model.diff, new_observed; kwargs...)
        )
end

function update_observed(loss::SemLoss, new_observed; kwargs...)
    new_functions = Tuple(update_observed(lossfun, new_observed; kwargs...) for lossfun in loss.functions)
    return SemLoss(
        new_functions,
        loss.weights,
        loss.objective,
        loss.gradient,
        loss.hessian)
end


text:
mutable struct SimpleTask{T,T1} <: AbstractTask
    scheduler::T
    visible::Bool
    realtime::Bool
    speed::Float64
    screen::Vector{VisualObject}
    data::T1
end

function SimpleTask(;
    scheduler,
    visible=false,
    realtime=false,
    speed=1.0,
    screen=Vector{VisualObject}(),
    )
    SimpleTask(scheduler, visible, realtime, speed, screen, DataFrame())
end

function start!(task::SimpleTask, model)
    isi = 2.0
    description = "Present Stimulus"
    register!(task, present_stimulus, after, isi, task, model;
        description)
end

function present_stimulus(task, model)
    vo = VisualObject(;text="hello")
    add_to_visicon!(model, vo; stuff=true)
    push!(task.screen, vo)
end

function press_key!(task::SimpleTask, actr, key)
    empty!(task.screen)
    push!(task.data, (key=key,))
end


text:
module ExtendedVariantModeExtendedVariantController

using ..Ahorn, Maple

@mapdef Entity "ExtendedVariantMode/ExtendedVariantController" ExtendedVariantController(x::Integer, y::Integer, variantChange::String="Gravity", enable::Bool=true, newValue::Integer=10, revertOnLeave::Bool=false)

const placements = Ahorn.PlacementDict(
	"Extended Variant Controller (Extended Variant Mode)" => Ahorn.EntityPlacement(
		ExtendedVariantController,
		"rectangle"
	)
)

Ahorn.editingOptions(entity::ExtendedVariantController) = Dict{String, Any}(
	"variantChange" => Ahorn.ExtendedVariantDictionary.Variants
)

function Ahorn.selection(entity::ExtendedVariantController)
	x, y = Ahorn.position(entity)

	return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ExtendedVariantController, room::Maple.Room) = Ahorn.drawImage(ctx, "ahorn/ExtendedVariantMode/extended_variant_controller", -12, -12)

end


text:
import Base: reshape, adjoint, transpose, *, sum

export matmul, linear, reshape, transpose, adjoint, broadcastto, sumto

"""
    MatMul <: Func
"""
@func mutable struct MatMul end

forward(f::MatMul, w, x) = w * x

function backward(f::MatMul, gy)
    w, x = f.args
    gw = matmul(gy, transpose(x))
    gx = matmul(transpose(w), gy)
    return gw, gx
end

matmul(w, x) = MatMul()(w, x)

Base.:*(A::Var, B::Var) = matmul(A, B)
Base.:*(A::Var, B) = matmul(A, B)
Base.:*(A, B::Var) = matmul(A, B)

"""
    Linear_F <: Func
"""
@func mutable struct Linear_F end

forward(f::Linear_F, w, x, b=nothing) = b isa Nothing ? w * x : w * x .+ b

function backward(f::Linear_F, gy)
    w, x, b = length(f.args) == 3 ? f.args : (f.args..., nothing)
    gb = b isa Nothing ? nothing : sumto(gy, size(b))
    gw = matmul(gy, transpose(x))
    gx = matmul(transpose(w), gy)
    return gw, gx, gb
end

linear(w, x, b=nothing) = b isa Nothing ? Linear_F()(w, x) : Linear_F()(w, x, b)

"""
    Reshape <: Func
"""
@func mutable struct Reshape
    shape::Tuple
    Reshape(shape) = new(shape)
end

forward(f::Reshape, x) = reshape(x, f.shape)

backward(f::Reshape, gy) = reshape(gy, size(f.args[1]))

function reshape(x, shape...)
    if length(shape) == 1 && shape[1] isa Union{Tuple,AbstractArray}
        shape = shape[1]
    end
    size(x) == shape && return asvar(x)
    return Reshape(shape)(x)
end

adjoint(x::Var) = transpose(x)

"""
    transpose <: Func
"""
@func mutable struct Transpose end

forward(f::Transpose, x) = transpose(x)

backward(f::Transpose, gy) = transpose(gy)

transpose(x::Var) = Transpose()(x)

"""
    Sum <: Func
"""
@func mutable struct Sum
    dims::Union{Int,Tuple,Nothing}
    Sum(dims) = new(dims)
end

function forward(f::Sum, x)
    if f.dims isa Nothing
        return sum(x)
    else
        return sum(x, dims=f.dims)
    end
end

backward(f::Sum, gy) = broadcastto(gy, size(args[1]))

sum(x::Var; dims=nothing) = Sum(dims)(x)

"""
    BroadcastTo <: Func
"""
@func mutable struct BroadcastTo
    shape::Tuple
    BroadcastTo(shape) = new(shape)
end

forward(f::BroadcastTo, x) = x .* ones(f.shape)

backward(f::BroadcastTo, gy) = sumto(gy, size(f.args[1]))

broadcastto(x::Var, shape) = size(x) == shape ? asvar(x) : BroadcastTo(shape)(x)

"""
    SumTo <: Func
"""
@func mutable struct SumTo
    shape::Tuple
    SumTo(shape) = new(shape)
end

forward(f::SumTo, x) = _sumto(x, f.shape)

backward(f::SumTo, gy) = broadcastto(gy, size(f.args[1]))

sumto(x::Var, shape) = size(x) == shape ? asvar(x) : SumTo(shape)(x)


text:
gcs_prefix = "gs://gresearch/deep-ope/d4rl"
folder_prefix = "deep-ope-d4rl"
policies = D4RL_POLICIES

function deep_ope_d4rl_init()
    for policy in policies
        gcs_policy_folder = policy["policy_path"]
        local_policy_folder = chop(split(gcs_policy_folder, "/")[end], head=0, tail=4)
        register(
            DataDep(
                "$(folder_prefix)-$(local_policy_folder)",
                """
                Policy: Benchmarks for Deep Off-Policy Evaluation
                Credits: https://openreview.net/forum?id=kWSeGEeHvF8
                Url: https://github.com/google-research/deep_ope
                Authors: Justin Fu, Mohammad Norouzi, Ofir Nachum, George Tucker, ziyu wang, Alexander Novikov, Mengjiao Yang, Michael R Zhang,
                Yutian Chen, Aviral Kumar, Cosmin Paduraru, Sergey Levine, Thomas Paine
                Year: 2021

                Deep OPE contains:
                Policies for the tasks in the D4RL, DeepMind Locomotion and Control Suite datasets.
                Policies trained with the following algorithms (D4PG, ABM, CRR, SAC, DAPG and BC) and snapshots along the training trajectory. This facilitates
                benchmarking offline model selection.

                D4RL:
                A subset of the tasks within the D4RL (Fu et. al. 2020) for offline reinforcement learning is included. These tasks include maze navigation
                with different robot morphologies, hand manipulation tasks (Rajeswaran et. al. 2017), and tasks from the OpenAI Gym bechmark (Brockman et. al. 2016).
                Each task includes a variety of datasets in order to study the interaction between dataset distributions and policies. For further information on
                what datasets are available, please refer to D4RL: Datasets for Deep Data-Driven Reinforcement Learning.
                """,
                "$(gcs_prefix)/$(gcs_policy_folder)";
                fetch_method=fetch_gc_file
            )
        )
    end
end


text:
#

"""
    const CODE128

`CODE128` is a DataFrame with columns `value`, `code128a`, `code128b`, `code128c`, and
`pattern`, containing the values and binary patterns for all three modes (or subtypes)
`code128a`, `code128b` and `code128c` encodings.
"""
const CODE128 = DataFrame(
    value = 0:106,
    code128a = [
        string.(Char(32):Char(95)) # ' ' to '_'
        # string.(Char(0):Char(31)) # NUL to US
        [
            "NUL" # Null - Char(0) = \x00
            "SOH" # Start of Header
            "STX" # Start of Text
            "ETX" # End of Text
            "EOT" # End of Transmission
            "ENQ" # Enquiry
            "ACK" # Acknowledge
            "BEL" # Bell
            "BS" # Backspace
            "HT" # Horizontal Tab
            "LF" # Line Feed
            "VT" # Vertical Tab
            "FF" # Form Feed
            "CR" # Carriage Return
            "SO" # Shift Out
            "SI" # Shift In
            "DLE" # Data Link Escape
            "DC1" # Device Control 1
            "DC2" # Device Control 2
            "DC3" # Device Control 3
            "DC4" # Device Control 4
            "NAK" # Negative Acknowledge
            "SYN" # Synchronize
            "ETB" # End of Transmission Block
            "CAN" # Cancel
            "EM" # End of Medium
            "SUB" # Substitute
            "ESC" # Escape
            "FS" # File Separator
            "GS" # Group Separator
            "RS" #	Record Separator
            "US" #	Unit Separator - Char(31)
        ] # Char(0):Char(31)
        [
            "FNC 3"
            "FNC 2"
            "SHIFT B"
            "CODE C"
            "CODE B"
            "FNC 4"
            "FNC 1"
            "START A"
            "START B"
            "START C"
            "STOP"
        ]
    ],
    code128b = [
        string.(Char(32):Char(126)) # SPACE to ~
        ["DEL", "FNC 3", "FNC 2", "SHIFT A", "CODE C", "FNC 4", "CODE A", "FNC 1"]
        ["START A", "START B", "START C", "STOP"]
    ],
    code128c = [
        string.(0:99, pad = 2)
        ["CODE B", "CODE A", "FNC 1", "START A", "START B", "START C", "STOP"]
    ],
    pattern = [
        "11011001100"
        "11001101100"
        "11001100110"
        "10010011000"
        "10010001100"
        "10001001100"
        "10011001000"
        "10011000100"
        "10001100100"
        "11001001000"
        "11001000100"
        "11000100100"
        "10110011100"
        "10011011100"
        "10011001110"
        "10111001100"
        "10011101100"
        "10011100110"
        "11001110010"
        "11001011100"
        "11001001110"
        "11011100100"
        "11001110100"
        "11101101110"
        "11101001100"
        "11100101100"
        "11100100110"
        "11101100100"
        "11100110100"
        "11100110010"
        "11011011000"
        "11011000110"
        "11000110110"
        "10100011000"
        "10001011000"
        "10001000110"
        "10110001000"
        "10001101000"
        "10001100010"
        "11010001000"
        "11000101000"
        "11000100010"
        "10110111000"
        "10110001110"
        "10001101110"
        "10111011000"
        "10111000110"
        "10001110110"
        "11101110110"
        "11010001110"
        "11000101110"
        "11011101000"
        "11011100010"
        "11011101110"
        "11101011000"
        "11101000110"
        "11100010110"
        "11101101000"
        "11101100010"
        "11100011010"
        "11101111010"
        "11001000010"
        "11110001010"
        "10100110000"
        "10100001100"
        "10010110000"
        "10010000110"
        "10000101100"
        "10000100110"
        "10110010000"
        "10110000100"
        "10011010000"
        "10011000010"
        "10000110100"
        "10000110010"
        "11000010010"
        "11001010000"
        "11110111010"
        "11000010100"
        "10001111010"
        "10100111100"
        "10010111100"
        "10010011110"
        "10111100100"
        "10011110100"
        "10011110010"
        "11110100100"
        "11110010100"
        "11110010010"
        "11011011110"
        "11011110110"
        "11110110110"
        "10101111000"
        "10100011110"
        "10001011110"
        "10111101000"
        "10111100010"
        "11110101000"
        "11110100010"
        "10111011110"
        "10111101110"
        "11101011110"
        "11110101110"
        "11010000100"
        "11010010000"
        "11010011100"
        "11000111010"
    ],
)


text:
using Inventory
using LightGraphs
using Test

include("items.jl")
include("relationships.jl")

house = House("166 rue Jean Bart", "Sanary-sur-Mer", "France")
bedroom = Room(house)
gd = GraphDatabase()
add_item!(gd, house, tags=["house"], description="Parents' house")
add_item!(gd, bedroom, description="Cedric's bedroom in his parents' house")
add_item!(gd, "Plug extension x3", tags=["plugs", "electric", "accessories"])
add_item!(gd, "Side power strip - Power strip with plugs on the sides (left, right, top/bottom)", tags=["plugs", "electric", "accessories"])

add_relationship!(gd, Contains(house, bedroom))
gd

@testset "Items" begin
    @test find_vert(gd, house) == 1
    @test find_vert(gd, bedroom) == 2
    @test props(gd, 1)[:tags] == ["house"]
    @test props(gd, 1)[:quantity] == 1
    @test props(gd, 1)[:description] == "Parents' house."
    @test isempty(props(gd, 2)[:tags])
end

@testset "Relationships" begin
    @test find_edges(gd, Contains) == [Edge(1, 2)]
end

@testset "Queries" begin
    @test find_vertices(gd, Dict(:tags => ["house"])) == [1]
    @test find_vertices(gd, [:description]) == [1, 2, 3, 4]
end


text:
using Pkg, DrWatson

begin
	using ParetoSmooth
	using ParetoSmoothedImportanceSampling
	using StanSample
	using StatisticalRethinking
	using StatisticalRethinkingPlots
end

begin
	df = CSV.read(sr_datadir("rugged.csv"), DataFrame)
	dropmissing!(df, :rgdppc_2000)
	dropmissing!(df, :rugged)
	df.log_gdp = log.(df[:, :rgdppc_2000])
	df.log_gdp_s = df.log_gdp / mean(df.log_gdp)
	df.rugged_s = df.rugged / maximum(df.rugged)
	df.cid = [df.cont_africa[i] == 1 ? 1 : 2 for i in 1:size(df, 1)]

	data = (N = size(df, 1), K = length(unique(df.cid)),
		G = df.log_gdp_s, R = df.rugged_s, cid=df.cid)

	PRECIS(df[:, [:rgdppc_2000, :log_gdp, :log_gdp_s, :rugged, :rugged_s, :cid]])
end

stan8_1 = "
data {
	int N;
	vector[N] G;
	vector[N] R;
}

parameters {
	real a;
	real b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	mu = a + b * (R - 0.125);
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

# ╔═╡ 8e815d5c-6a66-11eb-21dc-99734c31f4e1
begin
	m8_1s = SampleModel("m8.1s", stan8_1)
	rc8_1s = stan_sample(m8_1s; data)
end

stan8_2 = "
data {
	int N;
	int K;
	vector[N] G;
	vector[N] R;
	int cid[N];
}

parameters {
	vector[K] a;
	real b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	for (i in 1:N)
		mu[i] = a[cid[i]] + b * (R[i] - 0.215);
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

begin
	m8_2s = SampleModel("m8.2s", stan8_2)
	rc8_2s = stan_sample(m8_2s; data)
end

stan8_3 = "
data {
	int N;
	int K;
	vector[N] G;
	vector[N] R;
	int cid[N];
}

parameters {
	vector[K] a;
	vector[K] b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	for (i in 1:N)
		mu[i] = a[cid[i]] + b[cid[i]] * (R[i] - 0.215);
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

begin
	m8_3s = SampleModel("m8.3s", stan8_3)
	rc8_3s = stan_sample(m8_3s; data)
end

begin
	loo8_1s, loos8_1s, pk8_1s = psisloo(m8_1s)
	pk_qualify(pk8_1s)
end

begin
	loo8_2s, loos8_2s, pk8_2s = psisloo(m8_2s)
	pk_qualify(pk8_2s)
end

begin
	loo8_3s, loos8_3s, pk8_3s = psisloo(m8_3s)
	pk_qualify(pk8_3s)
end

include(joinpath(@__DIR__, "PsisLooCompare.jl"))

models = [m8_1s, m8_2s, m8_3s]

psis_loo_comparison = psis_loo_compare(models)

psis_loo_comparison.df |> display
println()

compare([m8_1s, m8_2s, m8_3s], :psis) |> display


text:
using Documenter, juliaAllotter

makedocs(;
    modules=[juliaAllotter],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/saurabhkm/juliaAllotter.jl/blob/{commit}{path}#L{line}",
    sitename="juliaAllotter.jl",
    authors="Saurabh Kumar, Indian Institute of Technology Bombay",
    assets=String[],
)

deploydocs(;
    repo="github.com/saurabhkm/juliaAllotter.jl",
)

max_stars_repo_name: saurabhkm/juliaAllotter

max_issues_repo_name: saurabhkm/juliaAllotter

max_forks_repo_name: saurabhkm/juliaAllotter


text:
using Documenter, FreeRay

makedocs(sitename="FreeRay")

deploydocs(
    repo    = "github.com/ducphucnguyen/FreeRay.jl.git",
    target  = "build",
    deps    = nothing,
    make    = nothing
)

#serve(dir="docs/build")


text:
using Documenter
using ManpowerPlanning

makedocs(
    sitename = "ManpowerPlanning",
    authors = "Johan Van Kerckhoven"
)


text:
using TimerOutputs, Printf
using Oceananigans

include("benchmark_utils.jl")

#####
##### Benchmark setup and parameters
#####

const timer = TimerOutput()

Ni = 2  # Number of iterations before benchmarking starts.
Nt = 5  # Number of iterations to use for benchmarking time stepping.

# Model resolutions to benchmarks. Focusing on 3D models for GPU benchmarking.
            Ns = [(32, 32, 32), (64, 64, 64), (128, 128, 128), (256, 256, 256)]
   float_types = [Float32, Float64]  # Float types to benchmark.
         archs = [CPU()]             # Architectures to benchmark on.
@hascuda archs = [CPU(), GPU()]      # Benchmark GPU on systems with CUDA-enabled GPUs.

#####
##### Run benchmarks
#####

for arch in archs, FT in float_types, N in Ns
    Nx, Ny, Nz = N
    Lx, Ly, Lz = 250e3, 250e3, 1000

    model = ChannelModel(architecture = arch, float_type = FT,
			 grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Lx, Lz)),
			 closure = ConstantIsotropicDiffusivity(ν=1e-2, κ=1e-2))

    # Eddying channel model setup.
    Ty = 4e-5  # Meridional temperature gradient [K/m].
    Tz = 2e-3  # Vertical temperature gradient [K/m].

    # Initial temperature field [°C].
    T₀(x, y, z) = 10 + Ty*y + Tz*z + 1e-4*randn()
    set!(model; T=T₀)

    time_step!(model, Ni, 1)  # First 1~2 iterations are usually slower.

    bname =  benchmark_name(N, "", arch, FT)
    @printf("Running eddying channel benchmark: %s...\n", bname)
    for i in 1:Nt
        @timeit timer bname time_step!(model, 1, 1)
    end
end

#####
##### Print benchmark results
#####

println()
print_benchmark_info()

print_timer(timer, title="Eddying channel benchmarks", sortby=:name)

println("\n\nCPU Float64 -> Float32 speedup:")
for N in Ns
    bn32 = benchmark_name(N, "", CPU(), Float32)
    bn64 = benchmark_name(N, "", CPU(), Float64)
    t32  = TimerOutputs.time(timer[bn32])
    t64  = TimerOutputs.time(timer[bn64])
    @printf("%s: %.3f\n", benchmark_name(N), t64/t32)
end

@hascuda begin
    println("\nGPU Float64 -> Float32 speedup:")
    for N in Ns
        bn32 = benchmark_name(N, "", GPU(), Float32)
        bn64 = benchmark_name(N, "", GPU(), Float64)
        t32  = TimerOutputs.time(timer[bn32])
        t64  = TimerOutputs.time(timer[bn64])
        @printf("%s: %.3f\n", benchmark_name(N), t64/t32)
    end

    println("\nCPU -> GPU speedup:")
    for N in Ns, ft in float_types
        bn_cpu = benchmark_name(N, "", CPU(), ft)
        bn_gpu = benchmark_name(N, "", GPU(), ft)
        t_cpu  = TimerOutputs.time(timer[bn_cpu])
        t_gpu  = TimerOutputs.time(timer[bn_gpu])
        @printf("%s: %.3f\n", benchmark_name(N, ft), t_cpu/t_gpu)
    end
end


text:
mutable struct Point{T}
  x::T
  y::T
end

Point{Int64} #> Point{Int64} (constructor with 1 method)
Point{Float64} #> Point{Float64} (constructor with 1 method)
Point{String} #> Point{String} (constructor with 1 method)

Point{Int64} <: Point #> true
Point{String} <: Point #> true
Point{Float64} <: Point{Real} #> false

p = Point{Int64}(2, 5) #> Point{Int64}(2,5)
p = Point(2, 5) #> Point{Int64}(2,5)
p = Point("London", "Great-Britain") #> Point{String}("London", "Great-Britain")

mutable struct PointP{T <: Real}
  x::T
  y::T
end

# p = PointP("London", "Great-Britain") #> ERROR: MethodError: no method matching PointP(::String, ::String)

add(x::T, y::T) where T = x + y
add(2, 3) #> 5
# add(2, 3.0) #> ERROR: MethodError: `add` has no method matching add(::Int64, ::Float64)
add(x::T, y::T) where T <: Number = x + y

function vecfloat(x::Vector{T}) where T <: AbstractFloat
# code
end


text:
#=
Created on Thurday 02 January 2019
Last update: Tuesday 28 January 2020

@author: Michiel Stock
michielfmstock@gmail.com

General stuff.
=#

# COLORS
# ------

myblue = "#304da5"
mygreen = "#2a9d8f"
myyellow = "#e9c46a"
myorange = "#f4a261"
myred = "#e76f51"
myblack = "#50514F"

mycolors = [myblue, myred, mygreen, myorange, myyellow]

# TRACING
# -------

abstract type Tracker end

struct NoTrack <: Tracker
end

notrack = NoTrack()

struct PathTrack{T<:Any} <: Tracker
    xsteps::Array{T,1}
    PathTrack(x::T) where {T} = new{T}(copy(T[x]))
end

trace(::NoTrack, x) = nothing
trace(tracker::PathTrack, x) = push!(tracker.xsteps, copy(x))

nsteps(tracker::PathTrack) = length(tracker.xsteps) - 1

function getx(tracker::PathTrack)
    n = length(tracker.xsteps)
    return [x[1] for x in tracker.xsteps], [x[2] for x in tracker.xsteps]
end

# PLOTS
# -----

path(tracker::PathTrack; kwargs...) = plot(getx(tracker)...; color=mygreen, lw=2, kwargs...)
path!(tracker::PathTrack; kwargs...) = plot!(getx(tracker)...; color=myorange, lw=2, kwargs...)

plotobj(f::Function, tracker::PathTrack; kwargs...) = plot(0:nsteps(tracker), f.(tracker.xsteps),
                                                 lw=2, color=myorange, xlabel="iteration"; kwargs...)
plotobj!(f::Function, tracker::PathTrack; kwargs...) = plot!(0:nsteps(tracker), f.(tracker.xsteps);
                                                lw=2,color=mygreen, kwargs...)

using Colors
colorscatter(colors; kwargs...) = scatter(red.(colors), green.(colors), blue.(colors),
                        xlabel="red", ylabel="green", zlabel="blue", color=colors, label="")

# FUNCTIONS
# --------
"""
Compute Euclidean distance between two vectors.
"""
dist(x::AbstractVector, y::AbstractVector) = sqrt(sum((x .- y).^2))
"""
Compute Euclidean distance matrix between two matrices.
"""
dist(X::AbstractMatrix, Y::AbstractMatrix) = [dist(X[i,:], Y[j,:]) for i in 1:size(X,1), j in 1:size(Y,1)]
"""
Compute Euclidean distance matrix.
"""
dist(X::AbstractMatrix) = dist(X::AbstractMatrix, X::AbstractMatrix)

# GRAPHS
# -----

# special types for graphs
EdgeList{T} = Array{Tuple{T,T},1}
WeightedEdgeList{R,T} = Array{Tuple{R,T,T},1}
Vertices{T} = Array{T,1}
AdjList{R,T} = Dict{T,Array{Tuple{R,T},1}}

"""
Turns a list of weighted edges in an adjacency matrix (implemented as a Dict).
If the keyword `double` is set to `true`, every edge is added twice: `(w, u, v)`
and `(w, v, u)`. This is the default behaviour.
"""
function edges2adjlist(edges::WeightedEdgeList{R,T}; double=true) where {R<:Real,T}
    adjlist = AdjList{R,T}()
    for (w, i, j) in edges
        if !haskey(adjlist, i); adjlist[i] = [] end
        if !haskey(adjlist, j); adjlist[j] = [] end
        push!(adjlist[i], (w, j))
        double && push!(adjlist[j], (w, i))
    end
    return adjlist
end

"""
Turns an adjacency list (implemented as a Dict) into an edge list.
"""
adjlist2edges(adjlist::AdjList{R,T}) where {R<:Real, T} =
            [(w, v, n) for (v, neighbors) in adjlist for (w, n) in neighbors]

"""
Returns the number of vertices in a graph.
"""
nvertices(adjlist::AdjList) = length(adjlist)

"""
Returns the vertices of a graph.
"""
vertices(adjlist::AdjList) = Set(keys(adjlist))

"""
Returns the vertices of a graph.
"""
function vertices(edgelist::WeightedEdgeList{R,T}) where {R<:Real,T}
    vertices = Set{T}()
    for (w, u, v) in edgelist
        push!(vertices, u)
        push!(vertices, v)
    end
    return vertices
end

"""
Returns the number of vertices in a graph.
"""
nvertices(edgelist::WeightedEdgeList) = length(vertices(edgelist))

function isconnected(adjlist::AdjList{R,T}) where {R<:Real, T}
    visited = Set{T}()
    to_explore = [first(keys(adjlist))]
    while length(to_explore) > 0
        u = pop!(to_explore)
        push!(visited, u)
        for (w, n) in adjlist[u]
            n ∉ visited && push!(to_explore, n)
        end
    end
    return length(visited) == nvertices(adjlist)
end

isconnected(edges::WeightedEdgeList) = isconnected(edges2adjlist(edges))

#= TODO: write routines for

- [ ] number of edges
- [x] number of vertices
- [x] is connected
- [ ] is tree
=#

max_stars_repo_name: mferrodo/STMO

max_issues_repo_name: mferrodo/STMO

max_forks_repo_name: mferrodo/STMO


text:
export Pinger, stop,
       Ponger, pong_received, attach, ping_sent

mutable struct Pinger <: AbstractPinger
    timer::Union{Timer, Nothing}
    interval::Float64

    function Pinger(interval::Float64)
        new(nothing, interval)
    end
end

function attach(pinger::Pinger, logic::AbstractClientProtocol)
    send_ping = x -> handle(logic, ClientPingRequest())
    pinger.timer = Timer(send_ping, pinger.interval, interval=pinger.interval)
end

function stop(p::Pinger)
    if p.timer != nothing
        close(p.timer)
        p.timer = nothing
    end
end

mutable struct Ponger <: AbstractPonger
    timeout::Float64
    pong_missed::Function
    pongs_received::UInt64
    misses::Int
    current_misses::Int

    Ponger(timeout::Float64; misses::Int=1) = new(timeout, x -> nothing, 0, misses, 0)
end

attach(ponger::Ponger, logic::AbstractClientProtocol) =
    ponger.pong_missed = () -> handle(logic, PongMissed())

function pong_received(ponger::Ponger)
    ponger.pongs_received += 1
    ponger.current_misses = 0
end

function ping_sent(ponger::Ponger)
    pongs_received_at_send = ponger.pongs_received
    fun = x -> begin
        ponger.current_misses += 1
        if ponger.pongs_received == pongs_received_at_send && ponger.current_misses >= ponger.misses
            ponger.pong_missed()
            ponger.current_misses = 0
        end
    end
    Timer(fun, ponger.timeout)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AvailabilitySetUpdate <: SwaggerModel
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags

    function AvailabilitySetUpdate(;tags=nothing)
        o = new()
        validate_property(AvailabilitySetUpdate, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        o
    end
end # type AvailabilitySetUpdate

const _property_map_AvailabilitySetUpdate = Dict{Symbol,Symbol}(Symbol("tags")=>Symbol("tags"))
const _property_types_AvailabilitySetUpdate = Dict{Symbol,String}(Symbol("tags")=>"Dict{String, String}")
Base.propertynames(::Type{ AvailabilitySetUpdate }) = collect(keys(_property_map_AvailabilitySetUpdate))
Swagger.property_type(::Type{ AvailabilitySetUpdate }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AvailabilitySetUpdate[name]))}
Swagger.field_name(::Type{ AvailabilitySetUpdate }, property_name::Symbol) =  _property_map_AvailabilitySetUpdate[property_name]

function check_required(o::AvailabilitySetUpdate)
    true
end

function validate_property(::Type{ AvailabilitySetUpdate }, name::Symbol, val)
end


text:
@rule BIFM(:in, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_zprev::ProdFinal{<:MultivariateNormalDistributionsFamily}, m_znext::MultivariateNormalDistributionsFamily, meta::BIFMMeta) = begin

    # fetch information of meta data
    A                   = getA(meta)
    B                   = getB(meta)
    μu                  = getμu(meta)
    Σu                  = getΣu(meta)
    ξztilde             = getξztilde(meta)
    Λztilde             = getΛztilde(meta)

    # fetch statistics from messages
    μ_zprev, Σ_zprev    = mean_cov(m_zprev)

    # calculate intermediate variables (with dual parameterization)
    # see Wadehn 2016 - On sparsity by NUV-EM ...
    ξtildex = Λztilde * (A * μ_zprev) - ξztilde
    Λtildex = Λztilde - Λztilde * A * Σ_zprev * A' * Λztilde
    tmp = B * Σu

    # calculate marginals of input WRONG
    μ_in = μu - (Σu * (B' * ξtildex))
    Σ_in = Σu - tmp' * Λtildex * tmp

    # Actual return type depends on meta object as well, so we explicitly cast the result here
    # Should be noop if type matches
    T = promote_type(eltype(m_out), eltype(m_zprev), eltype(m_znext))

    # return input marginal
    return ProdFinal(convert(MvNormalMeanCovariance{T}, MvNormalMeanCovariance(μ_in, Σ_in)))

end


text:
#=

  Balanced brackets in ConstraintSolver.jl

  This model all balanced brackets of size m*2.

  See https://en.wikipedia.org/wiki/Catalan_number

  The number of generated solutions for m:
   m        #
   ----------
    1       1
    2       2
    3       5
    4      14
    5      42
    6     132
    7     429
    8    1430
    9    4862
   10   16796
   11   58786
   12  208012
   13  742900

  Which - of course - is the Catalan numbers.

  http://oeis.org/search?q=1%2C2%2C5%2C14%2C42%2C132%2C429%2C1430%2C4862%2C16796%2C58786%2C208012&language=english&go=Search
  http://oeis.org/A000108

  This model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function balanced_brackets(m=4,print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            "lp_optimizer" => ipopt_optimizer,
                                        ))
    n = m*2
    t = [1,-1] # +1 for "[", -1 for "]"
    @variable(model, 1 <= x[1:n] <= 2, Int) # 1="[",  2="]"

    # Counter (cumulative sum)
    @variable(model, 0 <= c[1:n] <= n, Int)

    # Leading "["
    @constraint(model, x[1] == 1)
    @constraint(model, c[1] == 1)

    for i in 2:n
        tt = @variable(model, [1:1], CS.Integers([-1,1]))
        my_element(model,x[i], t, tt[1])
        @constraint(model, c[i] == c[i-1] + tt[1])
    end

    # Concluding "]"
    @constraint(model, x[n] == 2)
    @constraint(model, c[n] == 0)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                if print_solutions
                    x_val = convert.(Integer,JuMP.value.(x; result=sol))
                    c_val = convert.(Integer,JuMP.value.(c; result=sol))
                    println("x:$x_val")
                    println("c:$c_val")
                    println()
                end

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

begin
    @time local status, num_sols = balanced_brackets(4,true,true)
    println("num_sols: $num_sols")

    for m in 1:20
        println("m:$m")
        @time local status, num_sols = balanced_brackets(m,false,true)
        if status == MOI.TIME_LIMIT
            break
        end
    end
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module triangleInterpolator
    using Images
    export rasterizationBBOX

    function pointLine(x::Float64,
                       y::Float64,
                       line::Array{Float64},
                       linex::Float64,
                       liney::Float64
            )::Float64
        return line[2]*x - line[1]*y - line[2]*linex + line[1]*liney
    end

    function swap_points(ax::Float64,
                         ay::Float64,
                         bx::Float64,
                         by::Float64,
                         colorA::RGB{N0f8},
                         colorB::RGB{N0f8}
            )::Tuple{Float64, Float64, Float64, Float64, RGB{N0f8}, RGB{N0f8}}
        buf_x = ax
        buf_y = ay
        ax = bx
        ay = by
        bx = buf_x
        by = buf_y

        colorBuf = colorA
        colorA = colorB
        colorB = colorBuf

        return ax, ay, bx, by, colorA, colorB
    end
    function validate_entry(ax::Float64,
                            ay::Float64,
                            bx::Float64,
                            by::Float64,
                            cx::Float64,
                            cy::Float64,
                            colorA::RGB{N0f8},
                            colorB::RGB{N0f8},
                            colorC::RGB{N0f8}
            )::Tuple{Float64, Float64, Float64, Float64, Float64, Float64, RGB{N0f8}, RGB{N0f8}, RGB{N0f8}}
        ab = [bx-ax, by-ay]
        bc = [cx-bx, cy-by]
        ca = [ax-cx, ay-cy]

        if pointLine(cx, cy, ab, ax, ay) < 0
            ax, ay, cx, cy, colorA, colorC = swap_points(ax, ay, cx, cy, colorA, colorC)

        elseif pointLine(ax, ay, bc, bx, by) < 0
            ax, ay, bx, by, colorA, colorB = swap_points(ax, ay, bx, by, colorA, colorB)

        elseif pointLine(bx, by, ca, cx, cy) < 0
            bx, by, cx, cy, colorB, colorC = swap_points(bx, by, cx, cy, colorB, colorC)
        end

        ax, ay, bx, by, cx, cy, colorA, colorB, colorC
    end
    function interpolateColors(position::Array{Float64},
                               ax::Float64,
                               ay::Float64,
                               bx::Float64,
                               by::Float64,
                               cx::Float64,
                               cy::Float64,
                               colorA::RGB{N0f8},
                               colorB::RGB{N0f8},
                               colorC::RGB{N0f8}
            )::RGB{N0f8}
        ab = [bx-ax, by-ay]
        bc = [cx-bx, cy-by]
        ca = [ax-cx, ay-cy]

        i, j = position

        alpha = pointLine(j, i, bc, bx, by) / pointLine(ax, ay, bc, bx, by)
        beta = pointLine(j, i, ca, cx, cy) / pointLine(bx, by, ca, cx, cy)
        phi = pointLine(j, i, ab, ax, ay) / pointLine(cx, cy, ab, ax, ay)

        return alpha * colorA + beta * colorB + phi * colorC
    end
    function setupBBOX(A::Array{Float64},
                       B::Array{Float64},
                       C::Array{Float64},
                       maxX::Int64,
                       maxY::Int64
            )::NTuple{4, Float64}
        maxHeight = max(A[1], B[1], C[1])
        minHeight = min(A[1], B[1], C[1])
        if maxHeight > maxY
            maxHeight = maxY
        end
        if minHeight < 1
            minHeight = 1
        end
        maxWidth = max(A[2], B[2], C[2])
        minWidth = min(A[2], B[2], C[2])
        if maxWidth > maxX
            maxWidth = maxX
        end
        if minWidth < 1
                minWidth = 1
        end
        return maxHeight, minHeight, maxWidth, minWidth
    end
    function rasterizationBBOX(img::Array{RGB{N0f8}, 2},
                               A::Array{Float64},
                               B::Array{Float64},
                               C::Array{Float64},
                               colorA::RGB{N0f8},
                               colorB::RGB{N0f8},
                               colorC::RGB{N0f8}
            )

        maxHeight, minHeight, maxWidth, minWidth = setupBBOX(A, B, C, size(img)[2], size(img)[1])
        ax, ay, bx, by, cx, cy = A[2], A[1], B[2], B[1], C[2], C[1]
        ax, ay, bx, by, cx, cy, colorA, colorB, colorC = validate_entry(ax, ay, bx, by, cx, cy, colorA, colorB, colorC)

        ab = [bx-ax, by-ay]
        bc = [cx-bx, cy-by]
        ca = [ax-cx, ay-cy]

        for i=floor(Int, minHeight):ceil(Int, maxHeight)
            for j=floor(Int, minWidth):ceil(Int, maxWidth)
                floatj = Float64(j)
                floati = Float64(i)
                alpha = pointLine(floatj, floati, bc, bx, by)
                beta = pointLine(floatj, floati, ca, cx, cy)
                phi = pointLine(floatj, floati, ab, ax, ay)

                if beta >= 0 && alpha >= 0 && phi >= 0
                    img[i, j] = interpolateColors([floati, floatj], ax, ay, bx, by, cx, cy, colorA, colorB, colorC)
                end
            end
        end
    end
end


text:
using Games
using Test

include("test_pure_nash.jl")
include("test_repeated_game.jl")
include("test_normal_form_game.jl")
include("test_random.jl")
include("test_support_enumeration.jl")

include("generators/runtests.jl")