text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(Base, Symbol("#@threads")) && precompile(Tuple{getfield(Base.Threads, Symbol("#@threads")), LineNumberNode, Module, Int})
    isdefined(Base, Symbol("#@views")) && precompile(Tuple{getfield(Base, Symbol("#@views")), LineNumberNode, Module, Int})
    precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf), typeof(Base._views), Expr})
    precompile(Tuple{typeof(Base.Broadcast.copyto_nonleaf!), Array{Expr, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Base.OneTo{Int64}, Int64, Int64})
    precompile(Tuple{typeof(Base.Broadcast.copyto_nonleaf!), Array{Symbol, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Base.OneTo{Int64}, Int64, Int64})
    precompile(Tuple{typeof(Base.Broadcast.instantiate), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(Base._views), Tuple{Array{Any, 1}}}})
    precompile(Tuple{typeof(Base.Broadcast.restart_copyto_nonleaf!), Array{Any, 1}, Array{Expr, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Int64, Int64, Base.OneTo{Int64}, Int64, Int64})
    precompile(Tuple{typeof(Base.Broadcast.restart_copyto_nonleaf!), Array{Any, 1}, Array{Symbol, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Expr, Int64, Base.OneTo{Int64}, Int64, Int64})
    precompile(Tuple{typeof(Base.Broadcast.restart_copyto_nonleaf!), Array{Any, 1}, Array{Symbol, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Int64, Int64, Base.OneTo{Int64}, Int64, Int64})
    precompile(Tuple{typeof(Base.Docs.docm), LineNumberNode, Module, Int, Int})
    precompile(Tuple{typeof(Base.Docs.docstr), Int, Int})
    precompile(Tuple{typeof(Base.Docs.moduledoc), Int, Int, Int, Int, Expr})
    precompile(Tuple{typeof(Base.Docs.objectdoc), Int, Int, Int, Int, Int, Int})
    precompile(Tuple{typeof(Base.MainInclude.include), String})
    precompile(Tuple{typeof(Base.Threads._threadsfor), Expr, Expr, Symbol})
    precompile(Tuple{typeof(Base.Threads.resize_nthreads!), Array{Base.MPFR.BigFloat, 1}, Base.MPFR.BigFloat})
    precompile(Tuple{typeof(Base._views), Expr})
    precompile(Tuple{typeof(Base.hashindex), Tuple{Int64, Int32, UInt64}, Int64})
    precompile(Tuple{typeof(Base.hashindex), Tuple{Int64, Nothing, UInt64}, Int64})
    precompile(Tuple{typeof(Base.ht_keyindex), Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
    precompile(Tuple{typeof(Base.hvcat), Tuple{Int64, Int64, Int64, Int64, Int64, Int64, Int64, Int64, Int64, Int64}, Float64, Float64})
    precompile(Tuple{typeof(Base.isassigned), Core.SimpleVector, Int64})
    precompile(Tuple{typeof(Base.push!), Array{Function, 1}, typeof(identity)})
    precompile(Tuple{typeof(Base.rehash!), Base.Dict{Int64, Tuple{Function, Int64, Vararg{Int64, N} where N}}, Int64})
    precompile(Tuple{typeof(Base.rehash!), Base.Dict{String, Tuple{Any, Any, Int64}}, Int64})
    precompile(Tuple{typeof(Base.rehash!), Base.Dict{Tuple{Int64, Any, UInt64}, DataType}, Int64})
    precompile(Tuple{typeof(Base.require), Module, Symbol})
    precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Base.Dict{String, Any}, Pkg.BinaryPlatforms.FreeBSD})
    precompile(Tuple{typeof(Base.sqrt), Float64})
end


text:
module FSimBase

using Reexport
import DiffEqBase
using DiffEqBase: DECallback, DEIntegrator, ODEProblem, DiscreteProblem, CallbackSet, isinplace, init, solve
using DiffEqCallbacks: SavedValues, SavingCallback
@reexport using SimulationLogger
using Transducers: Map
using DataFrames: DataFrame
using NamedTupleTools: namedtuple
using DataFrames

export AbstractEnv, State, Params, Dynamics, Dynamics!
export AbstractController, Command
export apply_inputs
export Simulator
export solve, reinit!, log!, step!, step_until!

include("APIs/APIs.jl")

end

ext: jl


text:
module MoreBasisFunctions

    using GridArrays: ScatteredGrid
    using BasisFunctions

    export Bernstein, BernsteinIndex, BernsteinInterval
    export Lagrange, LagrangeIndex, LagrangeInterval
    export LagrangeGaussLegendre, LagrangeLobattoLegendre
    export LagrangeGaussChebyshev, LagrangeLobattoChebyshev

    export eval_element, eval_element_derivative, eval_element_antiderivative
    export nodes, nnodes, degree, support

    include("bases/generic/dictionary.jl")
    include("bases/generic/nodes.jl")
    include("bases/generic/vandermonde_matrix.jl")

    include("bases/poly/bernstein.jl")
    include("bases/poly/lagrange.jl")

end # module


text:
# This file is a part of SearchModels.jl

export AbstractSolutionSpace, search_models, SearchParams, rand, combine_select, combine, mutate
using Distributed, Random, StatsBase, Parameters

abstract type AbstractSolutionSpace end
## NOTE: all configuratons must define hash, and isequal

"""
    config_type(::T)

Config type identifier, it may or not be a type
"""
config_type(::T) where T = Symbol(Base.typename(T))

#function rand(space::AbstractSolutionSpace) end
#function combine(a, b) end

"""
    rand(space::AbstractSolutionSpace)

Creates a random configuration sampling the given space
"""

"""
    compatible_space(space::AbstractVector, c)

Selects the compatible space for configuration c
"""
function compatible_space(space::AbstractVector, c)
    for s in space
        if c isa eltype(s)
            return s
        end
    end
    error("uncompatible space for $(typeof(c)) in space list $(typeof(space))")
end

compatible_space(space::AbstractSolutionSpace, c) = space

function compatible_config(c, L::AbstractVector)
    t_ = config_type(c)

    for p in L
        x = p.first
        if config_type(x) == t_
            return x
        end
    end

    error("uncompatible configuration for $(typeof(c))")
end

"""
    combine(space::AbstractSolutionSpace, c1, c2)
    combine(c1, c2)

Combines two configurations into a single configuration. The three argument functions defaults to two argument function.

"""
function combine end

"""
    combine_select(a, L::AbstractVector)

Combines `a` configuration with some compatible configuation in the given list of pairs (config => score),
The `a` config is always at the end of `L` and also `L` is always shuffled.
If you are in doubt, use the higher level interface `combine(c1, c2)`.
"""
function combine_select(a::T, L::AbstractVector) where T
    # L is a vector of pairs config => score
    # L should be shuffled before combining
    combine(a, compatible_config(a, L))
end

function mutate end
"""
    mutate(space::AbstractSolutionSpace, c, iter)
    mutate(space::AbstractVector, c, iter)

Mutates configuration. If space is a list of spaces, then the proper space is determined.
"""
function mutate(space::AbstractVector, c, iter)
    mutate(compatible_space(space, c), c, iter)
end

"""
    evaluate_queue(errfun::Function, evalqueue, population, parallel, tmp)

Evaluates queue of using errfun; the resulting `config => performances` are stored in `config_and_perf`

"""
function evaluate_queue(errfun::Function, evalqueue, population, parallel, tmp)
    if parallel === :none
        if population === nothing
            population = [c => errfun(c) for c in evalqueue]
        else
            for c in evalqueue
                push!(population, c => errfun(c))
            end
        end

        empty!(evalqueue)
        return population
    end

    empty!(tmp)
    while length(evalqueue) > 0
        c = pop!(evalqueue)
        if parallel === :distributed
            perf = Distributed.@spawn errfun(c)
            push!(tmp, c => perf)
        elseif parallel === :threads
            perf = Threads.@spawn errfun(c)
            push!(tmp, c => perf)
        else
            error("Unknown parallel option")
        end
    end

    if population === nothing
        [k => fetch(v) for (k, v) in tmp]
    else
        for (k, v) in tmp
            push!(population, k => fetch(v))
        end

        population
    end
end

"""
    queue_config!(accept_config, conf, evalqueue, observed)

Pushes `conf` into the evaluation queue. It checks if was already observed and if it is accepted by the `accept_config` predicate.
"""
function queue_config!(accept_config, conf, evalqueue, observed)
    if !(conf in observed) && accept_config(conf)
        push!(evalqueue, conf)
        push!(observed, conf)
    end
end

"""
    SearchParams(;
        maxpopulation::Int = 32
        bsize::Int = 16
        mutbsize::Int = 8
        crossbsize::Int = 8
        maxiters::Int = 300
        tol::Float64 = 0.001
        verbose::Bool = true
    )

reates a mutable list of search parameters, it can be changed online via `inspect_population`.

- `maxpopulation`: the maximum number of configurations to be kept at each iteration (based on its minimum error)
- `bsize`: number of best items to be used for mutation and crossing procedures
- `mutbsize`: number of new configurations per iteration from a mutation procedure
- `crossbsize`: number of new configurations per iteration from a crossing procedure
- `maxiters`: maximum iterations of the search procedure
- `tol`: stop search tolerance to population errors on the best `maxpopulation` configurations; negative values will force to evaluate `maxiters`
- `verbose`: controls if the verbosity of the search iterations

"""
@with_kw mutable struct SearchParams
    maxpopulation::Int = 16
    bsize::Int = 8
    mutbsize::Int = 8
    crossbsize::Int = 8
    maxiters::Int = 300
    tol::Float64 = 0.001
    verbose::Bool = true
end

"""
    search_models(
        errfun::Function,  # error function, it can be used with `do conf ... end` block
        space::AbstractSolutionSpace,
        initialpopulation=32,
        params::SearchParams=SearchParams();
        geterr::Function=identity,
        accept_config::Function=config->true,
        inspect_population::Function=(space, params, population) -> nothing,
        parallel=:none, # :none, :threads, :distributed
    )

Explores the search space trying to minimize the given error function. The procedure consists on an iterative stochastic method
based on an evolutionary algorithm. It is starts with `initialpopulation` configurations, evaluate them mutate and cross them;
it selects at most `maxpopulation` configurations at any iteration (best ones).
- `errfun`: the function to minimize (receives the configuration as argument)
- `space`: Search space definition
- `params`: search parameters
- `geterr` a function to compute or fetch the cost (a single floating point) on the output of `errfun`
- `initialpopulation`: initial number of configurations to conform the population, or an initial list of configurations
- `accept_config`: an alternative way to deny some configuration to be evaluated (receives the configuration as argument)
- `inspect_population`: observes the population after evaluating a beam of solutions
- `parallel`: controls if the search is made with some kind of parallel strategy (only used for evaluation errors). Valid values are:
  - `:none`: there is no parallelization, the default value
  - `:threads`: evaluates error functions using threads
  - `:distributed`: evaluates error functions using a distributed environment (using the available workers)
"""
function search_models(
        errfun::Function,
        space::AbstractSolutionSpace,
        initialpopulation=32, # it can alos be a list of config seeds
        params::SearchParams=SearchParams();
        geterr::Function=identity,
        accept_config::Function=config -> true,
        inspect_population::Function=(space, params, population) -> nothing,
        parallel=:none # :none, :threads, :distributed
    )

    ConfigType = eltype(space)
    evalqueue = ConfigType[]
    observed = Set{ConfigType}() # solutions should override hash and isqual functions

    if initialpopulation isa Integer
        for i in 1:initialpopulation
            queue_config!(accept_config, rand(space), evalqueue, observed)
        end
    else
        for c in initialpopulation
            queue_config!(accept_config, c, evalqueue, observed)
        end
    end
    prev = 0.0
    iter = 0

    population = nothing
    tmp = Pair[]
    params.verbose && println(stderr, "SearchModels> search params iter=$iter, tol=$(params.tol), initialpopulation=$initialpopulation, maxpopulation=$(params.maxpopulation), bsize=$(params.bsize), mutbsize=$(params.mutbsize), crossbsize=$(params.crossbsize)")

    while iter <= params.maxiters
        iter += 1
        population = evaluate_queue(errfun, evalqueue, population, parallel, tmp)
        sort!(population, by=x->geterr(x.second))
        inspect_population(space, params, population)

        if params.maxpopulation < length(population)
            resize!(population, params.maxpopulation)
        end

        if iter >= params.maxiters
            params.verbose && println("SearchModels> reached maximum number of iterations $(params.maxiters)")
            return population
        end

        curr = geterr(population[end].second)
        if abs(curr - prev) <= params.tol
            params.verbose && println("SearchModels> stop by convergence error=$curr, tol=$(params.tol)")
            return population
        end

        prev = curr
        ## params.verbose && println(stderr, "SearchModels> *** generating extra indivuals bsize=$bsize, mutbsize=$mutbsize, crossbsize=$crossbsize")
        best_error, worst_error = population[1].second, population[end].second ## the top configurations will be shuffled
        L = @view population[1:min(params.bsize, length(population))]
        for i in 1:params.mutbsize
            conf = mutate(space, rand(L).first, iter)
            queue_config!(accept_config, conf, evalqueue, observed)
        end

        for i in 1:params.crossbsize
            shuffle!(L) # the way this procedure is designed is to support heterogeneous options
            i = rand(1:length(L))
            L[end], L[i] = L[end], L[i]
            conf = combine_select(L[end].first, L)
            queue_config!(accept_config, conf, evalqueue, observed)
        end

        params.verbose && println(stderr, "SearchModels iteration $iter> population: $(length(population)), bsize: $(params.bsize), queue: $(length(evalqueue)), observed: $(length(observed)), best-error: $(best_error) worst-error: $(worst_error)")
    end

    sort!(population, by=x->geterr(x.second))
end


text:
const _graph_funcs = Dict{Symbol,Any}(
    :spectral => spectral_graph,
    :sfdp => sfdp_graph,
    :circular => circular_graph,
    :shell => shell_graph,
    :spring => spring_graph,
    :stress => by_axis_local_stress_graph,
    :tree => tree_graph,
    :buchheim => buchheim_graph,
    :arcdiagram => arc_diagram,
    :chorddiagram => chord_diagram
)

const _graph_inputs = Dict{Symbol,Any}(
    :spectral => :adjmat,
    :sfdp => :adjmat,
    :circular => :adjmat,
    :shell => :adjmat,
    :stress => :adjmat,
    :spring => :adjmat,
    :tree => :sourcedestiny,
    :buchheim => :adjlist,
    :arcdiagram => :sourcedestiny,
    :chorddiagram => :sourcedestiny
)

function prepare_graph_inputs(method::Symbol, inputs...)
    input_type = get(_graph_inputs, method, :sourcedestiny)
    if input_type == :adjmat
        (get_adjacency_matrix(inputs...),)
    elseif input_type == :sourcedestiny
        get_source_destiny_weight(inputs...)
    elseif input_type == :adjlist
        (get_adjacency_list(inputs...),)
    end
end

# -----------------------------------------------------

function get_source_destiny_weight(mat::AbstractArray{T,2}) where T
    nrow, ncol = size(mat)    # rows are sources and columns are destinies

    nosymmetric = !issymmetric(mat) # plots only triu for symmetric matrices
    nosparse = !issparse(mat) # doesn't plot zeros from a sparse matrix

    L = length(mat)

    source  = Array{Int}(undef, L)
    destiny = Array{Int}(undef, L)
    weights  = Array{T}(undef, L)

    idx = 1
    for i in 1:nrow, j in 1:ncol
        value = mat[i, j]
        if !isnan(value) && ( nosparse || value != zero(T) ) # TODO: deal with Nullable

            if i < j
                source[idx]  = i
                destiny[idx] = j
                weights[idx]  = value
                idx += 1
            elseif nosymmetric && (i > j)
                source[idx]  = i
                destiny[idx] = j
                weights[idx]  = value
                idx += 1
            end

        end
    end

    resize!(source, idx-1), resize!(destiny, idx-1), resize!(weights, idx-1)
end

function get_source_destiny_weight(source::AbstractVector, destiny::AbstractVector)
    if length(source) != length(destiny)
        throw(ArgumentError("Source and destiny must have the same length."))
    end
    source, destiny, Float64[ 1.0 for i in source ]
end

function get_source_destiny_weight(source::AbstractVector, destiny::AbstractVector, weights::AbstractVector)
    if !(length(source) == length(destiny) == length(weights))
        throw(ArgumentError("Source, destiny and weights must have the same length."))
    end
    source, destiny, weights
end

function get_source_destiny_weight(adjlist::AbstractVector{V}) where V<:AbstractVector{Int}
    source = Int[]
    destiny = Int[]
    for (i,l) in enumerate(adjlist)
        for j in l
            push!(source, i)
            push!(destiny, j)
        end
    end
    get_source_destiny_weight(source, destiny)
end

# -----------------------------------------------------

function get_adjacency_matrix(mat::AbstractMatrix)
    mat
end

function get_adjacency_matrix(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector)
    n = max(maximum(source), maximum(destiny))
    Matrix(sparse(source, destiny, weights, n, n))
end

function get_adjacency_matrix(adjlist::AbstractVector{V}) where V<:AbstractVector{Int}
    s,d,w = get_source_destiny_weight(adjlist)
    get_adjacency_matrix(s, d, w)
end

# -----------------------------------------------------

function get_adjacency_list(mat::AbstractMatrix)
    get_adjacency_list(get_source_destiny_weight(mat))
end

function get_adjacency_list(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector)
    n = max(maximum(source), maximum(destiny))
    adjlist = [Int[] for i=1:n]
    for (s,d) in zip(source,destiny)
        push!(adjlist[s], d)
    end
    adjlist
end

function get_adjacency_list(adjlist::AbstractVector{V}) where V<:AbstractVector{Int}
    adjlist
end

# -----------------------------------------------------

function make_symmetric(A::AbstractMatrix)
    A = copy(A)
    for i=1:size(A,1), j=(i+1):size(A,2)
        A[i,j] = A[j,i] = A[i,j]+A[j,i]
    end
    A
end

function compute_laplacian(adjmat::AbstractMatrix, node_weights::AbstractVector)
    n, m = size(adjmat)
    # @show size(adjmat), size(node_weights)
    @assert n == m == length(node_weights)

    # scale the edge values by the product of node_weights, so that "heavier" nodes also form
    # stronger connections
    adjmat = adjmat .* sqrt(node_weights * node_weights')

    # D is a diagonal matrix with the degrees (total weights for that node) on the diagonal
    deg = vec(sum(adjmat; dims=1)) - diag(adjmat)
    D = diagm(0 => deg)

    # Laplacian (L = D - adjmat)
    L = eltype(adjmat)[i == j ? deg[i] : -adjmat[i,j] for i=1:n,j=1:n]

    L, D
end

import LightGraphs

# TODO: so much wasteful conversion... do better
function estimate_distance(adjmat::AbstractMatrix)
    source, destiny, weights = get_source_destiny_weight(sparse(adjmat))

    g = LightGraphs.Graph(adjmat)
    dists = convert(Matrix{Float64}, hcat(map(i->LightGraphs.dijkstra_shortest_paths(g, i).dists, LightGraphs.vertices(g))...))
    tot = 0.0; cnt = 0
    for (i,d) in enumerate(dists)
        if d < 1e10
            tot += d
            cnt += 1
        end
    end
    avg = cnt > 0 ? tot / cnt : 1.0
    for (i,d) in enumerate(dists)
        if d > 1e10
            dists[i] = 3avg
        end
    end
    dists
end

function get_source_destiny_weight(g::LightGraphs.AbstractGraph)
    source = Vector{Int}()
    destiny =  Vector{Int}()
    sizehint!(source, LightGraphs.nv(g))
    sizehint!(destiny, LightGraphs.nv(g))
    for e in LightGraphs.edges(g)
      push!(source, LightGraphs.src(e))
      push!(destiny, LightGraphs.dst(e))
    end
    get_source_destiny_weight(source, destiny)
end

function get_adjacency_matrix(g::LightGraphs.AbstractGraph)
    adjacency_matrix(g)
end

function get_adjacency_matrix(source::AbstractVector{Int}, destiny::AbstractVector{Int})
    get_adjacency_matrix(source, destiny, ones(length(source)))
end

function get_adjacency_list(g::LightGraphs.AbstractGraph)
    g.fadjlist
end

# -----------------------------------------------------

# a graphplot takes in either an (N x N) adjacency matrix
#   note: you may want to pass node weights to markersize or marker_z
# A graph has N nodes where adj_mat[i,j] is the strength of edge i --> j.  (adj_mat[i,j]==0 implies no edge)

# NOTE: this is for undirected graphs... adjmat should be symmetric and non-negative

@userplot GraphPlot

@recipe function f(g::GraphPlot;
                   dim = 2,
                   free_dims = nothing,
                   T = Float64,
                   curves = true,
                   curvature_scalar = 0.05,
                   root = :top,
                   node_weights = nothing,
                   names = [],
                   fontsize = 7,
                   nodeshape = :hexagon,
                   nodesize = 0.1,
                   nodecolor = 1,
                   x = nothing,
                   y = nothing,
                   z = nothing,
                   method = :stress,
                   func = get(_graph_funcs, method, by_axis_local_stress_graph),
                   shorten = 0.0,
                   axis_buffer = 0.2,
                   layout_kw = Dict{Symbol,Any}(),
                   edgewidth = (s,d,w)->1,
                   edgelabel = nothing,
                   edgelabel_offset = 0.0,
                   self_edge_size = 0.1,
                   edge_label_box = true,
                  )
    @assert dim in (2, 3)
    _3d = dim == 3
    adj_mat = get_adjacency_matrix(g.args...)
    isdirected = (g.args[1] isa DiGraph || !issymmetric(adj_mat)) && !in(method, (:tree, :buchheim))

    source, destiny, weights = get_source_destiny_weight(g.args...)
    if !(eltype(source) <: Integer)
        names = unique(sort(vcat(source,destiny)))
        source = Int[findfirst(names, si) for si in source]
        destiny = Int[findfirst(names, di) for di in destiny]
    end
    n = max(maximum(source), maximum(destiny))

    if isnothing(node_weights)
        node_weights = ones(n)
    end
    @assert length(node_weights) == n

    xyz = _3d ? (x,y,z) : (x,y)
    numnothing = count(isnothing, xyz)

    # do we want to compute coordinates?
    if numnothing > 0
        if isnothing(free_dims)
            # compute free_dims
            free_dims = findall(isnothing, xyz)
        end
        x, y, z = func(
            prepare_graph_inputs(method, source, destiny, weights)...;
            node_weights = node_weights,
            dim = dim,
            free_dims = free_dims,
            root = root,
            layout_kw...
        )
    end

    # reorient the points after root
    if root in (:left,:right)
        x,y = y,-x
    end
    if root == :left
        x,y = -x, y
    end
    if root == :bottom
        x,y = x, -y
    end

    if method == :arcdiagram
        xl, yl = arcdiagram_limits(x, source, destiny)
        xlims --> xl
        ylims --> yl
        ratio --> :equal
    elseif axis_buffer < 0 # equal axes
        ahw = 1.2 * 0.5 * maximum(v -> maximum(v)-minimum(v), xyz)
        xcenter = mean(extrema(x))
        #xlims --> (xcenter-ahw, xcenter+ahw)
        ycenter = mean(extrema(y))
        #ylims --> (ycenter-ahw, ycenter+ahw)
        if _3d
            zcenter = mean(extrema(z))
            #zlims --> (zcenter-ahw, zcenter+ahw)
        end
    else
        xlims = ignorenan_extrema(x)
        if method != :chorddiagram && numnothing > 0
            x .-= mean(x)
            x /= (xlims[2] - xlims[1])
            y .-= mean(y)
            ylims = ignorenan_extrema(y)
            y /= (ylims[2] - ylims[1])
        end
        xlims --> extrema_plus_buffer(x, axis_buffer)
        ylims --> extrema_plus_buffer(y, axis_buffer)
        if _3d
            if method != :chorddiagram && numnothing > 0
                zlims = ignorenan_extrema(z)
                z .-= mean(z)
                z /= (zlims[2] - zlims[1])
            end
            zlims --> extrema_plus_buffer(z, axis_buffer)
        end
    end
    # center and rescale to the widest of all dimensions
    xyz = _3d ? (x,y,z) : (x,y)
    # Get the coordinates for the edges of the nodes.
    node_vec_vec_xy = []
    nodewidth = 0.0
    nodewidth_array = Vector{Float64}(undef, length(x))
    if !_3d
        for i in 1:length(x)
            node_weight = isnothing(node_weights) ? 1 : (10 + 100node_weights[i]/sum(node_weights))/50
            xextent, yextent = if isempty(names)
                [x[i] .+ [-0.5nodesize*node_weight, 0.5nodesize*node_weight], y[i] .+ [-0.5nodesize*node_weight, 0.5nodesize*node_weight]]
            else
                annotation_extent(plotattributes,
                                  (x[i], y[i], names[ifelse(i % length(names) == 0, length(names),
                                                            i % length(names))], fontsize*nodesize*node_weight))
            end
            nodewidth = xextent[2] - xextent[1]
            nodewidth_array[i] = nodewidth
            if nodeshape == :circle
                push!(node_vec_vec_xy, partialcircle(0, 2π, [x[i], y[i]],
                                                     80, nodewidth/2))
            elseif (nodeshape == :rect) || (nodeshape == :rectangle)
                push!(node_vec_vec_xy, [(xextent[1],yextent[1]),
                                        (xextent[2],yextent[1]),
                                        (xextent[2],yextent[2]),
                                        (xextent[1],yextent[2]),
                                        (xextent[1],yextent[1])])
            elseif nodeshape == :hexagon
                push!(node_vec_vec_xy, partialcircle(0, 2π, [x[i], y[i]],
                                                     7, nodewidth/2))
            elseif nodeshape == :ellipse
                nodeheight = (yextent[2] - yextent[1])
                push!(node_vec_vec_xy, partialellipse(0, 2π, [x[i], y[i]],
                                                      80, nodewidth/2, nodeheight/2))
            else
                error("Unknown nodeshape: $(nodeshape). Choose from :circle, ellipse, :hexagon, :rect or :rectangle.")
            end
        end
    else
        @assert _3d # TODO Make 3d work.
    end
    # create a series for the line segments
    if get(plotattributes, :linewidth, 1) > 0
        # generate a list of colors, one per segment
        segment_colors = get(plotattributes, :linecolor, nothing)
        edge_label_array = Vector{Tuple}()
        edge_label_box_vertices_array = Vector{Array}()
        if !isa(edgelabel, Dict) && !isnothing(edgelabel)
            tmp = Dict()
            if length(size(edgelabel)) < 2
                matrix_size = round(Int, sqrt(length(edgelabel)))
                edgelabel = reshape(edgelabel, matrix_size, matrix_size)
            end
            for i in 1:size(edgelabel)[1]
                for j in 1:size(edgelabel)[2]
                    if islabel(edgelabel[i, j])
                        tmp[(i, j)] = edgelabel[i, j]
                    end
                end
            end
            edgelabel = tmp
        end
        # If the edgelabel dictionary is full of length two tuples, then make all of the
        # tuples length three with last element 1. (i.e. a multigraph that has no extra
        # edges).
        if edgelabel isa Dict
            for key in keys(edgelabel)
                if length(key) == 2
                    edgelabel[(key..., 1)] = edgelabel[key]
                end
            end
        end
        edge_has_been_seen = Dict()
        for edge in zip(source, destiny)
            edge_has_been_seen[edge] = 0
        end
        for (i, (si, di, wi)) in enumerate(zip(source, destiny, weights))
            edge_has_been_seen[(si, di)] += 1
            @series begin
                xseg = Vector{Float64}()
                yseg = Vector{Float64}()
                zseg = Vector{Float64}()
                l_wg = Vector{Float64}()

                # TO DO : Colouring edges by weight
                # add a line segment
                xsi, ysi, xdi, ydi = shorten_segment(x[si], y[si], x[di], y[di], shorten)
                # For directed graphs, shorten the line segment so that the edge ends at
                # the perimeter of the destiny node.
                θ = (edge_has_been_seen[(si, di)] - 1)*pi/8
                if isdirected && si != di && !_3d
                    α = atan(y[si] - y[di], x[si] - x[di])
                    if sign(si - di) < 0
                        α = x[si] < x[di] ? θ + α : α - θ
                    else
                        α = x[si] > x[di] ? θ + α : α - θ
                    end
                    xdi = x[di] + nodewidth_array[di]*cos(α)/2
                    ydi = y[di] + nodewidth_array[di]*sin(α)/2
                    if nodeshape != :circle
                        xsi, ysi, xdi, ydi = nearest_intersection(x[si], y[si], xdi, ydi,
                                                                  node_vec_vec_xy[di])
                    end
                end
                if curves
                    if method in (:tree, :buchheim)
                        # for trees, shorten should be on one axis only
                        # dist = sqrt((x[di]-x[si])^2 + (y[di]-y[si])^2) * shorten
                        dist = shorten * (root in (:left,:bottom) ? 1 : -1)
                        ishoriz = root in (:left,:right)
                        xsi, xdi = (ishoriz ? (x[si]+dist,x[di]-dist) : (x[si],x[di]))
                        ysi, ydi = (ishoriz ? (y[si],y[di]) : (y[si]+dist,y[di]-dist))

                        xpts, ypts = directed_curve(xsi, xdi, ysi, ydi,
                                    xview=get(plotattributes, :xlims, (0,1)), yview=get(plotattributes, :ylims, (0,1)), root=root)

                        append!(xseg, push!(xpts, NaN))
                        append!(yseg, push!(ypts, NaN))
                        append!(l_wg, [ wi for i in 1:length(xpts) ] )
                    elseif method == :arcdiagram
                        r  = (xdi - xsi) / 2
                        x₀ = (xdi + xsi) / 2
                        θ = range(0, stop=π, length=30)
                        xpts = x₀ .+ r .* cos.(θ)
                        ypts = r .* sin.(θ) .+ ysi # ysi == ydi
                        for x in xpts
                            push!(xseg, x)
                            push!(l_wg, wi)
                        end
                        push!(xseg, NaN)
                        for y in ypts
                            push!(yseg, y)
                        end
                        push!(yseg, NaN)
                    else
                        xpt, ypt = if method != :chorddiagram
                            control_point(xsi, x[di],
                                          ysi, y[di],
                                          edge_has_been_seen[(si, di)]*curvature_scalar*sign(si - di))
                        else
                            (0.0, 0.0)
                        end
                        xpts = [xsi, xpt, xdi]
                        ypts = [ysi, ypt, ydi]
                        t = range(0, stop=1, length=3)
                        A = hcat(xpts, ypts)

                        itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)

                        tfine = range(0, 1, length=30)
                        xpts, ypts = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
                        if !isnothing(edgelabel) && haskey(edgelabel, (si, di, edge_has_been_seen[(si, di)]))
                            q = control_point(xsi, x[di],
                                              ysi, y[di],
                                              (edgelabel_offset
                                              + edge_has_been_seen[(si, di)]*curvature_scalar)*sign(si - di))
                            push!(edge_label_array,
                                  (q...,
                                   string(edgelabel[(si, di, edge_has_been_seen[(si, di)])]), fontsize))
                            edge_label_box_vertices = (
                            annotation_extent(plotattributes,
                                              (q[1], q[2],
                                              edgelabel[(si, di, edge_has_been_seen[(si, di)])],
                                              0.05fontsize)))
                            if !any(isnan.(q))
                                push!(edge_label_box_vertices_array, edge_label_box_vertices)
                            end
                        end
                        if method != :chorddiagram && !_3d
                            append!(xseg, push!(xpts, NaN))
                            append!(yseg, push!(ypts, NaN))
                            push!(l_wg, wi)
                        else
                            push!(xseg, xsi, xpt, xdi, NaN)
                            push!(yseg, ysi, ypt, ydi, NaN)
                            _3d && push!(zseg, z[si], z[si], z[di], NaN)
                            push!(l_wg, wi)
                        end
                    end
                else
                    push!(xseg, xsi, xdi, NaN)
                    push!(yseg, ysi, ydi, NaN)
                    _3d && push!(zseg, z[si], z[di], NaN)
                end
            if si == di && !_3d
                inds = 1:n .!= si
                self_edge_angle = pi/8 + (edge_has_been_seen[(si, di)] - 1)*pi/8
                θ1 = unoccupied_angle(xsi, ysi, x[inds], y[inds]) - self_edge_angle/2
                θ2 = θ1 + self_edge_angle
                nodewidth = nodewidth_array[si]
                if nodeshape == :circle
                    xpts = [xsi + nodewidth*cos(θ1)/2,
                            NaN, NaN, NaN,
                            xsi + nodewidth*cos(θ2)/2]
                    xpts[2] = mean([xpts[1], xpts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*cos(θ1)
                    xpts[3] = mean([xpts[1], xpts[end]]) + edge_has_been_seen[(si, di)]*self_edge_size*cos((θ1 + θ2)/2)
                    xpts[4] = mean([xpts[1], xpts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*cos(θ2)
                    ypts = [ysi + nodewidth*sin(θ1)/2,
                            NaN, NaN, NaN,
                            ysi + nodewidth*sin(θ2)/2]
                    ypts[2] = mean([ypts[1], ypts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*sin(θ1)
                    ypts[3] = mean([ypts[1], ypts[end]]) + edge_has_been_seen[(si, di)]*self_edge_size*sin((θ1 + θ2)/2)
                    ypts[4] = mean([ypts[1], ypts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*sin(θ2)

                    t = range(0, stop=1, length=5)
                    A = hcat(xpts, ypts)

                    itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)

                    tfine = range(0, 1, length=50)
                    xpts, ypts = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
                else
                    _, _,
                    start_point1,
                    start_point2 = nearest_intersection(xsi, ysi,
                                                        xsi + 2nodewidth*cos(θ1),
                                                        ysi + 2nodewidth*sin(θ1),
                                                        node_vec_vec_xy[si])
                    _, _, end_point1,
                    end_point2 = nearest_intersection(xsi + edge_has_been_seen[(si, di)]*(nodewidth + self_edge_size)*cos(θ2),
                                                      ysi + edge_has_been_seen[(si, di)]*(nodewidth + self_edge_size)*sin(θ2),
                                                      xsi,
                                                      ysi,
                                                      node_vec_vec_xy[si])

                    xpts = [start_point1,
                            NaN, NaN, NaN,
                            end_point1]
                    xpts[2] = mean([xpts[1], xpts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*cos(θ1)
                    xpts[3] = mean([xpts[1], xpts[end]]) + edge_has_been_seen[(si, di)]*self_edge_size*cos((θ1 + θ2)/2)
                    xpts[4] = mean([xpts[1], xpts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*cos(θ2)
                    ypts = [start_point2,
                            NaN, NaN, NaN,
                            end_point2]
                    ypts[2] = mean([ypts[1], ypts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*sin(θ1)
                    ypts[3] = mean([ypts[1], ypts[end]]) + edge_has_been_seen[(si, di)]*self_edge_size*sin((θ1 + θ2)/2)
                    ypts[4] = mean([ypts[1], ypts[end]]) + 0.5*(0.5 + edge_has_been_seen[(si, di)])*self_edge_size*sin(θ2)

                    t = range(0, stop=1, length=5)
                    A = hcat(xpts, ypts)

                    itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)

                    tfine = range(0, 1, length=50)
                    xpts, ypts = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
                end
                append!(xseg, push!(xpts, NaN))
                append!(yseg, push!(ypts, NaN))
                mid_ind = div(length(xpts), 2)
                q = [xpts[mid_ind] + edgelabel_offset*cos((θ1 + θ2)/2),
                     ypts[mid_ind] + edgelabel_offset*sin((θ1 + θ2)/2)]
                if !isnothing(edgelabel) && haskey(edgelabel, (si, di, edge_has_been_seen[(si, di)]))
                    push!(edge_label_array,
                          (q...,
                           string(edgelabel[(si, di, edge_has_been_seen[(si, di)])]), fontsize))
                    edge_label_box_vertices = annotation_extent(plotattributes, (q...,
                                                                edgelabel[(si, di, edge_has_been_seen[(si, di)])],
                                                                0.05fontsize))
                    if !any(isnan.(q))
                        push!(edge_label_box_vertices_array, edge_label_box_vertices)
                    end
                end
            end

            if isa(segment_colors, ColorGradient)
                line_z := segment_colors[i]
            end
            linewidthattr = get(plotattributes, :linewidth, 1)
            seriestype := if method in (:tree, :buchheim, :chorddiagram)
                :curves
            else
                if _3d
                    if curves
                        :curves
                    else
                        :path3d
                    end
                else
                    :path
                end
            end
            linewidth --> linewidthattr * edgewidth(si, di, wi)
            markershape := :none
            markercolor := :black
            isdirected && (arrow --> :simple, :head, 0.3, 0.3)
            primary := false
            _3d ? (xseg, yseg, zseg) : (xseg, yseg)
            end

        end
    end
    # The boxes around edge labels are defined as another list of series that sits on top
    # of the series for the edges.
    edge_has_been_seen = Dict()
    for edge in zip(source, destiny)
        edge_has_been_seen[edge] = 0
    end
    if edge_label_box
        index = 0
        for (i, (si, di, wi)) in enumerate(zip(source, destiny, weights))
            edge_has_been_seen[(si, di)] += 1
            if !isnothing(edgelabel) && haskey(edgelabel, (si, di, edge_has_been_seen[(si, di)]))
                index += 1
                @series begin
                    seriestype := :shape
                    fillcolor --> get(plotattributes, :background_color, :white)
                    linewidth := 0
                    linealpha := 0
                    edge_label_box_vertices = edge_label_box_vertices_array[index]
                    ([edge_label_box_vertices[1][1], edge_label_box_vertices[1][2],
                      edge_label_box_vertices[1][2], edge_label_box_vertices[1][1],
                      edge_label_box_vertices[1][1]],
                     [edge_label_box_vertices[2][1], edge_label_box_vertices[2][1],
                      edge_label_box_vertices[2][2], edge_label_box_vertices[2][2],
                      edge_label_box_vertices[2][1]])
                end
            end
        end
    end

    framestyle := :none
    axis := nothing
    legend --> false
    if method == :chorddiagram
        seriestype := :scatter
        markersize := 0
        markeralpha := 0
        ratio --> :equal
        if length(names) == length(x)
            annotations := [(x[i], y[i], names[i]) for i in 1:length(x)]
        end
        @series begin
            seriestype := :shape
            N = length(x)
            angles = Vector{Float64}(undef, N)
            for i in 1:N
                if y[i] > 0
                    angles[i] = acos(x[i])
                else
                    angles[i] = 2pi - acos(x[i])
                end
            end
            δ = 0.4 * (angles[2] - angles[1])
            vec_vec_xy = [ arcshape(Θ-δ,Θ+δ) for Θ in angles ] # Shape
            [ [ xy[1] for xy in vec_xy ] for vec_xy in vec_vec_xy ], [ [ xy[2] for xy in vec_xy ] for vec_xy in vec_vec_xy ]
        end
    else
        if isempty(names)
            if _3d
                seriestype := :scatter3d
                linewidth := 0
                linealpha := 0
                series_annotations --> map(string,names)
                markersize --> (10 .+ (100 .* node_weights) ./ sum(node_weights))
            else
                for (i, vec_xy) in enumerate(node_vec_vec_xy)
                    @series begin
                        if !_3d
                            seriestype := :shape
                            tmp = nodecolor isa AbstractArray ? nodecolor[i] : nodecolor
                            fillcolor --> tmp
                            linewidth := 0
                            linealpha := 0
                            ([xy[1] for xy in vec_xy], [xy[2] for xy in vec_xy])
                        else  # TODO make 3d work.
                            seriestype := :volume
                            ([[xyz[1] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz],
                             [[xyz[2] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz],
                             [[xyz[3] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz])
                        end
                    end
                    seriestype := :scatter
                    markersize := 0
                    markeralpha := 0
                end
                !isnothing(edgelabel) && (annotation := edge_label_array)
            end
        else
            @assert !_3d  # TODO: make this work in 3D
            for (i, vec_xy) in enumerate(node_vec_vec_xy)
                @series begin
                    if !_3d
                        seriestype := :shape
                        tmp = nodecolor isa AbstractArray ? nodecolor[i] : nodecolor
                        fillcolor --> tmp
                        linewidth := 0
                        linealpha := 0
                        ([xy[1] for xy in vec_xy], [xy[2] for xy in vec_xy])
                    else  # TODO make 3d work.
                        seriestype := :volume
                        ([[xyz[1] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz],
                         [[xyz[2] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz],
                         [[xyz[3] for xyz in vec_xyz] for vec_xyz in node_vec_vec_xyz])
                    end
                end
            end
            seriestype := :scatter
            markersize := 0
            markeralpha := 0
            annotations := [edge_label_array ; [(x[i], y[i],
                                                names[ifelse(i % length(names) == 0,
                                                              length(names),
                                                              i % length(names))],
                                                fontsize) for i in 1:length(x)]]
        end
    end
    xyz
end

@recipe function f(g::AbstractGraph)
    GraphPlot(get_source_destiny_weight(get_adjacency_list(g)))
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

using Documenter, FEMCoupling

makedocs(modules=[FEMCoupling],
         format = :html,
         checkdocs = :all,
         sitename = "FEMCoupling.jl",
         pages = ["index.md"]
        )


text:
"""
Opaque data structure holding the whole set of parameters used to compute a routing.

# Internals

These are implementation details that could evolve at any time. Many parameters here are already documented in the main
public constructor for `RoutingData`.

The first four fields are the most dependable of all of them (i.e. they are likely to stay unchanged for a long period of time):

* `g`: topology.
* `k`: demands.
* `solver`.
* `name`.

Several fields configure the solving process.

* Some parameters are always used:
  * `model_type`.
  * `model_simplifications`. TODO: reintroduce them.
* Some parameters are sometimes used:
  * `sub_model_type`.
* Some parameters are only useful for oblivious routing:
  * `model_all_traffic_matrices`: only for cutting-plane implementation.
  * `model_exact_opt_d`: only for cutting-plane implementation.
  * `model_robust_reformulation_traffic_matrices`: only for dual-reformulation implementation.
* Some parameters are only useful for problems without uncertainty.
  * `traffic_matrix`.

Several fields configure the output process.

* `logmessage`
* `plot_final_results`: TODO: reintroduce.
* `plot_each_iteration`: TODO: reintroduce.
* `export_lps`
* `export_lps_on_error`
* `output_folder`

From these graphs, several data structures are derived. Some of them are lazily computed when required (hence the
`mutable` aspect of this structure).

* For path-based formulations, paths are stored and uniquely identified:
  * `paths_edges`: all the paths, represented as list of edges (Edges). Path indices are always consistent
    with `paths`. Used as follows: `paths_edges[path_id]` is a list of edges.
  * `demand_to_path_ids`: maps a demand (arc from source to destination) to its path IDs.
  * `path_id_to_demand`: maps a path ID to the corresponding demand.
* In order to draw the graph, node positions are computed at most once:
  * `locs_x`: for each node, its x coordinate when drawing the graph (lazily created).
  * `locs_y`: for each node, its y coordinate when drawing the graph (lazily created).
  * `locs_f`: function that computes the two previous fields when requested. Takes a graph as input and outputs
    two vectors, `locs_x` and `locs_y`.

The time to create the `RoutingData` object is memorised in `time_precompute_ms`.
"""
mutable struct RoutingData
    g::MetaDiGraph
    k::AbstractGraph
    solver
    name::String
    time_precompute_ms::Float64

    model_type::ModelType
    sub_model_type::ModelType
    model_simplifications::Bool
    model_all_traffic_matrices::Bool
    model_exact_opt_d::Bool
    model_robust_reformulation_traffic_matrices::Bool

    traffic_matrix::Dict{Edge{Int}, Float64}

    paths_edges::Vector{Vector{Edge}}
    demand_to_path_ids::Dict{Edge{Int}, Vector{Int}}
    path_id_to_demand::Dict{Int, Edge}

    locs_x::Union{Vector{Float64}, Nothing}
    locs_y::Union{Vector{Float64}, Nothing}
    locs_f

    logmessage
    plot_final_results
    plot_each_iteration
    export_lps
    export_lps_on_error
    output_folder
end

"""
Creates an opaque `RoutingData` object based on the topology of the network `g` and the demands to route `k`.
Calling this function may take a while, as it precomputes many things used later on, based on the type of model to use.
Not all parameters are used in all cases.

Parameters that are always used (only the first three are required and positional):

* `g` is a directed graph with annotation (hence the type `MetaDiGraph`). Each edge is supposed to have a property
  `capacity` giving its capacity (units are not imposed by the package, but must be consistent). Each node must
  have a `name` property, which will be used for generating the output.
* `k` is another directed graph, without necessarily annotations (hence the type `AbstractSimpleGraph`).
* `solver` is the optimization solver to use for this problem (like `CPLEX.Optimizer`).
* `model_type` is the type of optimisation model used to compute the oblivious routing. Its type is `ModelType`.
* `name` is an optional name for the model.

The parameters almost always apply, and may have tremendous impact on solving performance:

* `model_simplifications` enables simplifications in the model. They do not imply any kind of approximation.
  Computing what parts of the model can be simplified may take more time than what the simplifications may improve
  the running times of the optimalisation part (especially for small networks and those with very high connexity).
* `remove_unreachable_nodes` determines whether the nodes in the topology `g` that are isolated, i.e. those that
  have no neighbours, should be removed
* `remove_unsatisfiable_demands` determines whether the demands that are not routable (i.e. there is no path in `g`
  from their source to their destination, the topology `g` provides no connectivity between these two nodes) should
  be removed.

Some models may use these parameters:

* `sub_model_type` is only used by algorithms using subproblems, like oblivious routing with cutting planes.
  All the following parameters are reused for the subproblems: `model_simplifications`, `npaths`.
* `npaths` is the number of shortest paths that should be computed beforehand for path-based formulations, solved with
  column generation or not. This number of paths is considered *for each demand*: there will be
  probably more than `npaths` that get generated, but never more than `npaths` times the number of demands).
  By default, 20 paths at most are generated for each demand.
  For column-generation-based formulations, this provides a starting set of columns: this set may be increased in size
  during the solving process. Otherwise, this is the total number of paths that will be used for solving the problem.

These parameters are only used for oblivious routing:

* `model_all_traffic_matrices` indicates whether all traffic matrices that violate the oblivious routing constraints
  should be added or not, at each iteration of the master problem. (The original algorithm only adds one matrix,
  but adding more at a time may improve running times by reducing the total number of iterations to perform.)
* `model_exact_opt_d` indicates whether, in the added constraints, the optimal congestion is computed for the generated
  traffic matrix. In theory, this value is always 1.0.
* `model_robust_reformulation_traffic_matrices` enables computing traffic matrices for the robust reformulations.
  Using this options increases the number of constraints in the problem to solve in order to get the right dual values.

These parametesr are only used for problems without uncertainty.

* `traffic_matrix` is the only traffic matrix that is considered.

These parameters only control the output of the solving process:

* `locs_f` function that computes the position of the nodes when plotting the graph in 2D. It takes a graph as input and
  outputs two vectors, first the x positions, then the y positions.
* `logmessage` is a function that will be called whenever the package produces any textual output. By default,
  text is redirected to the shell (`println`). This function should take a single string argument, the message to print.
* `verbose` indicates whether this function (not the complete package) should be explicit about what it is doing
  (for instance, if it removes demands). This option does not include regular progress logging for long operations.
* `plot_final_results`: whether plots should be produced when the solving process is done.
* `plot_each_iteration`: for iterative algorithms, whether plots should be produced at each iteration.
* `export_lps`: export the LP files (if possible) for each and every (sub)problem that is solved. If the master problem
  is iteratively built, each iteration will be exported. As the LP format can only represent MILP problems, not all
  combinations of problem types and solving algorithms can benefit from this parameter.
* `export_lps_on_error`: export the LP files (if possible) of the (sub)problem being solved when an error occurs.
  As the LP format can only represent MILP problems, not all combinations of problem types and solving algorithms
  can benefit from this parameter.
* `output_folder`: the folder where all file output is supposed to be taking place. This folder must exist beforehand
  or be an empty string to disable all output.
"""
function RoutingData(g::AbstractMetaGraph, # Not MetaDiGraph, to show a more interesting error message if the graph is not directed.
                     k::AbstractGraph,
                     solver,
                     model_type::ModelType;
                     name::String="",
                     sub_model_type::Union{ModelType, Nothing}=nothing,
                     model_simplifications::Bool=false,
                     model_all_traffic_matrices::Bool=false,
                     model_exact_opt_d::Bool=false,
                     model_robust_reformulation_traffic_matrices::Bool=false,
                     traffic_matrix::Dict{Edge{Int}, Float64}=Dict{Edge{Int}, Float64}(),
                     npaths::Int=20,
                     remove_unreachable_nodes::Bool=true,
                     remove_unsatisfiable_demands::Bool=true,
                     locs_f=spring_layout,
                     logmessage=println,
                     verbose::Bool=true,
                     plot_final_results::Bool=false,
                     plot_each_iteration::Bool=false,
                     export_lps::Bool=false,
                     export_lps_on_error::Bool=false,
                     output_folder::String="")
    start = time_ns()

    if sub_model_type === nothing
        sub_model_type = copy(model_type)
    end

    # Check whether the right properties have been defined in the graphs.
    for v in vertices(g)
        if ! has_prop(g, v, :name)
            error("Node $(v) has no name!")
        end
    end
    for e in edges(g)
        if ! has_prop(g, e, :capacity)
            error("Edge $(e) has no capacity!")
        end
    end

    if ! is_directed(g)
        error("The demand graph is not directed!")
    end
    if ! is_directed(k)
        error("The demand graph is not directed!")
    end

    # Check whether the other parameters have easily-checkable properties, before more cumbersome computations.
    if output_folder != "" && ! isdir(output_folder)
        error("The output path $(output_folder) does not exist!")
    end

    # Copy the graphs locally, so that the caller does not see when they are modified internally
    # (for instance, due to remove_unsatisfiable_demands).
    g = copy(g)
    k = copy(k)

    # If need be, remove some vertices because they are isolated.
    if remove_unreachable_nodes
        nremoved = remove_unreachable_nodes!(g)
        if verbose && nremoved > 0
            logmessage("Removed $(nremoved) vertices because they were not reachable.")
        end
    end

    # If need be, remove some demands due to lack of connectivity.
    if remove_unsatisfiable_demands
        nremoved = remove_unsatisfiable_demands!(g, k)
        if verbose && nremoved > 0
            logmessage("Removed $(nremoved) demands because they were not routable.")
        end
    end

    # If this is a path-based formulation, it may require some path generation.
    if (model_type.type == PathFormulation() || sub_model_type.type == PathFormulation()) && npaths > 0
        # Precompute some shortest paths and directly assign them numbers.
        # A path is a sequence of integers, each being the index of a node.
        paths = Vector{Int}[]
        demand_to_path_ids = Dict{Edge{Int}, Vector{Int}}()
        path_id_to_demand = Dict{Int, Edge}()

        for d in edges(k)
            demand_paths = yen_k_shortest_paths(g, src(d), dst(d), weights(g), npaths).paths
            if ! (model_type.cg || sub_model_type.cg) && length(demand_paths) == npaths
                @warn("The limit of paths to generate ($(npaths)) has been reached; the result will not be guaranteed to be optimum. " *
                      "Switch to a column-generation formulation for exact results.")
            end

            demand_to_path_ids[d] = Int[]
            for path in demand_paths
                push!(paths, path)
                p_id = length(paths)
                push!(demand_to_path_ids[d], p_id)
                path_id_to_demand[p_id] = d
            end
        end

        # Paths are returned as a list of nodes, transform it into a list of edges.
        paths_edges = [map(t -> Edge(t[1], t[2]), zip(p[1:end-1], p[2:end])) for p in paths]
    else
        # Not a path-based formulation: create empty collections of the right type.
        paths_edges = Vector{Edge}[]
        demand_to_path_ids = Dict{Edge{Int}, Vector{Int}}()
        path_id_to_demand = Dict{Int, Edge}()
    end

    stop = time_ns()

    return RoutingData(g, k, solver, name, (stop - start) / 1_000_000.,
                       # Solving parameters.
                       model_type, sub_model_type, model_simplifications,
                       model_all_traffic_matrices, model_exact_opt_d,
                       model_robust_reformulation_traffic_matrices,
                       traffic_matrix,
                       # Precomputed things.
                       paths_edges, demand_to_path_ids, path_id_to_demand,
                       # Output parameters (graph plotting and logging).
                       nothing, nothing, locs_f, logmessage, plot_final_results,
                       plot_each_iteration, export_lps, export_lps_on_error, output_folder)
end

function remove_unreachable_nodes!(g::AbstractGraph)
    # First, determine which nodes are unreachable; then, delete them.
    # If deletion is made in the loop over all vertices, as soon as a deletion
    # is performed, the vertex indices may change: vertices(g) is no more
    # consistent with inneighbors(g, v) and outneighbors(g, v).
    to_remove = Int[]
    for v in vertices(g)
        if length(inneighbors(g, v)) == 0 && length(outneighbors(g, v)) == 0
            push!(to_remove, v)
        end
    end
    for v in to_remove
        rem_vertex!(g, v)
    end
    return length(to_remove)
end

function remove_unsatisfiable_demands!(g::AbstractGraph, k::AbstractGraph)
    to_remove = Edge{Int}[]
    for e in edges(k)
        paths = dijkstra_shortest_paths(g, src(e), weights(g))
        if length(paths.parents) < dst(e) || paths.parents[dst(e)] == 0
            push!(to_remove, e)
        end
    end
    for e in to_remove
        rem_edge!(k, e)
    end
    return length(to_remove)
end

function add_path!(rd::RoutingData, demand::Edge, path::Vector{Edge})
    # Update RoutingData with the new path.
    # All demands already have a path, by assumption. Otherwise, the
    # master problem should not be feasible. Hence, computing the new
    # path ID is easy: just take the length of the vector.
    push!(rd.paths_edges, path)
    path_id = length(rd.paths_edges)

    rd.path_id_to_demand[path_id] = demand
    push!(rd.demand_to_path_ids[demand], path_id)

    return path_id
end

graph(rd::RoutingData) = rd.g
edges(rd::RoutingData) = edges(rd.g)
vertices(rd::RoutingData) = vertices(rd.g)
demands(rd::RoutingData) = edges(rd.k)

n_nodes(rd::RoutingData) = nv(rd.g)
n_edges(rd::RoutingData) = ne(rd.g)
n_demands(rd::RoutingData) = ne(rd.k)
n_paths(rd::RoutingData) = length(rd.paths_edges)

capacity(rd::RoutingData, e::Edge{Int}) = get_prop(rd.g, e, :capacity)

function find_path_ids_with_edge(rd::RoutingData, e::Edge)
    if n_paths(rd) > 0
        return filter(1:n_paths(rd)) do p
            e in rd.paths_edges[p]
        end
    else
        error("Precomputed paths are not available. Are you trying to use this function in a flow-based formulation?")
    end
end


text:
using MTFun
using Test
using ApproxFun

@testset "Constructors, Evaluation, Transforms" begin
    λ = -0.7+0.9im
    MT = MalmquistTakenaka(λ)
    φ = (x,n) -> (-1.0im)^(n+1)*sqrt(abs(imag(λ))/π) * (λ - x)^n / (conj(λ) - x)^(n+1)
    @test coefficients(Fun(x->φ(x,0),MT,12)) ≈ [1;zeros(11)]
    @test coefficients(Fun(x->φ(x,5),MT,12)) ≈ [zeros(10);1;0]
    @test coefficients(Fun(x->φ(x,0),MT,11)) ≈ [1;zeros(10)]
    @test coefficients(Fun(x->φ(x,5),MT,11)) ≈ [zeros(10);1]

    f = x -> (1.0-2im)/(1+x^2)
    F = Fun(f,MT)
    x = 0.78
    @test F(x) ≈ f(x)

    f = x -> exp(-10x^2)
    F = Fun(f,MT)
    x = 0.78
    @test real(F(x)) ≈ f(x)

    vals = φ.(points(MT,10),0)
    cfs = [1.0+0.0im;zeros(9)]
    trans = MTFun.plan_transform(MT,vals)
    itrans = MTFun.plan_itransform(MT,cfs)
    @test cfs ≈ trans*vals
    @test vals ≈ itrans*cfs

    v = [0.83+0.0im;0.78;-1.13;0.65;1.11]
    trans = MTFun.plan_transform(MT,v)
    itrans = MTFun.plan_itransform(MT,v)
    @test v ≈ itrans*(trans*v)
    @test v ≈ trans*(itrans*v)

    W = Weideman(λ)
    φ = (x,n) -> (-1.0im)^n * (λ - x)^n / (conj(λ) - x)^n
    @test coefficients(Fun(x->φ(x,0),W,12)) ≈ [1;zeros(11)]
    @test coefficients(Fun(x->φ(x,5),W,12)) ≈ [zeros(10);1;0]
    @test coefficients(Fun(x->φ(x,0),W,11)) ≈ [1;zeros(10)]
    @test coefficients(Fun(x->φ(x,5),W,11)) ≈ [zeros(10);1]

    f = x -> (1.0-2im)/(1+x^2)
    F = Fun(f,W)
    x = 0.78
    @test F(x) ≈ f(x)

    f = x -> exp(-10x^2)
    F = Fun(f,W)
    x = 0.78
    @test real(F(x)) ≈ f(x)

    vals = φ.(points(W,10),0)
    cfs = [1.0+0.0im;zeros(9)]
    trans = MTFun.plan_transform(W,vals)
    itrans = MTFun.plan_itransform(W,cfs)
    @test cfs ≈ trans*vals
    @test vals ≈ itrans*cfs

    v = [0.83+0.0im;0.78;-1.13;0.65;1.11]
    trans = MTFun.plan_transform(W,v)
    itrans = MTFun.plan_itransform(W,v)
    @test v ≈ itrans*(trans*v)
    @test v ≈ trans*(itrans*v)

end

@testset "Operators" begin

    λ = -0.7+0.9im
    MT = MalmquistTakenaka(λ)

    f = x -> (1.0-2im)/(1+x^2)
    fprime = x -> -(1.0 - 2im)*2x/(1+x^2)^2
    F = Fun(f,MT)
    Fprime = differentiate(F)
    x = 0.78
    @test Fprime(x) ≈ fprime(x)
    @test (Derivative(MT,2)*F)(x) ≈ (differentiate(Fprime))(x)

    f = x -> exp(-10x^2)
    fprime = x -> -20*x*exp(-10x^2)
    F = Fun(f,MT)
    Fprime = differentiate(F)
    x = 0.78
    @test Fprime(x) ≈ fprime(x)
    @test (Derivative(MT,2)*F)(x) ≈ (differentiate(Fprime))(x)

    g = x -> sech(x)^2
    G = Fun(g,MT)
    M = Multiplication(G,MT)
    GF = M*F
    x = 0.78
    @test GF(x) ≈ g(x) * f(x)

    W = Weideman(λ)

    f = x -> (1.0-2im)/(1+x^2)
    F = Fun(f,MT)
    g = x -> sech(x)^2
    G = Fun(g,W)

    M = Multiplication(G,MT)
    GF = M*F
    x = 0.78
    @test GF(x) ≈ g(x) * f(x)

    F = Fun(f,W)
    M = Multiplication(G,W)
    GF = M*F
    x = 0.78
    @test GF(x) ≈ g(x) * f(x)

end


text:
include("activate.jl")

using LinearAlgebra
using Random
import GLPK
using JSON

n = 200
k = 3000
#k = 500

xpi = rand(n * n);
total = sum(xpi);
# next line needs to be commented out if we use the GLPK variants
xpi = reshape(xpi, n, n)
const xp = xpi # ./ total;

# better for memory consumption as we do coordinate-wise ops

function cf(x, xp)
    return LinearAlgebra.norm(x .- xp)^2 / n^2
end

function cgrad!(storage, x, xp)
    return @. storage = 2 * (x - xp) / n^2
end

# initial direction for first vertex
direction_vec = Vector{Float64}(undef, n * n)
randn!(direction_vec)
direction_mat = reshape(direction_vec, n, n)

lmo = FrankWolfe.BirkhoffPolytopeLMO()
x00 = FrankWolfe.compute_extreme_point(lmo, direction_mat)
target_accuracy = 1e-7

# modify to GLPK variant
# o = GLPK.Optimizer()
# lmo_moi = FrankWolfe.convert_mathopt(lmo, o, dimension=n)
# x00 = FrankWolfe.compute_extreme_point(lmo, direction_vec)

FrankWolfe.benchmark_oracles(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    () -> randn(n, n),
    lmo;
    k=100,
)

# vanllia FW

x0 = copy(x00)

x, v, primal, dual_gap, trajectoryFW = FrankWolfe.frank_wolfe(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    epsilon=target_accuracy,
    emphasis=FrankWolfe.memory,
    trajectory=true,
    verbose=true,
);

# arbitrary cache

x0 = copy(x00)

x, v, primal, dual_gap, trajectoryLCG = FrankWolfe.lazified_conditional_gradient(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    lmo,
    x0,
    max_iteration=k,
    epsilon=target_accuracy,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    trajectory=true,
    verbose=true,
);

# fixed cache size

x0 = copy(x00)

x, v, primal, dual_gap, trajectoryBLCG = FrankWolfe.lazified_conditional_gradient(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    epsilon=target_accuracy,
    emphasis=FrankWolfe.memory,
    trajectory=true,
    cache_size=500,
    verbose=true,
);

# AFW run

x0 = copy(x00)

x, v, primal, dual_gap, trajectoryLAFW = FrankWolfe.away_frank_wolfe(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    linesearch_tol=1e-9,
    epsilon=target_accuracy,
    emphasis=FrankWolfe.memory,
    lazy=true,
    trajectory=true,
    verbose=true,
);

# BCG run

x0 = copy(x00)

x, v, primal, dual_gap, trajectoryBCG = FrankWolfe.blended_conditional_gradient(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    linesearch_tol=1e-9,
    epsilon=target_accuracy,
    emphasis=FrankWolfe.memory,
    trajectory=true,
    verbose=true,
);

# BCG run (reference optimum)

x0 = copy(x00)

x, v, primal, dual_gap, trajectoryBCG_ref = FrankWolfe.blended_conditional_gradient(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    lmo,
    x0,
    max_iteration=2 * k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    epsilon=target_accuracy / 10.0,
    linesearch_tol=1e-9,
    emphasis=FrankWolfe.memory,
    trajectory=true,
    verbose=true,
);

open("lcg_expensive_data.json", "w") do f
    return write(
        f,
        JSON.json((
            FW=trajectoryFW,
            LCG=trajectoryLCG,
            BLCG=trajectoryBLCG,
            LAFW=trajectoryLAFW,
            BCG=trajectoryBCG,
            reference_BCG_primal=primal,
        )),
    )
end

data = [trajectoryFW, trajectoryLCG, trajectoryBLCG, trajectoryLAFW, trajectoryBCG]
label = ["FW", "L-CG", "BL-CG", "L-AFW", "BCG"]
FrankWolfe.plot_trajectories(data, label, xscalelog=true)


text:
# Test routine to compare performance for a number of SDP problems (max Eigenvalue problems)
workspace()
include("../../../src/Solver.jl")
include("../solverComparison/Compare.jl")

using OSQP, OSSDP, Base.Test, Compare

rng = Random.MersenneTwister(123734)

nn = 25
timestamp = Dates.format(now(), "yyddmm_HH-MM")
dataFile = "../resultDataFiles/SC_" * timestamp * "MaxEigSDP.jld"
problemType = "SDP-MaxEig"

sr1 = SolverResult(nn, problemType,"COSMO",timestamp,0,false,false)
sr2 = SolverResult(nn, problemType,"COSMO",timestamp,0,true,false)
sr3 = SolverResult(nn, problemType,"COSMO",timestamp,0,false,true)
sr4 = SolverResult(nn, problemType,"COSMO",timestamp,0,true,true)
resData = [sr1;sr2;sr3;sr4]
ws = 0
for iii =1:1:nn
 # generate symmetric test matrix A
  nA = rand(rng,2:35)
  A = Symmetric(randn(rng,nA,nA))

  # define max eigenvalue problem as SDP

  # create augmented matrices
  Aa = -vec(eye(nA))
  ba = -vec(A)
  q = [1]
  P = spzeros(1,1)
  # define cone membership
  K = Cone(0,0,[],[nA^2])
  pDims = [size(Aa,1);size(Aa,2);size(Aa,1)*size(Aa,2)]

  setOFF = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,check_termination=1,scaling = 0,eps_abs = 1e-3,eps_rel=1e-3,adaptive_rho=false)
  setON = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,check_termination=1,scaling = 10,eps_abs = 1e-3,eps_rel=1e-3,adaptive_rho=false)
  setOFFAd = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,check_termination=1,scaling = 0,eps_abs = 1e-3,eps_rel=1e-3,adaptive_rho=true)
  setONAd = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,check_termination=1,scaling = 10,eps_abs = 1e-3,eps_rel=1e-3,adaptive_rho=true)
  # Solve with OSSDP
  resOSSDP_unscaled,nothing = OSSDP.solve(P,q,Aa,ba,K,setOFF);
  print("\n.")
  resOSSDP_scaled,nothing = OSSDP.solve(P,q,Aa,ba,K,setON);
  print(".")
   resOSSDP_unscaledAd,nothing = OSSDP.solve(P,q,Aa,ba,K,setOFFAd);
  print("\n.")
  resOSSDP_scaledAd,nothing = OSSDP.solve(P,q,Aa,ba,K,setONAd);
  print(".")

  println("Correct max eigenvalue: $(maximum(eig(A)[1]))")

  resArray = [resOSSDP_unscaled;resOSSDP_scaled;resOSSDP_unscaledAd;resOSSDP_scaledAd]
  updateResults!(dataFile,resData,resArray,pDims,"SDP-$(iii)",0.,true)
  printStatus(iii,nn,"SDP-$(iii)",resData)
end


text:
# Autogenerated wrapper script for SHTOOLS_jll for powerpc64le-linux-gnu-libgfortran5
export libSHTOOLS, libSHTOOLS_mp

using CompilerSupportLibraries_jll
using FFTW_jll
using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("SHTOOLS")
JLLWrappers.@declare_library_product(libSHTOOLS, "libSHTOOLS.so")
JLLWrappers.@declare_library_product(libSHTOOLS_mp, "libSHTOOLS-mp.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, FFTW_jll, OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libSHTOOLS,
        "lib/libSHTOOLS.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libSHTOOLS_mp,
        "lib/libSHTOOLS-mp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module KmeansAlgorithmTask

using Plots, Statistics

export kmeansmy
export graficaclusters

include("solving_Kmeans.jl")

end # module


text:
module Cylinder

import SimpleMolecule: Molecule
import MolecularBuild: FibonacciCylinder, BuildCylinder

function assemble(systems::Array{Molecule,1}, params::Dict)
    return BuildCylinder.build(
        FibonacciCylinder,
        systems;
        number=convert(Int, params["number"]),
        radius=float(params["radius"]),
        zmin=float(params["zmin"]),
        zmax=float(params["zmax"]),
        center=Float64[float(x) for x in params["center"]],
        aligned=true,
        old_orientation=Float64[float(x) for x in params["orientation"]],
        inverted=params["inverted"]
    )
end

end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
using Pkg; Pkg.activate(".")
using Franklin, MLJ, PyPlot
FOLDER = @__DIR__


text:
const Spectral = [RGB{Float32}(0.62,0.004,0.259),
				RGB{Float32}(0.835,0.243,0.31),
				RGB{Float32}(0.957,0.427,0.263),
				RGB{Float32}(0.992,0.682,0.38),
				RGB{Float32}(0.996,0.878,0.545),
				RGB{Float32}(1.0,1.0,0.749),
				RGB{Float32}(0.902,0.961,0.596),
				RGB{Float32}(0.671,0.867,0.643),
				RGB{Float32}(0.4,0.761,0.647),
				RGB{Float32}(0.196,0.533,0.741),
				RGB{Float32}(0.369,0.31,0.635)]

const RdYlGn =   [RGB{Float32}(0.647,0.0,0.149),
				RGB{Float32}(0.843,0.188,0.153),
				RGB{Float32}(0.957,0.427,0.263),
				RGB{Float32}(0.992,0.682,0.38),
				RGB{Float32}(0.996,0.878,0.545),
				RGB{Float32}(1.0,1.0,0.749),
				RGB{Float32}(0.851,0.937,0.545),
				RGB{Float32}(0.651,0.851,0.416),
				RGB{Float32}(0.4,0.741,0.388),
				RGB{Float32}(0.102,0.596,0.314),
				RGB{Float32}(0.0,0.408,0.216)]

const RdBu =     [RGB{Float32}(0.404,0.0,0.122),
				RGB{Float32}(0.698,0.094,0.169),
				RGB{Float32}(0.839,0.376,0.302),
				RGB{Float32}(0.957,0.647,0.51),
				RGB{Float32}(0.992,0.859,0.78),
				RGB{Float32}(0.969,0.969,0.969),
				RGB{Float32}(0.82,0.898,0.941),
				RGB{Float32}(0.573,0.773,0.871),
				RGB{Float32}(0.263,0.576,0.765),
				RGB{Float32}(0.129,0.4,0.675),
				RGB{Float32}(0.02,0.188,0.38)]

const PiYG =     [RGB{Float32}(0.557,0.004,0.322),
				RGB{Float32}(0.773,0.106,0.49),
				RGB{Float32}(0.871,0.467,0.682),
				RGB{Float32}(0.945,0.714,0.855),
				RGB{Float32}(0.992,0.878,0.937),
				RGB{Float32}(0.969,0.969,0.969),
				RGB{Float32}(0.902,0.961,0.816),
				RGB{Float32}(0.722,0.882,0.525),
				RGB{Float32}(0.498,0.737,0.255),
				RGB{Float32}(0.302,0.573,0.129),
				RGB{Float32}(0.153,0.392,0.098)]

const YG = [RGB{Float32}(0.153,0.392,0.098),RGB{Float32}(0.153,0.392,0.098)]

const PRGn =     [RGB{Float32}(0.251,0.0,0.294),
				RGB{Float32}(0.463,0.165,0.514),
				RGB{Float32}(0.6,0.439,0.671),
				RGB{Float32}(0.761,0.647,0.812),
				RGB{Float32}(0.906,0.831,0.91),
				RGB{Float32}(0.969,0.969,0.969),
				RGB{Float32}(0.851,0.941,0.827),
				RGB{Float32}(0.651,0.859,0.627),
				RGB{Float32}(0.353,0.682,0.38),
				RGB{Float32}(0.106,0.471,0.216),
				RGB{Float32}(0.0,0.267,0.106)]

const RdYlBu =  [RGB{Float32}(0.647,0.0,0.149),
				RGB{Float32}(0.843,0.188,0.153),
				RGB{Float32}(0.957,0.427,0.263),
				RGB{Float32}(0.992,0.682,0.38),
				RGB{Float32}(0.996,0.878,0.565),
				RGB{Float32}(1.0,1.0,0.749),
				RGB{Float32}(0.878,0.953,0.973),
				RGB{Float32}(0.671,0.851,0.914),
				RGB{Float32}(0.455,0.678,0.82),
				RGB{Float32}(0.271,0.459,0.706),
				RGB{Float32}(0.192,0.212,0.584)]

const BrBG =     [RGB{Float32}(0.329,0.188,0.02),
				RGB{Float32}(0.549,0.318,0.039),
				RGB{Float32}(0.749,0.506,0.176),
				RGB{Float32}(0.875,0.761,0.49),
				RGB{Float32}(0.965,0.91,0.765),
				RGB{Float32}(0.961,0.961,0.961),
				RGB{Float32}(0.78,0.918,0.898),
				RGB{Float32}(0.502,0.804,0.757),
				RGB{Float32}(0.208,0.592,0.561),
				RGB{Float32}(0.004,0.4,0.369),
				RGB{Float32}(0.0,0.235,0.188)]

const RdGy =     [RGB{Float32}(0.404,0.0,0.122),
				RGB{Float32}(0.698,0.094,0.169),
				RGB{Float32}(0.839,0.376,0.302),
				RGB{Float32}(0.957,0.647,0.51),
				RGB{Float32}(0.992,0.859,0.78),
				RGB{Float32}(1.0,1.0,1.0),
				RGB{Float32}(0.878,0.878,0.878),
				RGB{Float32}(0.729,0.729,0.729),
				RGB{Float32}(0.529,0.529,0.529),
				RGB{Float32}(0.302,0.302,0.302),
				RGB{Float32}(0.102,0.102,0.102)]

const PuOr =     [RGB{Float32}(0.498,0.231,0.031),
				RGB{Float32}(0.702,0.345,0.024),
				RGB{Float32}(0.878,0.51,0.078),
				RGB{Float32}(0.992,0.722,0.388),
				RGB{Float32}(0.996,0.878,0.714),
				RGB{Float32}(0.969,0.969,0.969),
				RGB{Float32}(0.847,0.855,0.922),
				RGB{Float32}(0.698,0.671,0.824),
				RGB{Float32}(0.502,0.451,0.675),
				RGB{Float32}(0.329,0.153,0.533),
				RGB{Float32}(0.176,0.0,0.294)]

const Paired =   [RGB{Float32}(0.651,0.808,0.89),
				RGB{Float32}(0.122,0.471,0.706),
				RGB{Float32}(0.698,0.875,0.541),
				RGB{Float32}(0.2,0.627,0.173),
				RGB{Float32}(0.984,0.604,0.6),
				RGB{Float32}(0.89,0.102,0.11),
				RGB{Float32}(0.992,0.749,0.435),
				RGB{Float32}(1.0,0.498,0.0),
				RGB{Float32}(0.792,0.698,0.839),
				RGB{Float32}(0.416,0.239,0.604),
				RGB{Float32}(1.0,1.0,0.6)]

const Purples =  [RGB{Float32}(0.988,0.984,0.992),
				RGB{Float32}(0.937,0.929,0.961),
				RGB{Float32}(0.855,0.855,0.922),
				RGB{Float32}(0.737,0.741,0.863),
				RGB{Float32}(0.62,0.604,0.784),
				RGB{Float32}(0.502,0.49,0.729),
				RGB{Float32}(0.416,0.318,0.639),
				RGB{Float32}(0.329,0.153,0.561),
				RGB{Float32}(0.247,0.0,0.49)]

const BkOr =  [RGB{Float32}(0.0,0.0,0.0),
				RGB{Float32}(174.0/255.0,	112.0/255.0,	0.)]

const BuYlRd =  [RGB{Float32}(0.02,0.188,0.38),
				RGB{Float32}(204.0/255.,204.0/255.,0.),
				RGB{Float32}(0.404,0.0,0.122)]

const BlPu = [	RGB{Float32}(0.0,0.0,0.0),
				RGB{Float32}(0.62,0.604,0.784)]

const GrBl = [RGB{Float32}(0.6,0.6,0.6),
		RGB{Float32}(0.,0.,0.)]

const distgclrs = [RGB{Float32}(0., 0., 128.0/255.),
				RGB{Float32}(0., 130.0/255., 200.0/255.),
				RGB{Float32}(70.0/255., 240.0/255., 240.0/255.),
				RGB{Float32}(60.0/255., 180.0/255., 75.0/255.),
				RGB{Float32}(210.0/255., 245.0/255., 60.0/255.),
				RGB{Float32}(255.0/255., 225.0/255., 25.0/255.),
				RGB{Float32}(245.0/255., 130.0/255., 48.0/255.),
				RGB{Float32}(230.0/255., 25.0/255., 75.0/255.),
				RGB{Float32}(128.0/255., 0., 0.)]

const bone = 	[RGB{Float32}(0.    ,    0.	  ,    0.0052),
    			RGB{Float32}(0.0972,    0.0972,    0.1389),
    			RGB{Float32}(0.1944,    0.1944,    0.2726),
    			RGB{Float32}(0.2917,    0.2917,    0.4062),
    			RGB{Float32}(0.3889,    0.4149,    0.5139),
    			RGB{Float32}(0.4861,    0.5486,    0.6111),
    			RGB{Float32}(0.5833,    0.6823,    0.7083),
    			RGB{Float32}(0.6962,    0.8056,    0.8056),
    			RGB{Float32}(0.8481,    0.9028,    0.9028),
    			RGB{Float32}(1.0000,    1.0000,    1.0000)]


text:
using SimpleTraits

# Type hierarchy for cheap implications
abstract type MeasureTrait{X} <: SimpleTraits.Trait end
abstract type IsRepresentative{X} <: MeasureTrait{X} end
abstract type IsPrimitive{X} <: IsRepresentative{X} end

abstract type IsScalar{X} <: MeasureTrait{X} end

"""
    @scalar T

Declare that every measure of type `T` are scalar.
"""
macro scalar(T)
    MT = MeasureTheory
    ST = MeasureTheory.SimpleTraits

    # Really the same as `@traitimpl IsScalar{T}`
    # But maybe easier for beginners
    quote
        function $ST.trait(::$ST.Type{$MT.IsScalar{S}}) where {S <: $T}
            return $MT.IsScalar{S}
        end
    end
end

"""
    isscalar(::AbstractMeasure)
    isscalar(::Type{M}) where {M<:AbstractMeasure}

Tests whether a given measure has been declared scalar. See `@scalar`.
"""
function isscalar end

isscalar(::M) where {M} = isscalar(M)
isscalar(::Type{M}) where {M} = istrait(IsScalar{M})

"""
    isprimitive(::AbstractMeasure)
    isprimitive(::Type{M}) where {M<:AbstractMeasure}

Most measures are defined in terms of other measures, for example using a
density or a pushforward. Those that are not are considered (in this library,
it's not a general measure theory thing) to be _primitive_. The canonical
example of a primitive measure is `Lebesgue(X)` for some `X`.
"""
function isprimitive end

isprimitive(::M) where {M} = isprimitive(M)
isprimitive(::Type{M}) where {M} = istrait(IsPrimitive{T})

"""
    @primitive T

Declare that every measure of type `T` should be considered "primitive". A
measure is primitive is it is not defined in terms of another measure. Common
examples are Lebesgue and Counting measures.

Note that this is not a general measure-theoretic term, but is specific to the
MeasureTheory.jl implementation.
"""
macro primitive(T)
    MT = MeasureTheory
    ST = MeasureTheory.SimpleTraits

    # Really the same as `@traitimpl IsPrimitive{T}`
    # But maybe easier for beginners
    quote
        function $ST.trait(::$ST.Type{$MT.IsPrimitive{S}}) where {S <: $T}
            return $MT.IsPrimitive{S}
        end
    end
end

export basemeasure

"""
    basemeasure(μ)

Many measures are defined in terms of a logdensity relative to some base
measure. This makes it important to be able to find that base measure.

For measures not defined in this way, we'll typically have `basemeasure(μ) == μ`.
"""
function basemeasure end


text:
@testset "Localhost" begin
    wmic_command = "wmic path win32_computersystemproduct get uuid"

    windows_patch = @patch Sys.iswindows() = return true

    wmic_patch = @patch read(command::Cmd, ::Type{String}) =
        return "UUID\nEC2D1284-E32E-FB5E-20E4-F43F6E01CA7A"

    @testset "EC2 - Windows" begin
        apply([windows_patch, wmic_patch]) do
            @test localhost_is_ec2()
        end
    end
end

if Sys.iswindows()
    @testset "Windows - Is EC2" begin
        @test localhost_is_ec2()
    end
end


text:
function GBlackScholes(CallPutFlag, S, X, T, r, b, v)

  # b = r, give the BS stock option model
  # b = r - q, gives Merton 73, stock option with dividend
  # b = 0, give black 76
  # b = 0 && r = 0, gives Asay for margined futures options
  # b = r- r_f, gives garman kohlhagen currency options model

  d1 = (log(S / X) + (b + v^2 / 2) * T) / (v * sqrt(T))
  d2 = d1 - v * sqrt(T)

  if CallPutFlag == "c"
    GBS = S * exp((b - r) * T) * CND(d1) - X * exp(-r * T) * CND(d2)
  elseif CallPutFlag == "p"
    GBS = X * exp(-r * T) * CND(-d2) - S * exp((b - r) * T) * CND(-d1)
  end

end

function KirkSpreadOption(CallPutFlag, S1, S2, X, T, r, b1, b2, v1, v2, cor)

        S = (S1 * exp((b1 - r) * T)) / (S2 * exp((b2 - r) * T) + X * exp((-r) * T))
        F = (S2 * exp((b2 - r) * T)) / (S2 * exp((b2 - r) * T) + X * exp((-r) * T))
        sig = sqrt(v1^2 + (v2 * F)^2 - 2 * cor * v1 * v2 * F)

        d1 = (log(S) + (sig^2 / 2) * T) / (sig * sqrt(T))
        d2 = d1 - sig * sqrt(T)

        if CallPutFlag == "c"

                return (S2 * exp((b2 - r) * T) + X * exp((-r) * T)) * (S * CND(d1) - CND(d2))

        elseif CallPutFlag == "p"

                return (S2 * exp((b2 - r) * T) + X * exp((-r) * T)) * (CND(-d2) - S * CND(-d1))

        else

                println(STDERR, "Falsche Eingabe! Abbruch\n")
                return 0.0 #eventuell nochmal checken ob das sinn macht

        end

end

function GeoAveRate(CallPutFlag, S, X, T, r, b, v)

        sigAdj = v / sqrt(3)
        bAdj = (b - v^2 / 6) / 2

        d1 = (log(S / X) + (bAdj + (sigAdj^2 / 2)) * T) / (sigAdj * sqrt(T))
        d2 = d1 - sigAdj * sqrt(T)

        if CallPutFlag == "c"

                return S * exp((bAdj - r) * T) * CND(d1) - X * exp((-r) * T) * CND(d2)

        elseif CallPutFlag == "p"

                return X * exp((-r) * T) * CND((-d2)) - S * exp((bAdj - r) * T) * CND((-d1))

        else

                println(STDERR, "Falsche Eingabe! Abbruch\n")
                return 0.0 #eventuell nochmal checken ob das sinn macht
        end

end


text:
const range_h = """
struct Range {
    const int start;
    const int step;
public:
    __device__ int operator()(const int idx) { return start + step*idx; }
};

template<int N>
struct Ranges {
    Range data[N];
public:
    __device__ Range& operator[](const int idx) { return data[idx]; }

    __device__ void convert(int *subs) {
        for (int i = 0; i < N; i++) {
            subs[i] = data[i](subs[i]);
        }
    }
};
"""


text:
using Dates
const SRC_PATH = "../../../src"

include(joinpath(SRC_PATH, "MLPClassifier/MLPClassifier.jl"))
include(joinpath(SRC_PATH, "datasets/dataset_reader.jl"))
include(joinpath(SRC_PATH, "datasets/Wine-Classification/dataset3.jl"))
include(joinpath(SRC_PATH, "FileInterface/fileInterface.jl"))
#include("../MLP_test_functions.jl")
# include(joinpath(SRC_PATH, "MutualInformation/information-plane.jl"))

t0 = time()
#topology = [22,10]
#topology = [11,22]
#topology = [22,10,15,11]
topology = [16]

Random.seed!(1234)
net1 = MLPClassifier(dataset3_nfeatures, dataset3_labels, topology)
net2 = deepcopy(net1)

path = joinpath(SRC_PATH, "datasets/Wine-Classification")
dt = get_dataset3_setup(path, normalise=true, separate=true)
dt_b = get_dataset3_setup(path, normalise=true, separate=true, balance=true)

epochs = 5000

println("Topology: ", topology)
println("Epochs: ", epochs)

#Random.seed!(1234)
println("-> Not balanced")
println()
println("--> Results before training")
show_results(dt, net1)
fit!(net1, dt.input_train, dt.output_train; max_epochs = epochs)
println()
println("--> Results after training")
show_results(dt, net1)
println()
t1 = time()

#Random.seed!(1234)
println("-> Balanced")
println()
println("--> Results before training")
show_results(dt_b, net2)
fit!(net2, dt_b.input_train, dt_b.output_train; max_epochs = epochs)
println()
println("--> Results after training")
show_results(dt_b, net2)

println()

t2 = time()
println("Training time (not balanced):", t1 - t0," seconds.")
println("Training time (balanced):", t2 - t1," seconds.")

current_date = Dates.now()
time_str = Dates.format(current_date, "HH:MM:SS")
date_str = Dates.format(current_date, "dd-mm-yyyy")

#filename = "MLP-nets-dt3.jld2"
filename = string("MLP-nets-dt3-", date_str, "-", time_str, ".jld2")
addToFile(filename, "Net-notbalanced", net1)
addToFile(filename, "Net-balanced", net2)

#=
data = zip(dt.input_train, dt.output_train)
epochs = 10000
number_nets = 20
step_len = 1
nbins = 5
plot_IP_average(dt.n_features, dt.labels, topology, number_nets, data, epochs, step_len, "IP-dt1-MLP.png")
points_IP = calculate_IP_average_points(dt.n_features, dt.labels, topology, number_nets, data, epochs, step_len, nbins)
generate_plot_IP_average(points_IP, number_nets, epochs, step_len, "IP-dt1-MLP-$(nbins)bins.png")
=#
# Criar uma função para salvar as redes em um arquivo

#=
path = joinpath(SRC_PATH, "datasets/Wine-Classification")
dt = get_dataset3_setup(path, normalise=true, separate=true, balance=true)

topology = [16]
data = zip(dt.input_train, dt.output_train)
data_test = zip(dt.input_test, dt.output_test)
epochs = 5000
number_nets = 10
step_len = 250
nbins = 10

points_IP = calculate_IP_average_points(dt.n_features, dt.labels, topology, number_nets, data, epochs, step_len, nbins)
generate_plot_IP_average(points_IP, number_nets, epochs, step_len, "IP-dt1-MLP-$(nbins)bins.png")

MLP_error_rate(dt.n_features, dt.labels, topology, number_nets, data, data_test, epochs)
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
module PrototypesCriticisms

using Clustering
using Distances
using KernelFunctions
using LinearAlgebra
using Statistics

export criticisms,
        prototypes,
        sqmmd, mmd²,
        witness

"""
    sqmmd(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}, k::Kernel)
    mmd²(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}, k::Kernel)

Compute the squared maximum mean discrepancy between `X` and `Y` using the kernel function `k`.

`X` and `Y` are expected to store observations in columns.
"""
function sqmmd(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}, k::Kernel)
    XX = kernelmatrix(k, X, obsdim=2)
    XY = kernelmatrix(k, X, Y, obsdim=2)
    YY = kernelmatrix(k, Y, obsdim=2)
    return sqmmd(XX, XY, YY)
end

"""
    sqmmd(XX::AbstractMatrix{<:Real}, XY::AbstractMatrix{<:Real}, YY::AbstractMatrix{<:Real})
    mmd²(XX::AbstractMatrix{<:Real}, XY::AbstractMatrix{<:Real}, YY::AbstractMatrix{<:Real})

Compute the squared maximum mean discrepancy using the kernel matrices `XX`, `XY`, and `YY`.

`XX` is the kernel matrix of the matrices `X` and `X` etc.
"""
sqmmd(XX::AbstractMatrix{<:Real}, XY::AbstractMatrix{<:Real}, YY::AbstractMatrix{<:Real}) = mean(XX) - 2 * mean(XY) + mean(YY)

# Alias
const mmd² = sqmmd

"""
    prototypes(X::AbstractMatrix{<:Real}, k::Kernel, n::Int)

Return the indices of `n` prototypes in `X` using the kernel function `k`.

`X` is expected to store observations in columns.
"""
function prototypes(X::AbstractMatrix{<:Real}, k::Kernel, n::Int)
    K = kernelmatrix(k, X, obsdim=2)
    return prototypes(K, n)
end

"""
    prototypes(K::AbstractMatrix{<:Real}, n::Int)

Return the indices of `n` prototypes using the kernel matrix `K`.
"""
function prototypes(K::AbstractMatrix{<:Real}, n::Int)
    doubledkernelmeans = 2 * mean(K, dims=1)
    initialcandidates = 1:size(K, 2)
    protoids = Int[]
    while length(protoids) < n
        candidates = setdiff(initialcandidates, protoids)
        avgproximities1 = vec(doubledkernelmeans[candidates])
        if length(protoids) > 0
            proximities2 = vec(2 * sum(K[protoids, candidates], dims=1))
            avgproximities2 = proximities2 / (length(protoids) + 1)
            diffs = avgproximities1 - avgproximities2
        else
            diffs = avgproximities1
        end
        protoid = candidates[argmax(diffs)]
        push!(protoids, protoid)
    end
    return protoids
end

"""
    prototypes(c::KmedoidsResult, n::Int=1)

Return the indices of the `n` prototypes for every cluster of the k-medoids clustering `c`.
"""
function prototypes(c::KmedoidsResult, n::Int=1)
    n == 1 && return [[i] for i in c.medoids]
    return _instances(nclusters(c), assignments(c), c.costs, n)
end

"""
    prototypes(c::KmeansResult, n::Int=1)

Return the indices of the `n` prototypes for every cluster of the k-means clustering `c`.
"""
prototypes(c::KmeansResult, n::Int=1) = _instances(nclusters(c), assignments(c), c.costs, n)

"""
    prototypes(c::FuzzyCMeansResult, n::Int=1)

Return the indices of the `n` prototypes for every cluster of the fuzzy c-means clustering `c`.
"""
prototypes(c::FuzzyCMeansResult, n::Int=1) = _instances(c.weights, n, true)

"""
    prototypes(c::AffinityPropResult, X::AbstractMatrix{<:Real}, n::Int=1; dist::Metric=Euclidean())

Return the indices of the `n` prototypes for every cluster of the affinity propagation clustering `c`.
To find `n>1` prototypes, the observations `X` are used for comparing the distances between observations and cluster centers.

# Keyword arguments
- `dist`: the distance metric used for computing the distances.
"""
function prototypes(c::AffinityPropResult, X::AbstractMatrix{<:Real}, n::Int=1; dist::Metric=Euclidean())
    n == 1 && return [[i] for i in c.exemplars]
    D = pairwise(dist, X, dims=2)
    instances = Vector{Vector{Int}}()
    for i = 1:nclusters(c)
        v = view(D, assignments(c) .== i, c.exemplars[i])
        permutation = partialsortperm(v, 1:n)
        parentinstances = parentindices(v)[1]
        push!(instances, parentinstances[permutation])
    end
    return instances
end

"""
    witness(z::AbstractVector{<:Real}, X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}, k::Kernel)

Compute the witness function of `X` and `Y` at `z` using the kernel function `k`.

`X` and `Y` are expected to store observations in columns, and `z` is a single observation.
"""
witness(z::AbstractVector{<:Real}, X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}, k::Kernel) = mean(map(x -> k(z, x), eachcol(X))) - mean(map(y -> k(z, y), eachcol(Y)))

"""
    criticisms(X::AbstractMatrix{<:Real}, k::Kernel, protoids::AbstractVector{Int}, n::Int)

Return the indices of `n` criticisms in `X` using the prototype indices `protoids` and kernel function `k`.

`X` is expected to store observations in columns.
"""
function criticisms(X::AbstractMatrix{<:Real}, k::Kernel, protoids::AbstractVector{Int}, n::Int)
    K = kernelmatrix(k, X, obsdim=2)
    return criticisms(K, protoids, n)
end

"""
    criticisms(K::AbstractMatrix{<:Real}, protoids::AbstractVector{Int}, n::Int)

Return the indices of `n` criticisms using the prototype indices `indices` and the kernel matrix `K`.
"""
function criticisms(K::AbstractMatrix{<:Real}, protoids::AbstractVector{Int}, n::Int)
    kernelmeans = mean(K, dims=1)
    initialcandidates = setdiff(1:size(K, 2), protoids)
    critids = Int[]
    while length(critids) < n
        candidates = setdiff(initialcandidates, critids)
        avgproximities1 = kernelmeans[candidates]
        avgproximities2 = mean(K[protoids, candidates], dims=1)
        absws = abs.(vec(avgproximities1) - vec(avgproximities2)) .+ logdet(K[critids, critids])
        critid = candidates[argmax(absws)]
        push!(critids, critid)
    end
    return critids
end

"""
    criticisms(c::KmedoidsResult, n::Int=1)

Return the indices of the `n` criticisms for every cluster of the k-medoids clustering `c`.
"""
criticisms(c::KmedoidsResult, n::Int=1) = _instances(nclusters(c), assignments(c), c.costs, n, true)

"""
    criticisms(c::KmeansResult, n::Int=1)

Return the indices of the `n` criticisms for every cluster of the k-means clustering `c`.
"""
criticisms(c::KmeansResult, n::Int=1) = _instances(nclusters(c), assignments(c), c.costs, n, true)

"""
    criticisms(c::FuzzyCMeansResult, n::Int=1)

Return the indices of the `n` criticisms for every cluster of the fuzzy c-means clustering `c`.
"""
criticisms(c::FuzzyCMeansResult, n::Int=1) = _instances(c.weights, n)

"""
    criticisms(c::AffinityPropResult, X::AbstractMatrix{<:Real}, n::Int=1; dist::Metric=Euclidean())

Return the indices of the `n` criticisms for every cluster of the affinity propagation clustering `c`.
To find `n>1` criticisms, the observations `X` are used for comparing the distances between observations and cluster centers.

# Keyword arguments
- `dist`: the distance metric used for computing the distances.
"""
function criticisms(c::AffinityPropResult, X::AbstractMatrix{<:Real}, n::Int=1; dist::Metric=Euclidean())
    D = pairwise(dist, X, dims=2)
    instances = Vector{Vector{Int}}()
    for i = 1:nclusters(c)
        v = view(D, assignments(c) .== i, c.exemplars[i])
        permutation = partialsortperm(v, 1:n, rev=true)
        parentinstances = parentindices(v)[1]
        push!(instances, parentinstances[permutation])
    end
    return instances
end

# Return instances via their assignment costs
function _instances(k::Int, assignments::Vector{Int}, costs::Vector{<:Real}, n::Int, rev::Bool=false)
    clustercosts = [[] for _ in 1:k]
    for i = 1:length(assignments)
        push!(clustercosts[assignments[i]], (costs[i], i))
    end
    return map.(i -> i[2], partialsort!.(clustercosts, [1:n], rev=rev))
end

# Return instances via their assignment weights
function _instances(W::AbstractMatrix{<:Real}, n::Int=1, rev::Bool=false)
    instances = Vector{Vector{Int}}()
    hardassignments = vec(map(i -> i[2], argmax(W, dims=2)))
    for i = 1:size(W, 2)
        v = view(W, hardassignments .== i, i)
        permutation = partialsortperm(v, 1:n, rev=rev)
        parentinstances = parentindices(v)[1]
        push!(instances, parentinstances[permutation])
    end
    return instances
end

end # module


text:
using Test
using Psychro

include("test_utilities.jl")
include("test_hyland83a.jl")
include("test_hyland83.jl")
include("test_user.jl")
include("test_units.jl")