text:
export Solver, solver_startsolutions, solver, solve

###########
## STATS ##
###########
mutable struct SolveStats
    regular::Int
    regular_real::Int
    singular::Int
    singular_real::Int
end
SolveStats() = SolveStats(0, 0, 0, 0)

function init!(SS::SolveStats)
    SS.regular = SS.regular_real = SS.singular = SS.singular_real = 0
    SS
end

function update!(stats::SolveStats, R::PathResult)
    is_success(R) || return stats

    if is_singular(R)
        stats.singular_real += is_real(R)
        stats.singular += 1
    else
        stats.regular_real += is_real(R)
        stats.regular += 1
    end
    stats
end

function consolidated_stats(stats::Vector{SolveStats})
    (
        regular = sum(s.regular for s in stats),
        regular_real = sum(s.regular_real for s in stats),
        singular = sum(s.singular for s in stats),
        singular_real = sum(s.singular_real for s in stats),
    )
end

######################
## PathJumpingCheck ##
######################
function path_jumping_candidates(
    results::Vector{Union{Nothing,PathResult{V}}},
    tol::Float64,
) where {V}
    index_map = Int[]
    S = Vector{V}()
    for (i, r) in enumerate(results)
        if r !== nothing
            s = r.intermediate_solution
            if s !== nothing
                push!(S, s)
                push!(index_map, i)
            end
        end
    end

    clusters = multiplicities(S; tol = tol)
    indices = Int[]
    for cluster in clusters, i in cluster
        push!(indices, index_map[i])
    end
    indices
end

############
## SOLVER ##
############

"""
    Solver(pathtracker::PathTracker)

A `Solver` is a wrapper around a given `pathtracker` to track multiple paths.
It provides on top of the given `pathtracker` parallelization and an optional path jumping check.
To construct a  `Solver` it is convenient to use the [`solver`](@ref) or
[`solver_startsolutions`](@ref) functions. Given a solver one can use the [`solve`](@ref)
function to solve a system. This struct is constructed for any call to `solve` unless
already explicitly provided.

## Example

Assume we want to solve a polynomial system repeatedly for many different values of the
parameters `p`.
The following example shows how to use a `Solver` to avoid some computational overhead
compared to naively calling [`solve`](@ref).

```julia
@polyvar x y z p[1:3]
F = [
    x + 3 + 2y + 2 * y^2 - p[1],
    (x - 2 + 5y) * z + 4 - p[2] * z,
    (x + 2 + 4y) * z + 5 - p[3] * z,
]
q = randn(ComplexF64, 3)
S = solutions(solve(subs(F, p => q)))
# create some fake parameter values
params = [randn(3) for _ = 1:1000]
# create a `Solver` to reuse for the path tracking
F_solver = solver(F; parameters = p, generic_parameters = q)
# solve the system F for all paramaters p in params
params_solutions = map(params) do p
    solutions(solve(F_solver, S; target_parameters = p))
end
```
"""
struct Solver{PT<:AbstractPathTracker}
    trackers::Vector{PT}
    stats::Vector{SolveStats}
end
Solver(tracker::AbstractPathTracker) = Solver([tracker], [SolveStats()])
function Solver(prob::AbstractProblem, start_solutions; kwargs...)
    Solver(construct_tracker(prob, start_solutions; kwargs...))
end

"""
    solver_startsolutions(args...; kwargs...)

Create a [`Solver`](@ref) and start solutions. Takes almost the same arguments as [`solve`](@ref).
"""
function solver_startsolutions(args...; kwargs...)
    supported, rest = splitkwargs(kwargs, problem_startsolutions_supported_keywords)
    prob, start_solutions = problem_startsolutions(args...; supported...)
    Solver(prob, start_solutions; rest...), start_solutions
end

"""
    solver(args...; kwargs...)

Create a [`Solver`](@ref). Takes almost the same arguments as [`solve`](@ref).
"""
solver(args...; kwargs...) = first(solver_startsolutions(args...; kwargs...))

result_type(solver::Solver) = Result{typeof(solution(solver.trackers[1]))}

min_accuracy(T::PathTracker) = T.options.min_accuracy
min_accuracy(T::OverdeterminedTracker) = min_accuracy(T.tracker)
min_accuracy(T::PolyhedralTracker) = min_accuracy(T.generic_tracker)

accuracy(T::PathTracker) = T.default_ct_options.accuracy
accuracy(T::OverdeterminedTracker) = accuracy(T.tracker)
accuracy(T::PolyhedralTracker) = accuracy(T.generic_tracker)

max_corrector_iters(T::PathTracker) = T.default_ct_options.max_corrector_iters
max_corrector_iters(T::OverdeterminedTracker) = max_corrector_iters(T.tracker)
max_corrector_iters(T::PolyhedralTracker) = max_corrector_iters(T.generic_tracker)

const solve_supported_keywords = [
    :path_result_details,
    :path_jumping_check,
    :save_all_paths,
    :show_progress,
    :stop_early_cb,
    :threading,
            # deprecated
    :report_progress,
]

"""
    solve(args...; options...)::Result

The solve function takes many different arguments and options depending on your specific situation,
but in the end it always returns a [`Result`](@ref) containing the result of the computations.
Depending on the prodived arguments different kind of homotopies are constructed. In particular
it is possible to construct the following homotopies:

* Total degree homotopy
* Polyhedral homotopy
* Parameter homotopy
* Multi-homogeneous homotopy
* Start target homotopy

If the input is a *homogeneous* polynomial system, solutions in projective space are computed.
Otherwise affine solutions are computed.

## Options

The `solve` routines takes many different options. In particular all options to
[`CoreTracker`](@ref) and [`PathTracker`](@ref) are allowed.
Additionally the following options are allowed:

* `affine_tracking::Bool=true`: Indicate whether path tracking should happen in affine space.
* `early_stop_cb`: Here it is possible to provide a function (or any callable struct) which
  accepts a `PathResult` `r` as input and returns a `Bool`. If `early_stop_cb(r)` is `true`
  then no further paths are tracked and the computation is finished. This is only called
  for successfull paths unless `save_all_paths` is `true`. This is for example useful
  if you only want to compute one solution of a polynomial system.
  For this `early_stop_cb = _ -> true` would be sufficient.
* `homotopy::AbstractHomotopy`: A constructor to construct a [`AbstractHomotopy`](@ref) for
  the totaldegree and start target homotopy. The default is [`StraightLineHomotopy`](@ref).
  The constructor is called with `homotopy(start, target)` where `start` and `target` are
  homogeneous [`AbstractSystem`](@ref)s.
* `homvar::Union{Int,MultivariatePolynomials.AbstractVariable}`: This considers the
  *homogeneous* system `F` as an affine system which was homogenized by `homvar`.
  If `F` is an `AbstractSystem` `homvar` is the index (i.e. `Int`) of the homogenization
  variable. If `F` is an `AbstractVariables` (e.g. created by `@polyvar x`)
  `homvar` is the actual variable used in the system `F`.
* `path_jumping_check::Bool=true`: Enable a check whether one of
  the paths jumped to another one.
* `path_result_details=:default`: The amount of information computed in each path result.
  Possible values are `:minimal` (minimal details), `:default` (default) and `:extensive`.
* `projective_tracking::Bool=false`: Indicate whether path tracking should happen in
  projective space. The flag `affine_tracking` is dominant.
* `seed`: The random seed used during the computations. The seed is also reported in the
  result. For a given random seed the result is always identical.
* `show_progress` (default `true`): Whether a progress bar should be printed to report the
  progress of the current computation.
* `system::AbstractSystem`: A constructor to assemble a [`AbstractSystem`](@ref).
  The default is [`SPSystem`](@ref). This constructor is only applied to the input of `solve`.
  The constructor is called with `system(polynomials, variables)` where `polynomials` is a
  vector of `MultivariatePolynomials.AbstractPolynomial`s and `variables` determines the
  variable ordering. If you experience significant compilation times,
  consider to change system to `FPSystem`.
* `system_scaling` (default `:equations`) Whether to apply an automatic scaling of
  the equations (:equations), of the equations and variables (`:equations_and_variables`) or
  no scaling at all (`nothing`).
* `threading` (default `true`): Enable or disable multi-threading. The number of threads used
  is controlled by the environment variable `JULIA_NUM_THREADS`.
* `variable_ordering`: Provide a custom ordering of the variables.

# Examples

## Total Degree Homotopy

    solve(F; options...)

Solve the system `F` using a start system computed from the degrees of the entries of `F`.
The number of paths to track is equal to the total degree `d₁⋯dⱼ`, where `dᵢ` is
the degree of the `i`th entry of `F`. `F` can be

- `Vector{<:MultivariatePolynomials.AbstractPolynomial}` (e.g. constructed by using the
  exported `@polyvar`)
- A composition of polynomial systems constructed by [`compose`](@ref)
- Any [`AbstractSystem`](@ref)

We can solve the system ``F(x,y) = (x^2+y^2+1, 2x+3y-1)`` in the following way:

```julia
julia> @polyvar x y;
julia> solve([x^2+y^2+1, 2x+3y-1])
Result with 2 solutions
==================================
• 2 non-singular solutions (0 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 661766
```

## Polyhedral Homotopy

    solve(F; start_system = :polyhedral, only_torus=false, options...)

Solve the system `F` in two steps: first solve a generic system computed from
the Newton Polytopes of the entries of `F`, then perform a coefficient
homotopy towards `F`. The number of paths to track is equal to the mixed volume of the
Newton Polytopes of the entries of `F`. The mixed volume is at most the total degree of `F`.
`F` can be

- `Vector{<:MultivariatePolynomials.AbstractPolynomial}` (e.g. constructed by `@polyvar`)
- A composition of polynomial systems constructed by [`compose`](@ref). Note that the
  composition will not preserved.

If `only_torus == true` then only solutions in the algebraic torus
``(ℂ\\setminus \\{0\\})^n`` will be computed.

We can solve the system ``F(x,y) = (x^2+y^2+1, 2x+3y-1)`` in the following way:
```julia
julia> @polyvar x y;
julia> solve([x^2+y^2+1, 2x+3y-1]; start_system = :polyhedral)
Result with 2 solutions
==================================
• 2 non-singular solutions (0 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 222880
```

## Parameter Homotopy

    solve(F, startsolutions;
        parameters,
        start_parameters,
        target_parameters,
        start_gamma = nothing,
        target_gamma = nothing,
    )

Solve the parameter homotopy
```math
H(x, t) = F(x, \\frac{tγ₁p₁+(1-t)γ₀p₀}{tγ₁+(1-t)γ₀}),
```
where ``p₁`` (=`start_parameters`) and ``p₀`` (=`target_parameters`) are vectors of
parameter values for ``F`` and ``γ₁`` (=`start_gamma`) and ``γ₀`` (=`target_gamma`)
    are complex numbers.
If `start_parameters` or `target_parameters` is `nothing`, it is assumed that `γ₀=γ₁=1`.
The input `parameters` specifies the variables of `F` which should be considered as parameters.
Necessarily we have `length(parameters) == length(p₁) == length(p₀)`.

    solve(F, startsolutions; parameters, p₁, p₀, γ₁=nothing, γ₀=nothing)

This is a unicode variant where `γ₁=start_parameters`, `γ₀=target_parameters`,
    `γ₁=start_gamma`, γ₀=`target_gamma`.

We want to solve a parameter homotopy ``H(x,t) := F(x; t[1, 0]+(1-t)[2, 4])`` where
```math
F(x; a) := (x₁^2-a₁, x₁x₂-a₁+a₂)
```
and let's say we are only intersted in tracking of ``[1,1]``.
This can be accomplished as follows
```julia
@polyvar x[1:2] a[1:2]
F = [x[1]^2-a[1], x[1]*x[2]-a[1]+a[2]]
startsolutions = [[1, 1]]
p₁ = [1, 0]
p₀ = [3im, 0.5+2im]
solve(F, startsolutions; parameters=a, start_parameters=p₁, target_parameters=p₀)
# If you like unicode this is also possible
solve(F, startsolutions; parameters=a, p₁=p₁, p₀=p₀)
```

## Start Target Homotopy

    solve(G, F, start_solutions; options...)

This constructs the homotopy ``H(x,t) = tG(x)+(1-t)F(x)`` to compute solutions of the
system `F`. `start_solutions` is a list of solutions of `G` which are tracked to solutions
of `F`.

```julia
@polyvar x y
G = [x^2-1,y-1]
F = [x^2+y^2+z^2, 2x+3y-z]
solve(G, F, [[1, 1], [-1, 1]])
```

## Homogeneous Systems

If `F` consists of homogeneous equations, we return results in projective space

```julia
julia> @polyvar x y z;
julia> solve([x^2+y^2+z^2, 2x+3y-z])
Result{PVector{Complex{Float64},1}} with 2 solutions
====================================================
• 2 non-singular solutions (0 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 490575
```

It your polynomial system is not homogeneous, you can homogenize it as follows
```julia
@polyvar x y
g = [x^2+y^2+1, 2x+3y-1]
f = homogenize(g)
```
It is also possible to specify the homogenizing variable.
```julia
@polyvar x y z
g = [x^2+y^2+1, 2x+3y-1]
f = homogenize(g, z)
```

If your polynomial system is already homogeneous, but you would like to consider it as an affine system
you can do
```julia
@polyvar x y z
solve([x^2+y^2+z^2, 2x+3y-z], homvar=z)
```
This yields the same result as `solve([x^2+y^2+1, 2x+3y-1])`.

## Multi-homogeneous Systems

By exploiting the multi-homogeneous structure of a polynomial system it is possible
to decrease the number of paths necessary to track.

```julia
@polyvar x y
# Use variable groups to only track 2 paths instead of 4
solve([x*y - 6, x^2 - 5], variable_groups=[(x,), (y,)])
```
To check whether a certain variable grouping is beneficial you can use the [`bezout_number`](@ref)
function.

## Abstract Homotopy

    solve(H::AbstractHomotopy, start_solutions; options...)

Solve the homotopy `H` by tracking the each solution of
``H(⋅, t)`` (as provided by `start_solutions`) from ``t=1`` to ``t=0``.

## Many parameters

    solve(
        F, start_solutions;
        start_parameters = ...,
        target_parameters = ..., # multiple parameters
        options...,
        many_parameters_options...
    )

Solve the system paramterized system `F` for many different target parameters.
If no further options are specified (resp. only those which are explained above)
then the result of this funciton is similar to
```
map(params) do p
    res = solve(F, S; start_parameters = start_params, target_parameters = p, options...)
    (res, p)
end
```

Note if `threading = true` then the parallelization
is done on a parameter by parameter basis (and not as above on a path by path basis).
Note that order of the results is not guaranteed.
The `many_parameters_options` allow to apply transformations on the results

### Many parameter options

* `flatten::Bool`: Flatten the output of `transform_result`. This is useful for example if
  `transform_result` returns a vector of solutions, and you only want a single vector of
  solutions as the result (instead of a vector of vector of solutions).
* `transform_parameters::Function`: Transform a parameters values `p` before passing it to
  `target_parameters = ...`. By default this `identity`.
* `transform_result::Function`: A function taking two arguments, the `result` and the
  parameters `p`. By default this returns the tuple `(result, p)`.

### Many parameter examples

We start with some setup. We intersect a circle with a many different lines.

```julia
@polyvar x y
f = x^2 + y^2 - 1

@polyvar a b c
l = a * x + b * y + c
F = [f, l]

# Compute start solutions S₀ for given start parameters p₀
p₀ = randn(ComplexF64, 3)
S₀ = solutions(solve(subs(F, [a, b, c] => p₀)))
# The parameters we are intersted in
params = [rand(3) for i = 1:100]
```

Here are the some examples of the different options.
```julia-repl
julia> result1 = solve(
           F,
           S₀;
           parameters = [a, b, c],
           start_parameters = p₀,
           target_parameters = params,
       );

julia> typeof(result1)
Array{Tuple{Result{Array{Complex{Float64},1}},Array{Float64,1}},1}

julia> result1[1]
(Result{Array{Complex{Float64},1}} with 2 solutions
==================================================
• 2 non-singular solutions (2 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 814544
, [0.6323567622396351, 0.15381245892007867, 0.44070720179305134])

julia> # Only keep real solutions
       result2 = solve(
           F,
           S₀;
           parameters = [a, b, c],
           start_parameters = p₀,
           target_parameters = params,
           transform_result = (r, p) -> real_solutions(r),
       );

julia> typeof(result2)
Array{Array{Array{Float64,1},1},1}

julia> result2[1:3]
3-element Array{Array{Array{Float64,1},1},1}:
 [[-0.4840916699908168, -0.8750172884266356], [-0.831904725657279, 0.554918487193468]]
 [[-0.9979849536743425, 0.06345102236859917], [0.5252178600449688, -0.8509678016762935]]
 [[-0.8602211106475295, 0.5099212103809087], [0.7844127443972155, -0.6202391848530926]]

julia> # Now instead of an Array{Array{Array{Float64,1},1},1} we want to have an
       # Array{Array{Float64,1},1}
       result3 = solve(
           F,
           S₀;
           parameters = [a, b, c],
           start_parameters = p₀,
           target_parameters = params,
           transform_result = (r, p) -> real_solutions(r),
           flatten = true
       );

julia> typeof(result3)
Array{Array{Float64,1},1}

julia> result3[1:3]
3-element Array{Array{Float64,1},1}:
 [-0.4840916699908168, -0.8750172884266356]
 [-0.831904725657279, 0.554918487193468]
 [-0.9979849536743425, 0.06345102236859917]

julia> # The passed `params` do not directly need to be the target parameters.
       # Instead they can be some more concrete information (e.g. an index)
       # and we can them by using the `transform_parameters` method
       result4 = solve(
           F,
           S₀;
           parameters = [a, b, c],
           start_parameters = p₀,
           target_parameters = 1:100,
           transform_result = (r, p) -> (real_solutions(r), p),
           transform_parameters = _ -> rand(3)
       );

julia> typeof(result4)
Array{Tuple{Array{Array{Float64,1},1},Int64},1}

julia> result4[1:3]
3-element Array{Tuple{Array{Array{Float64,1},1},Int64},1}:
 ([[-0.9459225556889465, -0.3243925378902725], [0.39166252055983675, -0.9201089446303168]], 1)
 ([[-0.838341089344409, 0.5451460519134572], [0.4299213699880814, -0.9028663332008627]], 68)
 ([[-0.8261686468010706, -0.5634229024834612], [-0.19003485030903675, -0.9817773452611452]], 69)
```
"""
function solve(args...; kwargs...)
    solve_kwargs, rest = splitkwargs(kwargs, solve_supported_keywords)
    solver, start_solutions = solver_startsolutions(args...; rest...)
    solve(solver, start_solutions; solve_kwargs...)
end
function solve(
    F::Union{MPPolyInputs,AbstractHomotopy,AbstractSystem},
    starts;
    target_parameters = nothing,
    transform_parameters = nothing,
    kwargs...,
)
    if target_parameters === nothing || (
        target_parameters isa AbstractVector{<:Number} && transform_parameters === nothing
    )
        solve_kwargs, rest = splitkwargs(kwargs, solve_supported_keywords)
        if target_parameters !== nothing
            solver, start_solutions = solver_startsolutions(
                F,
                starts;
                target_parameters = target_parameters, rest...,
            )
        else
            solver, start_solutions = solver_startsolutions(F, starts; rest...)
        end
        solve(solver, start_solutions; solve_kwargs...)
    else
        many_parameters_solve(
            F,
            starts;
            target_parameters = target_parameters,
            transform_parameters = transform_parameters === nothing ? identity :
                                   transform_parameters,
            kwargs...,
        )
    end
end

function solve(
    solver::Solver,
    start_solutions;
    show_progress::Bool = true,
    stop_early_cb = nothing,
    kwargs...,
)
    if show_progress
        n = length(start_solutions)
        progress = ProgressMeter.Progress(
            n;
            dt = 0.1,
            desc = "Tracking $n paths... ",
            clear_output_ijulia = true,
            barlen = 40,
            delay = 0.3,
        )
    else
        progress = nothing
    end
    solve(solver, start_solutions, progress, stop_early_cb; kwargs...)
end

function solve(
    solver::Solver,
    start_solutions,
    progress::Union{Nothing,ProgressMeter.Progress},
    stop_early_cb = nothing;
    path_result_details::Symbol = :default,
    save_all_paths::Bool = false,
    path_jumping_check::Bool = true,
    threading::Bool = true,
    start_parameters = nothing,
    target_parameters = nothing,
)
    @unpack trackers, stats = solver

    Threads.resize_nthreads!(trackers)
    Threads.resize_nthreads!(stats)

    start_parameters !== nothing && start_parameters!(solver, start_parameters)
    target_parameters !== nothing && target_parameters!(solver, target_parameters)

    S = collect_startsolutions(start_solutions)
    n = length(S)

    for i = 1:Threads.nthreads()
        prepare!(trackers[i], start_solutions)
        init!(stats[i])
    end

    results = Vector{Union{Nothing,result_type(trackers[1])}}(undef, n)
    results .= nothing

    retracked_paths = 0
    n_blas_threads = single_thread_blas()
    try
        track_parallel!(
            results,
            trackers,
            S,
            1:n,
            stop_early_cb,
            stats,
            progress;
            path_result_details = path_result_details,
            threading = threading,
        )

        if path_jumping_check
            tol = accuracy(trackers[1])
            indices = path_jumping_candidates(results, tol)
            max_correctors = min(max_corrector_iters(trackers[1]) - 1, 2)
            retracked_paths = 0
            while !isempty(indices) && max_correctors > 1
                acc = min(accuracy(trackers[1]), 1e-6)
                track_parallel!(
                    results,
                    trackers,
                    S,
                    indices;
                    threading = threading,
                    max_corrector_iters = min(max_correctors - 1, 2),
                    accuracy = acc,
                )
                retracked_paths += length(indices)
                new_indices = path_jumping_candidates(results, tol)
                if isempty(new_indices) || new_indices ⊆ indices
                    break
                end
                indices = new_indices
            end
        end
    catch e
        if !isa(e, InterruptException)
            rethrow()
        end
    end

    n_blas_threads > 1 && set_num_BLAS_threads(n_blas_threads)
    ntracked = count(!isnothing, results)
    final_results = remove_nothings(results)
    if !save_all_paths
        filter!(
            r -> r.return_code == :success || r.return_code == :terminated_invalid_startvalue,
            final_results,
        )
    end
    Result(
        final_results,
        ntracked,
        seed(trackers[1]);
        retracked_paths = retracked_paths,
        multiplicity_tol = min_accuracy(trackers[1]),
    )
end

start_parameters!(solver::Solver, p::AbstractVector) =
    start_parameters!.(solver.trackers, Ref(p))
target_parameters!(solver::Solver, p::AbstractVector) =
    target_parameters!.(solver.trackers, Ref(p))

function track_parallel!(
    results,
    trackers,
    S,
    range,
    stop_early_cb = nothing,
    stats = nothing,
    progress = nothing;
    threading::Bool = true,
    max_corrector_iters::Union{Nothing,Int} = nothing,
    accuracy::Union{Nothing,Float64} = nothing,
    path_result_details::Symbol = :default,
)
    if threading && Threads.nthreads() > 1
        nthreads = Threads.nthreads()
        ntrackeds = fill(0, nthreads)
        interrupted = Ref(false)
        last_printed = Ref(0)

        @static if VERSION > v"1.3-"
            # create jobs channel
            jobs = Channel{Int}(length(range))
            for i in range
                put!(jobs, i)
            end
            close(jobs)

            # spawn workers which will track the paths
            @sync for _ = 1:nthreads
                Threads.@spawn begin
                    try
                        tid = Threads.threadid()
                        for i in jobs
                            results[i] = track(
                                trackers[tid],
                                S[i];
                                path_number = i,
                                details = path_result_details,
                                max_corrector_iters = max_corrector_iters,
                                accuracy = accuracy,
                            )
                            if is_success(results[i]) && stop_early_cb !== nothing
                                interrupted[] = stop_early_cb(results[i])
                            end

                            stats !== nothing && update!(stats[tid], results[i])
                            ntrackeds[tid] += 1

                            if progress !== nothing
                                ntracked = sum(ntrackeds)
                                if tid == 1 && ntracked - last_printed[] - 32nthreads > 0
                                    last_printed[] = ntracked
                                    update_progress!(
                                        progress,
                                        ntracked,
                                        consolidated_stats(stats),
                                    )
                                end
                            end
                            if interrupted[]
                                break
                            end
                        end
                    catch e
                        if isa(e, InterruptException)
                            interrupted[] = true
                        else
                            rethrow()
                        end
                    end
                end
            end
        else
            Threads.@threads for i in range
                try
                    tid = Threads.threadid()
                    results[i] = track(
                        trackers[tid],
                        S[i];
                        path_number = i,
                        max_corrector_iters = max_corrector_iters,
                        details = path_result_details,
                        accuracy = accuracy,
                    )
                    if is_success(results[i]) && stop_early_cb !== nothing
                        interrupted[] = stop_early_cb(results[i])
                    end

                    stats !== nothing && update!(stats[tid], results[i])
                    ntrackeds[tid] += 1

                    if progress !== nothing
                        ntracked = sum(ntrackeds)
                        if tid == 1 && ntracked - last_printed[] - 32nthreads > 0
                            last_printed[] = ntracked
                            update_progress!(progress, ntracked, consolidated_stats(stats))
                        end
                    end
                    if interrupted[]
                        break
                    end
                catch e
                    if isa(e, InterruptException)
                        interrupted[] = true
                    else
                        rethrow()
                    end
                end
            end
        end

        if progress !== nothing && last_printed != sum(ntrackeds)
            update_progress!(progress, sum(ntrackeds), consolidated_stats(stats))
        end
    else
        ntracked = 0
        for i in range
            results[i] = track(
                trackers[1],
                S[i];
                path_number = i,
                max_corrector_iters = max_corrector_iters,
                details = path_result_details,
                accuracy = accuracy,
            )
            stats !== nothing && update!(stats[1], results[i])
            ntracked += 1
            if ntracked % 32 == 0 && stats !== nothing
                update_progress!(progress, ntracked, stats[1])
            end

            if is_success(results[i]) &&
               stop_early_cb !== nothing && stop_early_cb(results[i])
                break
            end
        end
        if ntracked % 32 != 0 && stats !== nothing
            update_progress!(progress, ntracked, stats[1])
        end
    end
    results
end

prepare!(PT::PathTracker, S) = PT
collect_startsolutions(x::AbstractVector) = x
collect_startsolutions(x) = collect(x)

function update_progress!(progress, ntracked, stats)
    progress === nothing && return nothing

    nsols = stats.regular + stats.singular
    nreal = stats.regular_real + stats.singular_real

    showvalues = (
        ("# paths tracked", ntracked),
        ("# non-singular solutions (real)", "$(stats.regular) ($(stats.regular_real))"),
        ("# singular solutions (real)", "$(stats.singular) ($(stats.singular_real))"),
        ("# total solutions (real)", "$nsols ($nreal)"),
    )

    ProgressMeter.update!(progress, ntracked; showvalues = showvalues)
    nothing
end

function remove_nothings(v::Vector{Union{Nothing,T}}) where {T}
    w = T[]
    for vᵢ in v
        !isnothing(vᵢ) && push!(w, vᵢ)
    end
    w
end

#####################
## Many parameters ##
#####################

export ManyParametersSolver, many_parameters_solver, many_parameters_solve

"""
    ManyParametersSolver(solver)

A wrapper around a solver to prepare it for [`many_parameters_solve`](@ref).
"""
struct ManyParametersSolver{S<:Solver}
    solvers::Vector{S}
end
ManyParametersSolver(S::Solver) = ManyParametersSolver([S])

const many_parameters_solve_supported_keywords =
    [:transform_result, :transform_parameters, :flatten, :path_jumping_check, :threading]

# dispatch for the many parameters case
function solver(
    F::Union{MPPolyInputs,AbstractHomotopy,AbstractSystem},
    starts;
    target_parameters = nothing,
    transform_parameters = nothing,
    kwargs...,
)
    if target_parameters === nothing || (
        target_parameters isa AbstractVector{<:Number} && transform_parameters === nothing
    )
        solve_kwargs, rest = splitkwargs(kwargs, solve_supported_keywords)
        if target_parameters !== nothing
            solver, _ = solver_startsolutions(
                F,
                starts;
                target_parameters = target_parameters, rest...,
            )
        else
            solver, _ = solver_startsolutions(F, starts; rest...)
        end
    else
        solver = many_parameters_solver(
            F,
            starts;
            target_parameters = target_parameters,
            transform_parameters = transform_parameters === nothing ? identity :
                                       transform_parameters,
            kwargs...,
        )
    end
    solver
end

"""
    many_parameters_solver( F,
         S::Vector{<:AbstractVector},
         start_params::AbstractVector{<:Number},
         params::AbstractVector; options...)

Construct a `ManyParametersSolver` by reusing the [`solver`](@ref) function.
"""
function many_parameters_solver(
    F,
    starts::Vector{<:AbstractVector};
    generic_parameters = nothing,
    start_parameters = generic_parameters,
    target_parameters = throw(ArgumentError("`target_parameters = ...` needs to be specified")),
    kwargs...,
)
    if start_parameters === nothing
        throw(ArgumentError("`start_parameters = ...` needs to be provided."))
    end
    ManyParametersSolver(solver(
        F,
        starts;
        generic_parameters = start_parameters,
        kwargs...,
    ))
end

function many_parameters_solve(
    F,
    starts::Vector{<:AbstractVector};
    generic_parameters = nothing,
    start_parameters = generic_parameters,
    target_parameters = throw(ArgumentError("`target_parameters = ...` needs to be specified")),
    kwargs...,
)
    solve_kwargs, rest = splitkwargs(kwargs, many_parameters_solve_supported_keywords)
    solver = many_parameters_solver(
        F,
        starts;
        start_parameters = start_parameters, target_parameters = target_parameters, rest...,
    )
    many_parameters_solve(
        solver,
        starts,
        start_parameters,
        target_parameters;
        solve_kwargs...,
    )
end

function solve(
    MPS::ManyParametersSolver,
    starts;
    generic_parameters = nothing,
    start_parameters = generic_parameters,
    target_parameters = throw(ArgumentError("`target_parameters = ...` needs to be specified")),
    kwargs...,
)
    many_parameters_solve(MPS, starts, start_parameters, target_parameters; kwargs...)
end

function many_parameters_solve(
    MPS::ManyParametersSolver,
    starts::Vector{<:AbstractVector},
    start_params::Union{Nothing,AbstractVector{<:Number}},
    parameters::AbstractVector;
    transform_result::Function = tuple,
    transform_parameters::Function = identity,
    flatten::Bool = false,
    path_jumping_check::Bool = true,
    threading::Bool = true,
)
    Threads.resize_nthreads!(MPS.solvers)

    if !isnothing(start_params)
        start_parameters!.(MPS.solvers, Ref(start_params))
    end

    first_result = solve(
        MPS.solvers[1],
        starts;
        target_parameters = transform_parameters(parameters[1]),
        threading = false,
        show_progress = false,
        path_jumping_check = path_jumping_check,
    )

    if flatten
        results = transform_result(first_result, parameters[1])
        if !(results isa AbstractArray)
            throw(ArgumentError("Cannot flatten arguments of type `$(typeof(results))`"))
        end
    else
        results = [transform_result(first_result, parameters[1])]
    end
    many_parameters_solve_barrier!(
        results,
        MPS,
        starts,
        parameters,
        transform_result,
        transform_parameters,
        Val(flatten);
        path_jumping_check = path_jumping_check,
        threading = threading,
    )
    results
end

function many_parameters_solve_barrier!(
    results,
    MPS::ManyParametersSolver,
    starts::Vector{<:AbstractVector},
    parameters::AbstractVector,
    transform_result::Function,
    transform_parameters::Function,
    ::Val{Flatten};
    path_jumping_check::Bool = false,
    threading::Bool = true,
) where {Flatten}
    N = length(parameters)

    if threading && Threads.nthreads() > 1
        thread_results = [similar(results, 0) for _ = 2:Threads.nthreads()]
        push!(thread_results, results)

        Threads.@threads for i = 2:N
            tid = Threads.threadid()
            rᵢ = transform_result(
                solve(
                    MPS.solvers[tid],
                    starts;
                    target_parameters = transform_parameters(parameters[i]),
                    path_jumping_check = path_jumping_check,
                    threading = false,
                    show_progress = false,
                ),
                parameters[i],
            )

            if Flatten
                append!(thread_results[tid], rᵢ)
            else
                push!(thread_results[tid], rᵢ)
            end
        end

        for i = 1:(Threads.nthreads()-1)
            append!(results, thread_results[i])
        end
    else
        for i = 2:N
            rᵢ = transform_result(
                solve(
                    MPS.solvers[1],
                    starts;
                    target_parameters = transform_parameters(parameters[i]),
                    path_jumping_check = path_jumping_check,
                    threading = false,
                    show_progress = false,
                ),
                parameters[i],
            )

            if Flatten
                append!(results, rᵢ)
            else
                push!(results, rᵢ)
            end
        end
    end

    results
end


text:
using DataFrames

export
    LinearRegressionModel,
    fit!,
    predict

mutable struct LinearRegressionModel{T, S}
    y::Vector{T}
    xs::AbstractArray{S}
    argv::Vector{Float64}
    n::Integer
end

function LinearRegressionModel(
    df::DataFrame,
    label::Symbol,
    features::Vector{Symbol};
    β::Vector{<:Real}=rand(length(features)+1)
)
    n = nrow(df)
    y = df[!, label]
    xs = Matrix(df[!, features])
    if length(β) != length(features) + 1
        throw(DimensionMismatch("Number of features and arguments mismatch."))
    end

    return LinearRegressionModel(y, xs, β, n)
end

function LinearRegressionModel(
    df::DataFrame,
    label::Symbol,
    feature::Symbol;
    β::Vector{<:Real}=rand(2)
)
    n = nrow(df)
    y = df[!, label]
    x = Matrix(df[!, [feature]])
    if length(β) != 2
        throw(DimensionMismatch("Number of features and arguments mismatch."))
    end

    return LinearRegressionModel(y, x, β, n)
end

function y(model::LinearRegressionModel, xs::Matrix{T}) where {T<:Real}
    xs = hcat(ones(T, size(xs, 1)), xs)

    return xs * model.argv
end

residual(model::LinearRegressionModel) = model.y .- y(model, model.xs)

function loss(model::LinearRegressionModel)
    l = sum(x -> 0.5 * x^2, residual(model))

    return l/model.n
end

function ∇L(model::LinearRegressionModel)
    xs = hcat(ones(size(model.xs, 1)), model.xs)

    return vec(sum(-residual(model) .* xs, dims=1))
end

function gradient_descent(model::LinearRegressionModel, η::Real=1e-4, atol::Real=1e-6, show::Bool=false)
    β = model.argv
    while (l = loss(model)) > atol
        show && println("Loss: $l")
        β .-= η .* ∇L(model)
    end

    return β
end

function fit!(model::LinearRegressionModel; method=gradient_descent, η::Real=1e-4, atol::Real=1e-6, show=false)
    β = method(model, η, atol, show)
    model.argv .= β

    return model
end

function predict(model::LinearRegressionModel, xsᵢ::Vector{T}) where {T<:Real}
    xsᵢ = vcat(ones(T, 1), xsᵢ)'

    return xsᵢ * model.argv
end

function predict(model::LinearRegressionModel, xsᵢ::T) where {T<:Real}
    xsᵢ = vcat(ones(T, 1), xsᵢ)'

    return xsᵢ * model.argv
end

ext: jl


text:
# Update_CovHypers.jl

function llik_ScInvChiSq(y::Vector{T}, ν::T, s0::T) where T <: Real
    sum1divy = sum(1.0 ./ y)
    sumlogy = sum(log.(y))
    n = length(y)

    halfnu = 0.5 * ν
    a1 = halfnu * n * (log(ν) + log(s0) - log(2.0))
    a2 = n * SpecialFunctions.loggamma(halfnu)
    a3 = halfnu * sumlogy
    a4 = halfnu * s0 * sum1divy

    return a1 - a2 - a3 - a4
end

function rpost_ScInvChiSq_s0(y::Vector{T}, ν::T, prior::GammaParams) where T <: Real
    sum1divy = sum(1.0 ./ y)
    n = length(y)

    halfnu = 0.5 * ν
    a1 = prior.shape + halfnu * n
    b1 = prior.rate + halfnu * sum1divy

    return rand(Distributions.Gamma(a1, 1.0/b1))
end

function rpost_ScInvChiSq_ν(y::Vector{T}, s0::T, prior::Vector{T}) where T <: Real
    ## Here the prior is a vector of values representing a discrete uniform over them.
    lp = [ llik_ScInvChiSq(y, nu, s0) for nu in prior ]
    w = exp.( lp .- maximum(lp) )
    return StatsBase.sample(prior, Weights(w))
end

function update_cov_hypers!(state::State_GPMTD, prior::PriorCovHyper_Matern) where T <: Real

    L = length(state.mixcomps)

    nζ = StatsBase.counts(state.ζ, 1:L) # intercept does not involve these parameters
    use_indx = findall( nζ .> 0 )

    if length(use_indx) > 0

        κ_dat = [ state.mixcomps[j].κ for j in use_indx ]
        lenscale_dat = [ state.mixcomps[j].corParams.lenscale for j in use_indx ]

        κ_ν_out = rpost_ScInvChiSq_ν(κ_dat, state.mixcomps[1].κ_hypers.κ0, prior.κ_ν)
        κ0_out = rpost_ScInvChiSq_s0(κ_dat, κ_ν_out, prior.κ0)

        lenscale_ν_out = rpost_ScInvChiSq_ν(lenscale_dat, state.mixcomps[1].corHypers.lenscale0, prior.lenscale_ν)
        lenscale0_out = rpost_ScInvChiSq_s0(lenscale_dat, lenscale_ν_out, prior.lenscale0)

    else

        κ_ν_out = StatsBase.sample(prior.κ_ν)
        κ0_out = rand(Distributions.Gamma(prior.κ0.shape, 1.0/prior.κ0.rate))

        lenscale_ν_out = StatsBase.sample(prior.lenscale_ν)
        lenscale0_out = rand(Distributions.Gamma(prior.lenscale0.shape, 1.0/prior.lenscale0.rate))

    end

    # update these parameters in all mixcomps
    for ℓ = 1:L
        state.mixcomps[ℓ].κ_hypers.κ_ν = κ_ν_out
        state.mixcomps[ℓ].κ_hypers.κ0 = κ0_out
        state.mixcomps[ℓ].corHypers.lenscale_ν = lenscale_ν_out
        state.mixcomps[ℓ].corHypers.lenscale0 = lenscale0_out
    end

    return nothing
end

### tests passed 4/12/19
# ν = 7.5
# s0 = 10.0
# using Distributions
# using StatsBase
# using SpecialFunctions
# n = 1000
#
# xx = rand(InverseGamma(0.5*ν, 0.5*ν*s0), n)
#
# nsamp = 1000
# samp = ones(nsamp, 2)
# for i = 2:nsamp
#     samp[i,1] = rpost_ScInvChiSq_ν(xx, samp[i-1,2], [5.0, 7.5, 10.0, 25.0, 50.0])
#     samp[i,2] = rpost_ScInvChiSq_s0(xx, samp[i,1], GammaParams(1.0, 1.0))
# end
#
# using Plotly
# plot(samp[:,1])
# plot(samp[:,2])


text:
unsafeattribute(::Nothing, ::Core.AbstractString) = nothing
unsafeattribute(x::LightXML.XMLElement, name::Core.AbstractString) = attribute(x, name)

function parse_scalar(xel, name::String, T; default::Union{String,Nothing}=nothing)
    scalarstr = unsafeattribute(xel, name)
    if scalarstr === nothing
        if default === nothing
            @error "no parsable scalar found"
        else
            scalarstr = default
        end
    end

    return parse(T, scalarstr)
end

function parse_vector(xel, name::String, T; default::Union{String,Nothing}=nothing)
    vectorstr = unsafeattribute(xel, name)
    if vectorstr === nothing
        if default === nothing
            @error "no parsable vector found"
        else
            vectorstr = default
        end
    end

    return parse.(T,split(vectorstr))
end

function parse_inertiamatrix(xinertia, T)
    if xinertia === nothing
        J = zeros(T, 3, 3)
    else
        ixx = parse_scalar(xinertia, "ixx", T)
        ixy = parse_scalar(xinertia, "ixy", T, default = "0")
        ixz = parse_scalar(xinertia, "ixz", T, default = "0")
        iyy = parse_scalar(xinertia, "iyy", T)
        iyz = parse_scalar(xinertia, "iyz", T, default = "0")
        izz = parse_scalar(xinertia, "izz", T)
        J = [ixx ixy ixz; ixy iyy iyz; ixz iyz izz]
    end

    return J
end

function parse_pose(xpose, T)
    if xpose === nothing
        x, q = zeros(T, 3), one(UnitQuaternion{T})
    else
        x = parse_vector(xpose, "xyz", T, default = "0 0 0")
        rpy = parse_vector(xpose, "rpy", T, default = "0 0 0")
        q = UnitQuaternion(RotZYX(rpy[3], rpy[2], rpy[1]))
    end

    return x, q
end

function parse_inertia(xinertial, T)
    if xinertial === nothing
        x = zeros(T, 3)
        q = one(UnitQuaternion{T})
        m = zero(T)
        J = zeros(T, 3, 3)
    else
        x, q = parse_pose(find_element(xinertial, "origin"), T)
        J = parse_inertiamatrix(find_element(xinertial, "inertia"), T)
        m = parse_scalar(find_element(xinertial, "mass"), "value", T, default = "0")
    end

    return x, q, m, J
end

function parse_robotmaterials(xroot, T)
    xmaterials = get_elements_by_tagname(xroot, "material")

    mdict = Dict{String,Vector{T}}()

    for xmaterial in xmaterials
        name = attribute(xmaterial, "name")
        colornode = find_element(xmaterial, "color")
        colorvec = parse_vector(colornode, "rgba", T)
        mdict[name] = colorvec
    end

    return mdict
end

function parse_xmaterial(xmaterial, materialdict, T)
    if xmaterial === nothing
        color = RGBA(0.75, 0.75, 0.75)
    else
        name = unsafeattribute(xmaterial, "name")
        if haskey(materialdict, name)
            colorvec = materialdict[name]
        else
            colorvec = [0.75; 0.75; 0.75; 1.0]
        end
        colornode = find_element(xmaterial, "color")
        colorvec = parse_vector(colornode, "rgba", T, default = string(colorvec[1]," ",colorvec[2]," ",colorvec[3]," ",colorvec[4]))

        color = RGBA(colorvec...)
    end

    return color
end

function parse_shape(xvisual, materialdict, T)
    if xvisual === nothing
        shape = nothing
    else
        xgeometry = find_element(xvisual, "geometry")
        @assert xgeometry !== nothing

        color = parse_xmaterial(find_element(xvisual, "material"), materialdict, T)
        x, q = parse_pose(find_element(xvisual, "origin"), T)

        shapenodes = LightXML.XMLElement[]
        for node in child_nodes(xgeometry)  # node is an instance of XMLNode
            if is_elementnode(node)
                push!(shapenodes, XMLElement(node))
            end
        end

        if length(shapenodes) == 0
            shape = nothing
        else
            if length(shapenodes) > 1
                @info "Multiple geometries."
            end

            shapenode = shapenodes[1]
            if name(shapenode) == "box"
                xyz = parse_vector(shapenode, "size", T, default = "1 1 1")
                shape = Box(xyz..., zero(T), color = color, xoffset = x, qoffset = q)
            elseif name(shapenode) == "cylinder"
                r = parse_scalar(shapenode, "radius", T, default = "0.5")
                l = parse_scalar(shapenode, "length", T, default = "1")
                shape = Cylinder(r, l, zero(T), color = color, xoffset = x, qoffset = q)
            elseif name(shapenode) == "sphere"
                r = parse_scalar(shapenode, "radius", T, default = "0.5")
                shape = Sphere(r, zero(T), color = color, xoffset = x, qoffset = q)
            elseif name(shapenode) == "mesh"
                path = attribute(shapenode, "filename")
                shape = Mesh(path, zero(T), zeros(T, 3, 3), color = color, xoffset = x, qoffset = q)
            else
                @info "Unkown geometry."
                shape = nothing
            end
        end

    end

    return shape
end

function parse_link(xlink, materialdict, T)
    xvisual = find_element(xlink, "visual")

    x, q, m, J = parse_inertia(find_element(xlink, "inertial"), T)
    shape = parse_shape(find_element(xlink, "visual"), materialdict, T)
    name = attribute(xlink, "name")

    if shape === nothing
        link = Body(m, J, name=name)
    else
        shape.m = m
        shape.J = J
        link = Body(shape, name=name)
    end

    link.state.xc = x
    link.state.qc = q

    return link, shape
end

function parse_links(xlinks, materialdict, T)
    ldict = Dict{String,AbstractBody{T}}()
    shapes = Shape{T}[]

    for xlink in xlinks
        link, shape = parse_link(xlink, materialdict, T)
        ldict[link.name] = link
        shape !== nothing && push!(shapes, shape)
    end

    return ldict, shapes
end

function parse_joint(xjoint, origin, plink, clink, T)
    jointtype = attribute(xjoint, "type")
    x, q = parse_pose(find_element(xjoint, "origin"), T)
    axis = parse_vector(find_element(xjoint, "axis"), "xyz", T, default = "1 0 0")
    p1 = x
    p2 = zeros(T, 3)
    name = attribute(xjoint, "name")

    # TODO limits for revolute joint?
    if jointtype == "revolute" || jointtype == "continuous"
        joint = EqualityConstraint(Revolute(plink, clink, axis; p1=p1, p2=p2, qoffset = q), name=name)
    elseif jointtype == "prismatic"
        joint = EqualityConstraint(Prismatic(plink, clink, axis; p1=p1, p2=p2, qoffset = q), name=name)
    elseif jointtype == "planar"
        joint = EqualityConstraint(Planar(plink, clink, axis; p1=p1, p2=p2, qoffset = q), name=name)
    elseif jointtype == "fixed"
        joint = EqualityConstraint(Fixed(plink, clink; p1=p1, p2=p2, qoffset = q), name=name)
    elseif jointtype == "floating"
        joint = EqualityConstraint(Floating(plink, clink), name=name)
    end

    return joint
end

function parse_joints(xjoints, ldict, floating, T)
    origins = Origin{T}[]
    links = Body{T}[]
    joints = EqualityConstraint{T}[]
    floatingname = ""

    for name in keys(ldict)
        child = false
        for xjoint in xjoints
            xchild = find_element(xjoint, "child")
            childname = attribute(xchild, "link")
            if childname == name
                child = true
                break
            end
        end
        if child
            push!(links, ldict[name])
        else
            origin = Origin{T}()
            if floating # keep current link and create new origin
                push!(links, ldict[name])
                floatingname = name
            else # make current link origin
                origin.id = ldict[name].id
                origin.name = name
            end
            push!(origins, origin)
        end
    end

    @assert length(origins) == 1 "Multiple origins"
    origin = origins[1]

    for xjoint in xjoints
        xplink = find_element(xjoint, "parent")
        plink = ldict[attribute(xplink, "link")]
        if plink.id == origin.id
            plink = origin
        end
        xclink = find_element(xjoint, "child")
        clink = ldict[attribute(xclink, "link")]

        joint = parse_joint(xjoint, origin, plink, clink, T)
        push!(joints, joint)
    end

    if floating
        originjoint = EqualityConstraint(Floating(origin, ldict[floatingname]),name="autoorigincon")
        push!(joints, originjoint)
    end

    return origin, links, joints
end

function parse_urdf(filename, floating, ::Type{T}) where T
    xdoc = LightXML.parse_file(filename)
    xroot = LightXML.root(xdoc)
    @assert LightXML.name(xroot) == "robot"

    xlinks = get_elements_by_tagname(xroot, "link")
    xjoints = get_elements_by_tagname(xroot, "joint")

    materialdict = parse_robotmaterials(xroot, T)
    ldict, shapes = parse_links(xlinks, materialdict, T)
    origin, links, joints = parse_joints(xjoints, ldict, floating, T)

    free(xdoc)

    return origin, links, joints, shapes
end


text:
@testset "Convolution methods" begin
    function conv_test(psf, img, img_out, dims, s)
        otf = fft(psf, dims)
        otf_r = rfft(psf, dims)
        otf_p, conv_p = plan_conv_r(psf, dims)
        @testset "$s" begin
            @test img_out ≈ conv_psf(img, psf, dims)
            @test img_out ≈ conv_otf(img, otf, dims)
            @test img_out ≈ conv_otf_r(img, otf_r, dims)
            @test img_out ≈ conv_p(img, otf_p)
        end
    end

    N = 5
    psf = zeros((N, N))
    psf[1, 1] = 1
    img = randn((N, N))
    conv_test(psf, img, img, [1,2], "Convolution random image with delta peak")

    N = 5
    psf = abs.(randn((N, N, 2)))
    img = randn((N, N, 2))
    dims = [1, 2]
    img_out = conv_psf(img, psf, dims)
    conv_test(psf, img, img_out, dims, "Convolution with random 3D PSF and random 3D image over 2D dimensions")

    N = 5
    psf = abs.(randn((N, N, N, N, N)))
    img = randn((N, N, N, N, N))
    dims = [1, 2, 3, 4]
    img_out = conv_psf(img, psf, dims)
    conv_test(psf, img, img_out, dims, "Convolution with random 5D PSF and random 5D image over 4 Dimensions")

    N = 5
    psf = abs.(zeros((N, N, N, N, N)))
    for i = 1:N
        psf[1,1,1,1, i] = 1
    end
    img = randn((N, N, N, N, N))
    dims = [1, 2, 3, 4]
    img_out = conv_psf(img, psf, dims)
    conv_test(psf, img, img, dims, "Convolution with 5D delta peak and random 5D image over 4 Dimensions")

    function conv_test_complex(psf, img, img_out, dims, s)
        otf = fft(psf, dims)
        @testset "$s" begin
            @test img_out ≈ conv(img, psf, dims, real_res=false)
            @test img_out ≈ conv_v_ft(img, otf, dims, real_res=false)
        end
    end

    N = 5
    psf = abs.(zeros((N, N, N, N, N)))
    for i = 1:N
        psf[1,1,1,1, i] = 1
    end
    img = randn((N, N, N, N, N))
    img = img .+ exp.(1im .* img)
    dims = [1, 2, 3, 4]
    img_out = conv(img, psf, dims, real_res=false)
    conv_test_complex(psf, img, img, dims, "Complex Convolution with 5D delta peak and random 5D image over 4 Dimensions")

    N = 5
    psf = zeros((N, N))
    psf[1, 1] = 1
    img = randn((N, N))
    img = img .+ exp.(1im .* img)
    conv_test_complex(psf, img, img, [1,2], "Complex Convolution random image with delta peak")

    N = 5
    psf = abs.(randn((N, N, 2)))
    img = randn((N, N, 2))
    img = img .+ exp.(1im .* img)
    dims = [1, 2]
    img_out = conv(img, psf, dims, real_res=false)
    conv_test_complex(psf, img, img_out, dims, "Complex Convolution with random 3D PSF and random 3D image over 2D dimensions")

    N = 5
    psf = abs.(randn((N, N, N, N, N)))
    img = randn((N, N, N, N, N))
    img = img .+ exp.(1im .* img)
    dims = [1, 2, 3, 4]
    img_out = conv(img, psf, dims, real_res=false)
    conv_test_complex(psf, img, img_out, dims, "Complex Convolution with random 5D PSF and random 5D image over 4 Dimensions")

    @testset "dims argument nothing" begin
        N = 5
        psf = abs.(randn((N, N, N, N, N)))
        img = randn((N, N, N, N, N))
        dims = [1,2,3,4,5]
        @test conv(psf, img) ≈ conv(psf, img, dims)
        @test conv_v_ft(img, fft(psf, dims)) ≈ conv_v_ft(img, fft(psf, dims), dims)
    end

    @testset "Test shift argument" begin
        N = 10
        img = randn((N, N))
        psf = abs.(randn((N, N)))

        @test conv(img, psf, shift=true) ≈ conv(img, ifftshift(psf))
        @test conv_psf(img, psf, shift=true) ≈ conv(img, ifftshift(psf))
        @test conv_psf(img, psf, shift=true) ≈ conv(img, ifftshift(psf))
    end

    @testset "Check types" begin
        N = 10
        img = randn((N, N))
        psf = abs.(randn((N, N)))

        @test conv_psf(img, psf, shift=true) ≈ conv(img, ifftshift(psf))
        @test typeof(conv_psf(img, psf, shift=true)) != typeof(conv(img, ifftshift(psf)))
        @test typeof(conv_psf(img, psf, shift=true)) == typeof(conv(img, ifftshift(psf), real_res=true))
    end
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
using Replace
using Test
using MacroTools

module Sine
  sine(x) = sin(x)
end
@testset "Replace" begin

@testset "Basic" begin
  @test 1.0 == @replace sin cos sin(0.0)
  @test 0.0 == @replace cos sin cos(0.0)
  # make sure that we haven't clobbered the definition of sin and cos
  @assert cos(0.0) == 1.0 && sin(0.0) == 0.0

  replacement() = pi
  @test pi + exp(0) == @replace rand replacement rand() + exp(0)

  x, y = rand(2)
  mapping = Dict(cos=>sin, tan=>exp)
  f(x, y) = sin(x) + exp(y)
  g(x, y) = cos(x) + tan(y)
  @test f(x, y) == @replace mapping g(x, y)
end

@testset "In a function" begin
  function foo()
    x = @replace cos sin cos(0.0)
  end
  @test foo() == 0.0

  function bar()
    x = @replace Dict(tan=>exp) tan(1.0)
  end
  @test bar() == exp(1.0)
end

@testset "Function declared outside another" begin
  foo(x) = cos(x)
  function bar(x)
    return @replace cos sin foo(x)
  end
  @test bar(0.0) == 0.0

  function baz(x)
    return @replace Dict(cos=>sin) foo(x)
  end
  @test baz(0.0) == 0.0
end

@testset "Inside a module" begin
  @test 1.0 == @replace sin cos Sine.sine(0.0)
end

@testset "How Cassette is supposed to work" begin
  replacement() = true
  original() = false
  @eval Cassette.@context Ctx
  Cassette.overdub(::Ctx, fn::typeof(original), args...) = replacement(args...)
  function functionbarrier()
    x = Cassette.overdub(Ctx(), original)
    return x
  end
  @test functionbarrier()
end

@testset "Some things inside a function" begin
  replacement() = true
  original() = false
  function functionbarrier()
    x = @replace Dict(original=>replacement) original()
    return x
  end
  @test_broken functionbarrier()
  @test functionbarrier()
end

@testset "Everything inside a function" begin
  function functionbarrier()
    replacement() = true
    original() = false
    x = @replace Dict(original=>replacement) original()
    return x
  end
  @test_broken functionbarrier()
  @test functionbarrier()
end

end


text:
include("IOM/src/share/CyclicData.jl")

using .CyclicData
using CFTime
using NCDatasets
using Dates

test_file = "test.nc"
println("Making test data $(test_file)")

period = 2*365
t = (range(0, period, length=51) |> collect)[1:end-1]
t0 = 0.0
v = cos.( (t .- t0) / 365.0 * 2π) + t * 0.001

# v is of shape (x, y, time)
v = reshape(v, 1, 1, :)

Dataset(test_file, "c") do ds

    defDim(ds, "Nx", 1)
    defDim(ds, "Ny", 1)
    defDim(ds, "time", Inf)

    for (varname, vardata, vardim, attrib) in [
        ("v",     v, ("Nx", "Ny", "time",), Dict()),
        ("time",  t, ("time",), Dict(
            "calendar" => "noleap",
            "units"    => "days since 0000-01-01 00:00:00",
        )),
    ]
        println("Doing varname:", varname)
        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

println("Making CyclicDataManager")
_t = DateTimeNoLeap(3,1, 1)
cdm = CyclicDataManager(;
    timetype = DateTimeNoLeap,
    var_file_map = Dict("v" => test_file),
    beg_time = DateTimeNoLeap(0, 1, 1),
    end_time = DateTimeNoLeap(2, 1, 1),
    align_time = _t,
)

d = makeDataContainer(cdm)

println("Interpolating data")
dt = Second(86400)
interp_d = []
interp_t = []

for i = 1:365*5

    global _t

    interpData!(
        cdm,
        _t,
        d,
    )

    push!(interp_d, d["v"][1,1])
    push!(interp_t, _t)

    _t = _t + dt
end

interp_t_day = [ timeencode(_t, "days since 0001-01-01 00:00:00", DateTimeNoLeap) for (_, _t) in enumerate(interp_t) ]

println("Plotting")

using PyPlot
plt = PyPlot

plt.plot(interp_t_day, interp_d)
plt.scatter(t, v, s=20, marker="o")

for i=1:5
    plt.scatter(t .+ period*(i-1), v, s=2, marker="o", color="black")
end

plt.show()

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, Random
CairoMakie.activate!() # HIDE
let
    Random.seed!(123)
    x = -3:0.05:3
    y = exp.(-x.^2)
    n = 15
    y[1:n] = y[1:n] .+ 0.02*randn(n)
    fig = Figure(resolution = (700, 450))
    ax1 = Axis(fig, xlabel = "x", ylabel = "f(x)", xgridvisible = true,
              ygridvisible = true)
    ax2 = Axis(fig, bbox = BBox(140, 260, 260, 350), xticklabelsize = 12,
            yticklabelsize = 12, showgrid=false,
            title = "inset  at (100, 300, 100, 200)")
    line1 = lines!(ax1, x, y, color = :red)
    # inset
    lines!(ax2, x, y, color = :red)
    limits!(ax2, -3.1,-1.9,-0.05,0.05)
    ax2.yticks = [-0.05,0,0.05]
    ax2.xticks = [-3,-2.5,-2]
    translate!(ax2.scene, 0, 0, 10)
    translate!(ax2.elements[:background], 0, 0, 9)
    translate!(ax2.elements[:xgridlines], 0, 0, 9)
    translate!(ax2.elements[:ygridlines], 0, 0, 9)

    leg = Legend(fig, [line1], ["f(x)"], halign = :right, valign = :top,
    tellheight =false, tellwidth = false, margin = (10, 10, 10, 10),)
    fig[1, 1] = ax1
    fig[1, 1] = leg
    fig
    #save("Fig2Lines.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "LineWithInset.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using DataFrames, Gadfly
ind = collect(1:28)
mydata = DataFrame(hcat(ind, rmsefitted2, rmsefitted2_sigma10))
stack(mydata, [:x2, :x3])

myplot = plot(stack(mydata, [:x2, :x3]),
    x = :x1,
    y = :value,
    color = :variable,
    Geom.point,
    Geom.line)
draw(PNG("myplot.png", 3inch, 3inch), myplot)

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
#=
Temporary copy of resources from experimental code transformation package
https://github.com/perrutquist/CodeTransformation.jl until it becomes tagged.
=#

import Core: SimpleVector, svec, CodeInfo
import Base: uncompressed_ast, unwrap_unionall

jl_method_def(argdata::SimpleVector, ci::CodeInfo, mod::Module) =
    ccall(:jl_method_def, Cvoid, (SimpleVector, Any, Ptr{Module}), argdata, ci, pointer_from_objref(mod))

typevars(T::UnionAll) = (T.var, typevars(T.body)...)
typevars(T::DataType) = ()

@nospecialize # the below functions need not specialize on arguments

getmodule(F::Type{<:Function}) = F.name.mt.module
getmodule(f::Function) = getmodule(typeof(f))

makesig(f::Function, args) = Tuple{typeof(f), args...}

argdata(sig) = svec(unwrap_unionall(sig).parameters::SimpleVector, svec(typevars(sig)...))
argdata(sig, f::Function) = svec(svec(typeof(f), unwrap_unionall(sig).parameters[2:end]...), svec(typevars(sig)...))

addmethod!(f::Function, argtypes::Tuple, ci::CodeInfo) = addmethod!(makesig(f, argtypes), ci)

function addmethod!(sig::Type{<:Tuple{F, Vararg}}, ci::CodeInfo) where {F<:Function}
    jl_method_def(argdata(sig), ci, getmodule(F))
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("randExp.jl")

# we write an improoved function that uses Monte Carlo with exponential importance sampling
# to estimate the volume of the unit hypersphere in d dimensions.
# We choose our points x with independent Cartesian components and prob density exp(-lam.xi) / [ lam*((1-exp(-lam)) ].
# We achieve this by taking xi = (1/lam) log (1/y) where y is chosen uniformly in (exp(-lam),1)
function sphereVolMCImpExp(d, n, lam)
	count = 0.0
	for i = 1:n

		# we generate x according to our distribution
		x = randExp(d, lam)

		# we calculate the squared norm of x and check if x is inside the hypersphere
		norm2 = x' * x
		if norm2 < 1

			# we compute the density at x
			rho = 1
			for i = 1:d

				# note that y(i) = exp(-lam * x(i))
				rho *= exp(-lam*x[i]) * lam / (1-exp(-lam))
			end
			count += 1/rho
		end
	end
	vol =  2^d * count / n
	return vol
end


text:
using Test
using BenchmarkTools
using Plots
using SparseArrays

opts = TrajectoryOptimization.SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-4
opts.cost_tolerance_intermediate = 1e-4
opts.constraint_tolerance = 1e-3
opts.square_root = false
opts.active_constraint_tolerance = 0.0
opts.outer_loop_update_type = :default
opts.live_plotting = false

# Block move
model, obj = TrajectoryOptimization.Dynamics.double_integrator
u_min = -0.2
u_max = 0.2
obj_con = TrajectoryOptimization.ConstrainedObjective(obj, tf=5.0,use_xf_equality_constraint=true, u_min=u_min, u_max=u_max)#, x_min=x_min, x_max=x_max)

solver = TrajectoryOptimization.Solver(model,obj_con,integration=:rk4,N=7,opts=opts)
U = rand(solver.model.m, solver.N)

results, stats = TrajectoryOptimization.solve(solver,U)
@assert max_violation(results) < opts.constraint_tolerance

λ_update_default!(results,solver)

J_prev = cost(solver,results)
c_max_prev = stats["c_max"][end]

# Newton 2
results_new = copy(results)
# newton_solve!(results_new,solver)

# # newton_solve!(results_new,solver)
# @assert results_new.X == results.X
# @assert results_new.U == results.U
# @assert results_new.C == results.C
# @assert results_new.λ == results.λ

newton_results = NewtonResults(solver)
newton_active_set!(newton_results,results_new,solver)
update_newton_results!(newton_results,results_new,solver)
newton_step!(results_new,newton_results,solver,1.0)
c_max = max_violation(results_new)
J = cost(solver,results_new)

J-J_prev
# c_max-c_max_prev
#
# cost_constraints(solver,results_new) - cost_constraints(solver,results)
# _cost(solver,results_new) - _cost(solver,results)
#
# @assert results.active_set == results_new.active_set
#
# # project into feasible space
# backwardpass!(results_new,solver)
# rollout!(results_new,solver,1.0)
# results_new.X .= deepcopy(results_new.X_)
# results_new.U .= deepcopy(results_new.U_)
# update_constraints!(results_new,solver)
# JJ = cost(solver,results_new)
# cc_max = max_violation(results_new)
# update_newton_results!(newton_results,results_new,solver)
# norm(newton_results.c)
# norm(vcat(results_new.C...) .+ 0.5*newton_results.s.^2)
# # @assert vcat(results_new.active_set...) == newton_results.active_set
# JJ-J_prev
# cc_max
#
# @assert results_new.active_set == results.active_set
# newton_results.λ[findall(x->x < 0.0, newton_results.λ[newton_results.active_set_ineq])]
#
# 0.5*newton_results.z'*newton_results.∇²J*newton_results.z + newton_results.∇J'*newton_results.z + newton_results.λ'*newton_results.c


text:
# Adaptive MPP:
"""
    AdaptiveMPP{SymbolType}(context_length::Int [,α::Float64 ])

Creates an Adaptive MPP predictor with a context depth of `context_length` and mixing coefficient
``α`` for `SymbolType`. Default value of mixing parameter ``α = 0.10``. `SymbolType` can be
any valid type including `Char`, `Int` etc.,

## Examples
```julia-repl
julia> p = AdaptiveMPP{Char}(4)
DiscretePredictors.adaptiveMPP{Char}([*] (0)
, Char[], 4, 0.1, Inf, 0.0, [0.2, 0.2, 0.2, 0.2, 0.2], Dict{Char,Float64}[Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}()])

julia> p = AdaptiveMPP{Char}(4,0.13)
DiscretePredictors.adaptiveMPP{Char}([*] (0)
, Char[], 4, 0.13, Inf, 0.0, [0.2, 0.2, 0.2, 0.2, 0.2], Dict{Char,Float64}[Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}(), Dict{Char,Float64}()])
```

Reference:
V’yugin, Vladimir V. "Online Aggregation of Unbounded Signed Losses Using Shifting Experts." Conformal and Probabilistic Prediction and Applications. 2017.
"""
mutable struct AdaptiveMPP{T} <: BasePredictor{T}
    model::Trie{T,Int}
    context::Vector{T}
    cxt_length::Int
    α::Float64
    η::Float64
    Δ::Float64
    weights::Vector{Float64}
    last_prediction::Vector{Dict{T,Float64}}

    AdaptiveMPP{T}( _c::Int, _α = 0.10 ) where {T} = new( Trie{T,Int}(), Vector{T}(),
                                                _c, _α, +Inf, 0.00,
                                                1/(1+_c)*ones(Float64,_c+1), fill(Dict{T,Float64}(),_c+1) )
end

function add!( p::AdaptiveMPP{T}, sym::T ) where {T}
    p.model.value += 1;

    # Calculate expert loss
    expert_loss = zeros(Float64,p.cxt_length+1);
    for expert_idx = 1:p.cxt_length+1
        expert_best_symbol      = get_best_symbol( p.last_prediction[expert_idx] )
        expert_loss[expert_idx] = (expert_best_symbol==sym) ? 0 : 1
    end

    # Compute the aggregating algorithm loss
    h_t = sum( p.weights .* expert_loss )

    # Make loss update
    _w  = map( w->isnan(w) ? 1/(p.cxt_length+1) : w, p.weights.*exp.(-p.η*expert_loss) )
    _w  = _w ./ sum(_w)

    # Mixing Update
    # Based on example-1 with α_{t+1} = 0.1 ∀ t
    p.weights   = p.α * 1/(p.cxt_length) .+ (1-p.α) * _w

    # Learning Parameter Update
    m_t = -1/p.η * log( sum( p.weights .* exp.(-p.η*expert_loss) ) )
    δ_t = h_t - m_t
    p.Δ = p.Δ + δ_t
    p.η = 1 / p.Δ

    # Create node path
    push!( p.context, sym )
    buffer  = p.context[1:end]
    while !isempty(buffer)
        p.model[buffer] = (haskey(p.model,buffer)) ? p.model[buffer]+1 : 1
        popfirst!( buffer )
    end

    if length(p.context) > p.cxt_length
        popfirst!( p.context )
    end

    nothing
end

function predict( p::AdaptiveMPP{T} ) where {T}
    # Clear away last predictions
    p.last_prediction   = fill(Dict{T,Float64}(),p.cxt_length+1)

    # Calculate normalized weights
    norm_weights    = p.weights ./ sum(p.weights)
    # println( "    norm_weights: ", norm_weights )

    # Get prediction from 0-order expert
    symbols = Dict( k => (p.model[[k]]/p.model.value) for k ∈ keys(children(p.model,Vector{T}())) )

    # Save this for later
    p.last_prediction[1]    = deepcopy( symbols )

    # Apply 0-order expert weight
    for sym ∈ keys(symbols)
        symbols[sym] *= norm_weights[1]
    end
    # println( "After Order-0: ", symbols )

    # Do prediction from all other higher orders
    no_of_active_experts    = min( p.model.value, p.cxt_length )
    for expert_idx  = 1:no_of_active_experts
        # Predict
        this_expert_prediction = predict_from_subcontext( p, p.context[1:expert_idx] )
        # Normalize probability out from experts - jugaad
        denom   = 0;
        for (sym,prob) ∈ this_expert_prediction denom += prob; end
        for sym ∈ keys(this_expert_prediction) this_expert_prediction[sym] /= denom;    end

        # Save
        p.last_prediction[expert_idx+1] = deepcopy( this_expert_prediction )
        # Apply weight to expert and add to final prediction
        for sym ∈ keys( symbols )
            symbols[sym] += (this_expert_prediction[sym]*norm_weights[expert_idx+1])
        end
        # println( "After order-$expert_idx: ", symbols );
    end

    return symbols;
end

function info_string( p::AdaptiveMPP{T} ) where {T}
    return @sprintf( "Adaptive MPP(%d,\$\\alpha_t=%3.2f\$)", p.cxt_length,p.α )
end

function unique_string( p::AdaptiveMPP{T} ) where {T}
    return @sprintf( "adaptiveMPP_%02d_%03d", p.cxt_length, trunc(Int,100*p.α) )
end

function predict_from_subcontext( p::AdaptiveMPP{T}, sub_cxt::Vector{T} ) where {T}
    # Create a dictionary with symbols
    symbols = Dict( k => (p.model[[k]]/p.model.value) for k ∈ keys(children(p.model,Vector{T}())) )

    buffer  = Vector{T}()
    for i = length(sub_cxt):-1:1
        pushfirst!( buffer, sub_cxt[i] )
        # println( "    Buffer: ", buffer, " -> ", p.model[buffer] );
        # Apply escape probability
        esc_prob = 1/p.model[buffer]
        for symbol in keys(symbols)
            symbols[symbol] *= esc_prob
        end
        # Get each child probability
        list_of_children    = children( p.model, buffer )
        for k ∈ keys(list_of_children)
            # println( "        ", k, " -> ", list_of_children[k].value )
            symbols[k] += list_of_children[k].value/p.model[buffer]
        end
    end

    return symbols;
end


text:
module ContinuumAFS

using RCall
using RvSpectMLBase
using EchelleInstruments

function calc_continuum_model_order(spectrum::AbstractSpectra2D; order_idx::Integer )
    possible_pix = get_inst_module(get_inst(spectrum)).get_pixel_range(get_inst(spectrum),order_idx)
    bad_pix = get_inst_module(get_inst(spectrum)).bad_col_ranges(get_inst(spectrum),order_idx)
    pix_rng = EchelleInstruments.calc_complement_index_ranges(possible_pix,bad_pix)
    pix = mapreduce(p->collect(p),vcat,pix_rng)
    afs_inputs = zeros(Float64,length(pix),2)
    afs_inputs[:,1] .= spectrum.λ[pix,order_idx]
    afs_inputs[:,2] .= spectrum.flux[pix,order_idx]
    @assert !any(isnan.(afs_inputs))
    #=
    wv = mapreduce(p->spec.λ[p,order_idx],vcat,pix_rng)
    @assert !any(isnan.(wv))
    inten = mapreduce(p->convert(Vector{Float64},spec.flux[p,order_idx]),vcat,pix_rng)
    @assert !any(isnan.(inten))
    afs_inputs = hcat(wv,inten)
    =#
    #df = DataFrame("wv"=>wv,"intes"=>inten)
    afs_src = joinpath("src","AFS.R")
    R"source($afs_src)"
    afs_output_R = R"AFS($afs_inputs,0.95,0.25)"
    afs_output = rcopy(afs_output_R)
    continuum = zeros(eltype(spectrum.flux),size(spectrum.flux[:,order_idx]))
    continuum = fill(NaN,size(spectrum.flux[:,order_idx]))
    continuum[pix] .= afs_output
    return continuum
end

function calc_continuum_model(spectrum::AbstractSpectra2D )
    vec_of_orders = map(ord->calc_continuum_model_order(spectrum,order_idx=ord), min_order(get_inst(spectrum)):max_order(get_inst(spectrum)) )
    output = fill(NaN, size(spectrum.flux))
    for (i,ord) in enumerate(min_order(get_inst(spectrum)):max_order(get_inst(spectrum)))
        output[:,ord] .= vec_of_orders[i]
    end
    return output
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2e0f805a-68d3-11eb-35c6-95c237883fcd
using PlutoUI

# ╔═╡ e5eee32e-68de-11eb-3267-33a033f1dab1
begin
	using Printf

	macro timeit(expr)
		quote
			local t₀ = time()         ## avoid poluting client env.
			local res = $(esc(expr))  ## eval. expr
			local dt = time() - t₀
			print("Ellapsed time: ")
			@printf("%.9f\n", dt)

			res                       ## finally return res
		end
	end
end

# ╔═╡ f6867c32-68df-11eb-042e-1bd6b900f6fa
using Test

# ╔═╡ dc528154-68d2-11eb-2f60-f3499c34b532
md"""
## Macros in Julia
"""

# ╔═╡ eb0c80ff-9b03-4391-a2b1-ff56e3e4c8a9
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ f88b338e-68d2-11eb-1fa8-0bf5b252e9b6
md"""
### Example: assert macro
"""

# ╔═╡ 04400be6-68d3-11eb-35b3-79c07b358527
macro assert(expr)
    :($expr ? nothing : error("Assertion failed: ", $(string(expr))))
end

# ╔═╡ 3f24e7a4-68d3-11eb-22b9-f70cbdfa1b76
md"""
### Assert with message
"""

# ╔═╡ 51b16460-68d3-11eb-0b52-1b31d835485a
macro assert(expr, msg)
    quote
        if $expr
            nothing
        else
            local m = length($msg) == 0 ? "Assertion failed: " : string($msg, ' ')
            error(m, $(string(expr)))
        end
    end
end

# ╔═╡ 1148416e-68d3-11eb-24b3-b1a15aa812a0
begin
	@assert 1 == 1.0

	@assert 1 != 42.0
end

# ╔═╡ 282df856-68d3-11eb-3424-23d08b7af6a3
@assert 1 == 0.0

# ╔═╡ 66f4760a-68d3-11eb-1426-871a3354d9dd
@assert 1 == 42.0

# ╔═╡ 736424d0-68d3-11eb-1ada-750f7c8d622c
@assert(1 == 42.0, "How come?")

# ╔═╡ 90ab2322-68d3-11eb-21f0-35e49b201db4
@assert 1 == 42.0 "How come?"   ## no parentheses, no comma

# ╔═╡ ba35d962-68d3-11eb-1e10-05a84354a399
md"""
### Macro timeit
"""

# ╔═╡ db90b314-68df-11eb-02d5-69315dac74ad
function fact(n::Integer)::Integer
	n < 0 && throw(ArgumentError("n should be positive"))
    (n == 0 || n == 1) && return 1
    n * fact(n - 1) # remember no tail recursion!
end

# ╔═╡ fdcc1536-68df-11eb-34e1-55281df41ceb
begin
	@test fact(0) == 1
	@test fact(1) == 1
	@test_throws ArgumentError fact(-2)
	@test fact(2) == 2
	@test fact(3) == 6
	@test fact(4) == 24
end

# ╔═╡ 8237cd88-68e0-11eb-111c-3b0c6245d629
with_terminal() do
	@timeit fact(7)
end

# ╔═╡ c7e510ca-68e0-11eb-3bfe-d338d10081aa
with_terminal() do
	@timeit fact(15)
end

# ╔═╡ 0923f20e-68e1-11eb-37ae-19e4a24f48b4
md"""
### Macro Trace
"""

# ╔═╡ 5b0d1a14-68e1-11eb-2068-b5b4249c0d61
macro trace(expr)
	quote
		local fn_name = $(string(expr.args[1]))
		local args = $(string(expr.args[2]))

		println("Start of exec. function ", string(fn_name), " with args: ", args)
		local res = $(esc(expr))    ## call the actual expr ≡ function
		println("End of exec. function ", string(fn_name))

		res
	end
end

# ╔═╡ 3dab143e-68e2-11eb-100c-016893824ccd
begin
	fn_α = function (α::Integer=2)  ## function builder
    	## this is a closure
    	λ(x::Integer) = x * α
	end

	double_fn = fn_α()
	triple = fn_α(3)
end

# ╔═╡ e5397084-68e1-11eb-3322-5956fa8fd17a
# let see the content of this macro, using another macro
@macroexpand @trace double_fn(2)

# ╔═╡ 86936d9a-68e2-11eb-2217-47e57c631638
begin
	expr = :(double_fn(3))

	expr.args
end

# ╔═╡ a60c13b6-68e2-11eb-3fab-f73606ea715d
@test triple(11) == 33

# ╔═╡ b4e68d26-68e2-11eb-3bf9-bbba25301549
with_terminal() do
	@trace triple(11)
end

# ╔═╡ df8b48a8-68e2-11eb-1698-8f414b4a7fc4
md"""
### Macro unless
"""

# ╔═╡ ebe12cb4-68e2-11eb-216e-a3869afd2b1f
macro unless0(cond, body)
	quote
		if !$(esc(cond))
			$(esc(body))
		end
	end
end

# ╔═╡ 2599be12-68e3-11eb-07c4-376988e519d0
with_terminal() do
	x = 2

	@unless0 x < 0 begin
		println("we are done with x: $(x)")
	end
end

# ╔═╡ 6c19c67a-68e3-11eb-1c71-21ffa0892f7d
with_terminal() do
	x = -1

	@unless0 x < 0 begin
		println("we are done with x: $(x)")
	end
end

# Expect nothing

# ╔═╡ 78acc86c-68e3-11eb-2048-0969c36a0f1a
macro unless(cond, body)
	:(!$(esc(cond)) && $(esc(body)))

	# ≡ @eval $expr ? nothing : $body
end

# ╔═╡ a3511828-68e3-11eb-07ec-4de4d1836d95
with_terminal() do
	x = 2

	@unless x < 0 begin
		println("we are done with x: $(x)")
	end
end

# ╔═╡ a5bdc6f6-68e3-11eb-2248-fb359467846d
md"""
### Expression
"""

# ╔═╡ c3b46a5c-68e3-11eb-0edd-9f98a253f055
with_terminal() do
	dump(:(2 + a * b - c))
end

# ╔═╡ dee8a4f0-68e3-11eb-2336-abcb7af27dcf
with_terminal() do
	dump(:(-c + 2 - a * b))
end

# ╔═╡ a87243cb-3c70-4ecd-ad6e-80a876d6d261
html"""
<style>
  main {
        max-width: calc(800px + 25px + 6px);
  }
  .plutoui-toc.aside {
    background: linen;
  }
  h3, h4 {
        background: wheat;
  }
</style>
"""

# ╔═╡ Cell order:
# ╟─dc528154-68d2-11eb-2f60-f3499c34b532
# ╠═2e0f805a-68d3-11eb-35c6-95c237883fcd
# ╟─eb0c80ff-9b03-4391-a2b1-ff56e3e4c8a9
# ╟─f88b338e-68d2-11eb-1fa8-0bf5b252e9b6
# ╠═04400be6-68d3-11eb-35b3-79c07b358527
# ╠═1148416e-68d3-11eb-24b3-b1a15aa812a0
# ╠═282df856-68d3-11eb-3424-23d08b7af6a3
# ╠═3f24e7a4-68d3-11eb-22b9-f70cbdfa1b76
# ╠═51b16460-68d3-11eb-0b52-1b31d835485a
# ╠═66f4760a-68d3-11eb-1426-871a3354d9dd
# ╠═736424d0-68d3-11eb-1ada-750f7c8d622c
# ╠═90ab2322-68d3-11eb-21f0-35e49b201db4
# ╟─ba35d962-68d3-11eb-1e10-05a84354a399
# ╠═e5eee32e-68de-11eb-3267-33a033f1dab1
# ╠═db90b314-68df-11eb-02d5-69315dac74ad
# ╠═f6867c32-68df-11eb-042e-1bd6b900f6fa
# ╠═fdcc1536-68df-11eb-34e1-55281df41ceb
# ╠═8237cd88-68e0-11eb-111c-3b0c6245d629
# ╠═c7e510ca-68e0-11eb-3bfe-d338d10081aa
# ╟─0923f20e-68e1-11eb-37ae-19e4a24f48b4
# ╠═5b0d1a14-68e1-11eb-2068-b5b4249c0d61
# ╠═e5397084-68e1-11eb-3322-5956fa8fd17a
# ╠═3dab143e-68e2-11eb-100c-016893824ccd
# ╠═86936d9a-68e2-11eb-2217-47e57c631638
# ╠═a60c13b6-68e2-11eb-3fab-f73606ea715d
# ╠═b4e68d26-68e2-11eb-3bf9-bbba25301549
# ╟─df8b48a8-68e2-11eb-1698-8f414b4a7fc4
# ╠═ebe12cb4-68e2-11eb-216e-a3869afd2b1f
# ╠═2599be12-68e3-11eb-07c4-376988e519d0
# ╠═6c19c67a-68e3-11eb-1c71-21ffa0892f7d
# ╠═78acc86c-68e3-11eb-2048-0969c36a0f1a
# ╠═a3511828-68e3-11eb-07ec-4de4d1836d95
# ╟─a5bdc6f6-68e3-11eb-2248-fb359467846d
# ╠═c3b46a5c-68e3-11eb-0edd-9f98a253f055
# ╠═dee8a4f0-68e3-11eb-2336-abcb7af27dcf
# ╟─a87243cb-3c70-4ecd-ad6e-80a876d6d261

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Qt5Charts_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Qt5Charts")
JLLWrappers.@generate_main_file("Qt5Charts", UUID("dd720b4e-75c8-5196-993d-eac563881c8e"))
end  # module Qt5Charts_jll


text:
module IoInterfacesTests

end # module


text:
###################
# ForwardOptimize #
###################

for (M, f, arity) in DiffRules.diffrules()
    if arity == 1
        @eval @inline $M.$(f)(t::TrackedReal) = ForwardOptimize($M.$(f))(t)
    elseif arity == 2
        @eval @inline $M.$(f)(a::TrackedReal, b::TrackedReal) = ForwardOptimize($M.$(f))(a, b)
        for R in REAL_TYPES
            @eval begin
                @inline $M.$(f)(a::TrackedReal, b::$R) = ForwardOptimize($M.$(f))(a, b)
                @inline $M.$(f)(a::$R, b::TrackedReal) = ForwardOptimize($M.$(f))(a, b)
            end
        end
    end
end

################
# SkipOptimize #
################

# unary #
#-------#

for f in SKIPPED_UNARY_SCALAR_FUNCS
    @eval @inline Base.$(f)(t::TrackedReal) = SkipOptimize($(f))(t)
end

# binary #
#--------#

for f in SKIPPED_BINARY_SCALAR_FUNCS
    @eval @inline Base.$(f)(a::TrackedReal, b::TrackedReal) = SkipOptimize($(f))(a, b)
    for R in REAL_TYPES
        @eval begin
            @inline Base.$(f)(a::$R, b::TrackedReal) = SkipOptimize($(f))(a, b)
            @inline Base.$(f)(a::TrackedReal, b::$R) = SkipOptimize($(f))(a, b)
        end
    end
end

###########
# reverse #
###########

@noinline function scalar_reverse_exec!(instruction::ScalarInstruction{F,I,O,C}) where {F,I,O,C}
    f = instruction.func
    input = instruction.input
    output = instruction.output
    partials = instruction.cache[]
    if istracked(input)
        increment_deriv!(input, deriv(output) * partials)
    else
        a, b = input
        output_deriv = deriv(output)
        a_partial, b_partial = partials
        if istracked(a) && istracked(b)
            increment_deriv!(a, output_deriv * a_partial)
            increment_deriv!(b, output_deriv * b_partial)
        elseif istracked(a)
            increment_deriv!(a, output_deriv * a_partial)
        else
            increment_deriv!(b, output_deriv * b_partial)
        end
    end
    unseed!(output)
    return nothing
end

###########
# forward #
###########

@noinline function scalar_forward_exec!(instruction::ScalarInstruction{F,I,O,C}) where {F,I,O,C}
    f = instruction.func
    input = instruction.input
    output = instruction.output
    cache = instruction.cache
    if istracked(input)
        unary_scalar_forward_exec!(f, output, input, cache)
    else
        binary_scalar_forward_exec!(f, output, input, cache)
    end
    return nothing
end

@noinline function unary_scalar_forward_exec!(f::F, output::O, input, cache) where {F,O}
    pull_value!(input)
    result1 = DiffResult(zero(valtype(O)), zero(valtype(O)))
    result1 = ForwardDiff.derivative!(result1, f, value(input))
    value!(output, DiffResults.value(result1))
    cache[] = DiffResults.derivative(result1)
    return nothing
end

@noinline function binary_scalar_forward_exec!(f::F, output::O, input, cache) where {F,O}
    a, b = input
    pull_value!(a)
    pull_value!(b)
    if istracked(a) && istracked(b)
        result2 = DiffResults.GradientResult(SVector(zero(valtype(O)), zero(valtype(O))))
        result2 = ForwardDiff.gradient!(result2, x -> f(x[1], x[2]), SVector(value(a), value(b)))
        value!(output, DiffResults.value(result2))
        cache[] = DiffResults.gradient(result2)
    else
        result1 = DiffResult(zero(valtype(O)), zero(valtype(O)))
        if istracked(a)
            result1 = ForwardDiff.derivative!(result1, va -> f(va, b), value(a))
        else
            result1 = ForwardDiff.derivative!(result1, vb -> f(a, vb), value(b))
        end
        value!(output, DiffResults.value(result1))
        partial = DiffResults.derivative(result1)
        cache[] = SVector(partial, partial)
    end
    return nothing
end


text:
module Ex5_2
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.2 - EFDP has no solution
println("Example 5.2")
# define s as an improper transfer function
s = rtf('s');
# define Gu(s), Gd(s), Gf(s)
Gu = [(s+1)/(s+2); (s+2)/(s+3)];     # enter Gu(s)
Gd = [1/(s+2); 0];                   # enter Gd(s)
Gf = [(s+1)/(s+2) 0; 0 1];           # enter Gf(s)
p = 2; mu = 1; md = 1; mf = 2;       # set dimensions

# build model with additive faults
#sysf = fdimodset(ss([Gu Gd Gf]),struct('c',1:mu,'d',mu+(1:md),'f',mu+md+(1:mf)));
sysf = fdimodset(dss([Gu Gd Gf]),c = 1:mu, d = mu.+(1:md), f = (mu+md).+(1:mf))

# call of EFDSYN with default options
try
  Q, Rf, info = efdsyn(sysf)
  @test false
catch err
  @test true
  display(err)
  # compute achievable structure matrix
  println("S_achievable = $(fdigenspec(sysf))")
end

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
"""
	PerformanceIndicators

This module includes performance indicators to assess evolutionary multi-objective
optimization algorithms.

- `gd` Generational Distance.
- `igd` Inverted Generational Distance.
- `gd_plus` Generational Distance plus.
- `igd_plus` Inverted Generational Distance plus.
- `covering` Covering indicator (C-metric).
- `hypervolume` Hypervolume indicator.

### Example

```jldoctest
julia> import Metaheuristics: PerformanceIndicators, TestProblems

julia> A = [ collect(1:3) collect(1:3) ]
3×2 Array{Int64,2}:
 1  1
 2  2
 3  3

julia> B = A .- 1
3×2 Array{Int64,2}:
 0  0
 1  1
 2  2

julia> PerformanceIndicators.gd(A, B)
0.47140452079103173

julia> f, bounds, front = TestProblems.get_problem(:ZDT1);

julia> front
                          F space
         ┌────────────────────────────────────────┐
       1 │⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠈⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⠀⠉⠢⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   f_2   │⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠲⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⢤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠙⠢⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⠢⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠢⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠑⠢⢤⣀⠀⠀⠀⠀⠀│
       0 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠢⢄⣀│
         └────────────────────────────────────────┘
         0                                        1
                            f_1

julia> PerformanceIndicators.igd_plus(front, front)
0.0
```

"""
module PerformanceIndicators

import ..State, ..fvals, ..norm, ..xFgh_indiv, ..fval, ..compare, ..mean
import ..get_non_dominated_solutions

include("hypervolume.jl")

"""
	generational_distance(front, true_pareto_front; p = 1, inverted = false, plus=false)

Returns the Generational Distance.

### Parameters

- `front` and `true_pareto_front` can be:
	- `N×m` matrix where `N` is the number of points and `m` is the number of objectives.
	- `State`
	- `Array{xFgh_indiv}` (usually `State.population`)
- `true_pareto_front` is a `M×m` matrix.
- `p` is the power in for the ‖⋅‖_p
- `plus` if true then computes the GD+
- `inverted` if true then computes IGD
"""
function generational_distance(front, true_pareto_front; p = 1, inverted = false, plus=false)

	if inverted
		tmp = true_pareto_front
		true_pareto_front = front
		front = tmp
	end

	distances = zeros(size(true_pareto_front, 1))

	s = 0.0
	for i in 1:size(front,1)
		x = view(front, i, :)
		for j = 1:size(true_pareto_front, 1)
			y = view(true_pareto_front, j, :)
			if plus	&& inverted
				distances[j] = norm( max.(y - x, 0) )
			elseif plus && !inverted
				distances[j] = norm( max.(x - y, 0) )
			else
				distances[j] = norm( y - x)
			end

		end

		s += minimum(distances) ^ p
	end

	return (s^(1/p)) /  size(front, 1)

end

function generational_distance(front::Union{Array{xFgh_indiv}, State},
						       true_pareto_front::Union{Array{xFgh_indiv}, State};
								p = 1,
								inverted = false,
								plus = false
		)
	generational_distance(fvals(front), fvals(true_pareto_front); p=p,inverted=inverted,plus=plus)
end

function generational_distance(front::Union{Array{xFgh_indiv}, State},
						       true_pareto_front;
								p = 1,
								inverted = false,
								plus = false
		)
	generational_distance(fvals(front), (true_pareto_front); p=p,inverted=inverted,plus=plus)
end

"""
	gd(front, true_pareto_front; p = 1)

Returns the Generational Distance.

### Parameters

`front` and `true_pareto_front` can be:

- `N×m` matrix where `N` is the number of points and `m` is the number of objectives.
- `State`
- `Array{xFgh_indiv}` (usually `State.population`)
"""
gd(front, true_pareto_front; p = 1) = generational_distance(front, true_pareto_front, p=p)

"""
	igd(front, true_pareto_front; p = 1)

Returns the Inverted Generational Distance.

### Parameters

`front` and `true_pareto_front` can be:

- `N×m` matrix where `N` is the number of points and `m` is the number of objectives.
- `State`
- `Array{xFgh_indiv}` (usually `State.population`)

"""
igd(front, true_pareto_front; p = 1) = generational_distance(front, true_pareto_front, inverted=true, p=p)

"""
	gd_plus(front, true_pareto_front; p = 1)

Returns the Generational Distance Plus.

### Parameters

`front` and `true_pareto_front` can be:

- `N×m` matrix where `N` is the number of points and `m` is the number of objectives.
- `State`
- `Array{xFgh_indiv}` (usually `State.population`)
"""
gd_plus(front, true_pareto_front; p = 1) = generational_distance(front, true_pareto_front, p=p, plus=true)

"""
	igd_plus(front, true_pareto_front; p = 1)

Returns the Inverted Generational Distance Plus.

### Parameters

`front` and `true_pareto_front` can be:

- `N×m` matrix where `N` is the number of points and `m` is the number of objectives.
- `State`
- `Array{xFgh_indiv}` (usually `State.population`)

"""
igd_plus(front, true_pareto_front; p = 1) = generational_distance(front, true_pareto_front, inverted=true, p=p, plus=true)

"""
	spacing(A)
Computes the Schott spacing indicator. `spacing(A) == 0` means that vectors in `A` are
uniformly distributed.
"""
function spacing(A::Array{Array{Float64,1},1})
    n = length(A)

    d = zeros(n)
    for i = 1:n
    	@inbounds d[i] = minimum( [ sum( abs.(A[i] - A[j]) ) for j in deleteat!(collect(1:n), i) ] )
    end

    d_hat = mean(d)
    s = (1.0 / ( n - 1.0 )) * sum( (d_hat .- d) .^ 2 )
    return sqrt(s)

end

spacing(A::Array{xFgh_indiv}) = spacing( fval.(A) )
spacing(A::State{xFgh_indiv}) = spacing(A.population)
spacing(A::Matrix) = spacing([A[i,:] for i in 1:size(A,1)])

"""
	covering(A, B)
Computes the covering indicator (percentage of vectors in B that are dominated by vectors in A)
from two sets with non-dominated solutions.

A and B with size (n, m) where n is number of samples and m is the vector dimension.

Note that `covering(A, B) == 1` means that all solutions in B are dominated by those
in A. Moreover, `covering(A, B) != covering(B, A)` in general.

If `A::State` and `B::State`, then computes `covering(A.population, B.population)` after
ignoring dominated solutions in each set.
"""
function covering(A::Array{Vector{T}}, B::Array{Vector{T}}) where T <: Real
    n_a = length(A)
    n_b = length(B)
    s = 0.0
    for i = 1:n_b
        for j = 1:n_a
			# `compare` returns 1 if argument 1 dominates argument 2
            if compare(A[j], B[i]) == 1
                s += 1.0
                break
            end
        end

    end

    return s / n_b

end

covering(A::Vector{xFgh_indiv}, B::Vector{xFgh_indiv}) = covering(fval.(A), fval.(B))

function covering(A::State{xFgh_indiv}, B::State{xFgh_indiv})
	A_non_dominated = get_non_dominated_solutions(A.population)
	B_non_dominated = get_non_dominated_solutions(B.population)

	length(A_non_dominated) != length(A.population) &&
	@warn "Some solutions in B dominate each other (solutions ignored)."

	length(B_non_dominated) != length(B.population) &&
	@warn "Some solutions in B dominate each other (solutions ignored)."

	covering(A_non_dominated, B_non_dominated)
end

function covering(A::Matrix, B::Matrix)
	A_arr = [ A[i,:] for i in 1:size(A,1) ]
	B_arr = [ B[i,:] for i in 1:size(B,1) ]
	covering(A_arr, B_arr)
end

"""
	hypervolume(front, reference_point)
Computes the hypervolume indicator, i.e., volume between points in `front` and `reference_point`.

Note that each point in `front` must (weakly) dominates to `reference_point`. Also, `front`
is a non-dominated set.

If `front::State` and `reference_point::Vector`, then computes `hypervolume(front.population, reference_point)` after
ignoring solutions in `front` that do not dominate `reference_point`.
"""
function hypervolume(front::Array{Vector{T}}, reference_point::Vector) where T <: Real

	weaklyDominates(point, other) = begin
		for i in 1:length(point)
			if point[i] > other[i]
				return false
			end
		end
		return true
	end

	relevantPoints = Vector[]
	for point in front
		# only consider points that dominate the reference point
		if weaklyDominates(point, reference_point)
			push!(relevantPoints, point)
		end
	end

	if length(relevantPoints) != length(front)
		ign = length(front) - length(relevantPoints)
		rel = length(relevantPoints)
		@warn "Ignoring $ign points dominated by the reference point ($rel points are used)."
	end

	return HyperVolume.hv(relevantPoints, reference_point)
end

hypervolume(front::Vector{xFgh_indiv}, reference_point::Vector{xFgh_indiv}) = hypervolume(fval.(front), fval(reference_point))
hypervolume(front::Vector{xFgh_indiv}, reference_point::Vector) = hypervolume(fval.(front), reference_point)

function hypervolume(front::Matrix, reference_point::Vector)
	front_ = [ front[i,:] for i in 1:size(front,1) ]
	hypervolume(front_, reference_point)
end

"""
	deltap(front, true_pareto_front; p = 1)
	Δₚ(front, true_pareto_front; p = 1)

Returns the averaged Hausdorff distance indicator aka Δₚ (Delta p).

"`Δₚ`" can be typed as `\\Delta<tab>\\_p<tab>`.

### Parameters

`front` and `true_pareto_front` can be:

- `N×m` matrix where `N` is the number of points and `m` is the number of objectives.
- `Array{xFgh_indiv}` (usually `State.population`)

"""
function deltap(front::Array{Vector{T}}, true_pareto_front::Array{Vector{T}}; p = 1) where T <: Real
	distances = zeros(length(true_pareto_front), length(front))

	for i in eachindex(true_pareto_front)
		for j in eachindex(front)
			distances[i,j] = norm(true_pareto_front[i] - front[j])
		end
	end

	gd  = mean(minimum(distances, dims = 1) .^ p)
	igd = mean(minimum(distances, dims = 2) .^ p)
	return max(gd^(1/p), igd^(1/p))
end

deltap(front::Vector{xFgh_indiv}, true_pareto_front::Vector{xFgh_indiv}; kargs...) = deltap(fval.(front), fval.(true_pareto_front); kargs...)

function deltap(front::AbstractMatrix, true_pareto_front::AbstractMatrix; kargs...)
	front_ = [ front[i,:] for i in 1:size(front,1) ]
	true_pareto_front_ = [ true_pareto_front[i,:] for i in 1:size(true_pareto_front,1) ]
	deltap(front_, true_pareto_front_; kargs...)
end

const Δₚ = deltap
#Δₚ(args...;kargs...) = deltap(args...;kargs...)

end


text:
type SteadyStateSolution{uType,R,P,A} <: AbstractSteadyStateSolution
  u::uType
  resid::R
  prob::P
  alg::A
  retcode::Symbol
end

function build_solution(prob::AbstractSteadyStateProblem,
                        alg,u,resid;calculate_error = true,
                        retcode = :Default, kwargs...)
  SteadyStateSolution(u,resid,prob,alg,retcode)
end


text:
type Run{F <: AbstractFloat, I <: Integer, C <: Configuration}
    cost::Function
    cost_arguments::Dict{Symbol, Any}
    cost_evaluations::I
    cost_values::Array{F, 1}
    starting_point::C
    starting_cost::F
    report_after::I
    reporting_criterion::Function
    measurement_method::Function
    stopping_criterion::Function
    duration::I
    methods::Array{Any, 2}
    channel_size::I

    function call{F <: AbstractFloat, I <: Integer, C <: Configuration}(::Type{Run};
            cost::Function                    = (c(x) = 0),
            cost_arguments::Dict{Symbol, Any} = Dict{Symbol, Any}(),
            cost_evaluations::I               = 1,
            cost_values::Array{F, 1}          = [0.0],
            starting_point::C                 = Configuration("empty"),
            starting_cost::F                  = 0.0,
            report_after::I                   = 100,
            reporting_criterion::Function     = iterations_reporting_criterion,
            measurement_method::Function      = measure_mean!,
            stopping_criterion::Function      = iterations_criterion,
            duration::I                       = 1_000,
            methods::Array{Any, 2}            = [[:simulated_annealing 2];],
            channel_size::I                   = 4096)

        new{F, I, C}(cost,
                     cost_arguments,
                     cost_evaluations,
                     cost_values,
                     starting_point,
                     starting_cost,
                     report_after,
                     reporting_criterion,
                     measurement_method,
                     stopping_criterion,
                     duration,
                     methods,
                     channel_size)
    end
end


text:
# Functions supporting trendline identification (support/resistance, zigzag, elliot waves, etc.)

"""
```
maxima(x::Array{T}; threshold::T=0.0, order::Int=1) where {T<:Real}
```

Estimate local maxima of a time series
"""
function maxima(x::Array{T}; threshold::T=0.0, order::Int=1) where {T<:Real}
    @assert threshold >= 0.0 "threshold must be positive"
    @assert order > 0 "order must be a positive integer"
    n = size(x,1)
    crit = falses(n)
    @inbounds for i=2:n-1
        if (x[i]-x[i-1] >= threshold) && (x[i]-x[i+1] >= threshold)
            crit[i] = true
        end
    end
    while order > 1
        idx = findall(crit)
        crit[idx[!maxima(x[crit], threshold=threshold)]] .= false
        order -= 1
    end
    return crit
end

"""
```
minima(x::Array{T}; threshold::T=0.0, order::Int=1) where {T<:Real}
```

Estimate local minima of a time series
"""
function minima(x::Array{T}; threshold::T=0.0, order::Int=1) where {T<:Real}
    @assert threshold <= 0.0 "threshold must be negative"
    @assert order > 0 "order must be a positive integer"
    n = size(x,1)
    crit = falses(n)
    @inbounds for i=2:n-1
        if (x[i]-x[i-1] <= threshold) && (x[i]-x[i+1] <= threshold)
            crit[i] = true
        end
    end
    while order > 1
        idx = findall(crit)
        crit[idx[!minima(x[crit], threshold=threshold)]] .= false
        order -= 1
    end
    return crit
end

function interpolate(x1::Int, x2::Int, y1::T, y2::T) where {T<:Real}
	m = (y2-y1)/(x2-x1)
	b = y1 - m*x1
	x = collect(x1:1.0:x2)
	y = m*x .+ b
	return y
end

"""
```
resistance(x::Array{T}; order::Int=1, threshold::T=0.0) where {T<:Real}
```

Estimate resistance lines of a financial time series
"""
function resistance(x::Array{T}; order::Int=1, threshold::T=0.0) where {T<:Real}
    out = zeros(size(x))
    crit = maxima(x, threshold=threshold, order=order)
    out[.!crit] .= NaN
    idx = findall(crit)
    @inbounds for i=2:length(idx)
        out[idx[i-1]:idx[i]] .= interpolate(idx[i-1], idx[i], x[i-1], x[i])
    end
    return out
end

"""
```
support(x::Array{T}; order::Int=1, threshold::T=0.0) where {T<:Real}
```

Estimate support lines of a financial time series
"""
function support(x::Array{T}; order::Int=1, threshold::T=0.0) where {T<:Real}
    out = zeros(size(x))
    crit = minima(x, threshold=threshold, order=order)
    out[.!crit] .= NaN
    idx = findall(crit)
    @inbounds for i=2:length(idx)
        out[idx[i-1]:idx[i]] .= interpolate(idx[i-1], idx[i], x[i-1], x[i])
    end
    return out
end