text:
export solve, Result, nresults, nsolutions, nfinite, nsingular, nat_infinity,
    nfailed, nnonsingular, nreal, ntracked, finite, results, mapresults,
    failed, at_infinity, singular, nonsingular, seed,
    solutions, real_solutions, multiplicities, statistics, multiplicities!

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

The solve function takes many different arguments and options depending on your specific situation,
but in the it always returns a [`Result`](@ref) containing the result of the computations.
In the following we show the different inputs `solve` takes and at the end we list all possible options.

# 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 `d1⋯dn`, where `di` is the degree of the `i`th entry of `F`. `F` can be
- `Vector{<:MultivariatePolynomials.AbstractPolynomial}` (e.g. constructed by `@polyvar`)
- A composition of polynomial systems constructed by [`compose`](@ref).
- [`AbstractSystem`](@ref) (the system has to represent a **homogeneous** polynomial system.)

### Example
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` using a start system computed from the Newton Polytopes of the entries `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).
- [`AbstractSystem`](@ref) (the system has to represent a **homogeneous** polynomial system.)

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

### Example
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
```

# Homogeneous Systems

If `F` has is homogeneous, we return results in projective space

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

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])`.

# Multihomogeneous Systems

By exploiting the multi-homogenous 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.

# 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`.

### Example
```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]])
```

# 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, (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 `γ₁` and `γ₀` are ``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`.

### Example
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₀)
```

# 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``.
Note that `H` has to be a homotopy between *homogeneous* polynomial systems.
If it should be considered as an affine system indicate which is the index
of the homogenization variable, e.g. `solve(H, startsolutions, homvar=3)`
if the third variable is the homogenization variable.

# Options
### General options:

* `seed::Int`: The random seed used during the computations.
* `show_progress=true`: Whether a progress bar should be printed to standard out.
* `threading=true`: Enable or disable multi-threading.
* `path_result_details=:default`: The amount of information computed in each path result. Possible values are `:minimal` (minimal details), `:default` (default) and `:extensive` (all information possible).
* `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`.
* `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`.
* `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.
* `affine_tracking::Bool=true`: Indicate whether path tracking should happen in affine space.
* `projective_tracking::Bool=false`: Indicate whether path tracking should happen in projective space. The flag `affine_tracking` is dominant.
* `path_jumping_check::Bool=true`: Enable a check whether one of the paths jumped to another one.

### Path tracking specific options:
* `accuracy=1e-7`: The accuracy required during the path tracking.
* `corrector::AbstractCorrector`: The corrector used during in the predictor-corrector scheme. The default is [`NewtonCorrector`](@ref).
* `initial_step_size=0.1`: The size of the first step.
* `max_corrector_iters=2`: The maximal number of correction steps in a single step.
* `max_lost_digits::Real`: The tracking is terminated if we estimate that we loose more than `max_lost_digits` in the linear algebra steps. This threshold depends on the `precision` argument.
* `max_refinement_iters=5`: The maximal number of correction steps used to refine the final value.
* `max_steps=1_000`: The maximal number of iterations the path tracker has available. Note that this changes to `10_000` for parameter homotopies.
* `max_step_size=Inf`: The maximal step size.
* `min_step_size=1e-14`: The minimal step size.
* `precision::PrecisionOption=PRECISION_FIXED_64`: The precision used for evaluating the residual in Newton's method.
* `predictor::AbstractPredictor`: The predictor used during in the predictor-corrector scheme. The default is [`Heun`](@ref)()`.
* `refinement_accuracy=1e-8`: The precision required for the final value.
* `simple_step_size_alg=false`: Use a more simple step size algorithm.
* `steps_jacobian_info_update::Int=1`: Every n-th step a linear system will be solved using a QR factorization to obtain an estimate for the condition number of the Jacobian.
* `terminate_ill_conditioned::Bool=true`: Indicates whether the path tracking should be terminated for ill-conditioned paths. A path is considerd ill-conditioned if the condition number of the Jacobian is larger than ≈1e14 or if it is larger than 1e`max_lost_digits`.

### Endgame specific options:
* `accuracy_eg::Float64=min(accuracy, 1e-5))`: It is possible to change the accuracy during the path tracking. Usually you want lower the accuracy.
* `cond_eg_start::Float64=1e4`: The endgame is only started if the condition of the Jacobian is larger than this threshold.
* `max_winding_number::Int=12`: This limits the maximal number of loops taken in applying Cauchy's formula.
* `min_cond_at_infinity::Float64=1e7`: A path is declared as going to infinity only if it's Jacobian is also larger than this threshold.
* `samples_per_loop::Int=12`: To compute singular solutions Cauchy's integral formula is used. The accuracy of the solutions increases with the number of samples per loop.
* `t_eg_start::Float64=0.1`: The endgame starts only if `t` is smaller than this threshold.
* `tol_val_inf_accurate::Float64=1e-4`: A valuation which would result in a path declared as going to infinity is only accepted if the estimated accuracy of the valuation is less than this threshold.
* `tol_val_finite_accurate::Float64=1e-3`: A valuation which would result in a proper solution is only accepted if the estimated accuracy of the valuation is less than this threshold. This is only affects solutions where the path has at some point near 0 a condition number larger than `cond_eg_start`.
It is recommended to also take a look at the [`PathTracker`](@ref) documentation for some context.

### Overdetermined system specific options:
* `overdetermined_min_accuracy=1e-5`: The minimal accuracy a non-singular solution needs to have to be considered a solution of the original system.
* `overdetermined_min_residual=1e-3`: The minimal residual a singular solution needs to have to be considered a solution of the original system.
"""
function solve end

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

function solve(args...; kwargs...)
    solve_kwargs, rest = splitkwargs(kwargs, solve_supported_keywords)
    tracker, start_solutions = pathtracker_startsolutions(args...; rest...)
    solve(tracker, start_solutions; solve_kwargs...)
end

function solve(tracker::Union{<:PathTracker,<:PolyhedralTracker}, start_solutions;
        path_result_details=:default, save_all_paths=false,
        path_jumping_check=true, kwargs...)
    all_results, ntracked = track_paths(tracker, start_solutions; path_result_details=path_result_details, save_all_paths=save_all_paths, kwargs...)
    if path_jumping_check
        path_jumping_check!(all_results, tracker, path_result_details; finite_results_only=!save_all_paths)
    end

    Result(all_results, ntracked, seed(tracker); multiplicity_tol=100*accuracy(tracker))
end

accuracy(T::PathTracker) = T.options.accuracy
accuracy(T::PolyhedralTracker) = accuracy(T.generic_tracker)

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

function track_paths(tracker, start_solutions;
                threading=true, show_progress=true,
                path_result_details::Symbol=:default, save_all_paths=false,
                # deprecated
                report_progress=nothing)
    results = Vector{result_type(tracker)}()
    n = length(start_solutions)

    @deprecatekwarg report_progress show_progress

    if show_progress
        progress = ProgressMeter.Progress(n; dt=0.1, desc="Tracking $n paths... ",
                                    clear_output_ijulia=true, delay=0.3)
    else
        progress = nothing
    end

    stats = SolveStats()
    ntracked = 0
    try
        nthreads = Threads.nthreads()
        if threading && nthreads > 1
            batch_size = 32 * nthreads
            batches = BatchIterator(start_solutions, batch_size)
            batch_tracker = BatchTracker(tracker, batches, path_result_details, save_all_paths)
            for batch in batches
                prepare_batch!(batch_tracker, batch, ntracked)
                ccall(:jl_threading_run, Ref{Cvoid}, (Any,), batch_tracker)

                for R in batch_tracker.results
                    if R !== nothing
                        push!(results, R)
                        is_success(R) && update!(stats, R)
                    end
                end
                ntracked += length(batch)
                if batch_tracker.interrupted
                    return results, ntracked
                end

                update_progress!(progress, ntracked, stats)
            end
        else
            for (k, s) in enumerate(start_solutions)
                return_code = track!(tracker, s)
                if save_all_paths ||
                   return_code == PathTrackerStatus.success ||
                   return_code == PathTrackerStatus.terminated_invalid_startvalue

                    R = PathResult(tracker, s, k; details=path_result_details)
                    push!(results, R)

                    if return_code == PathTrackerStatus.success
                        update!(stats, R)
                    end
                    ntracked = k
                end
                k % 32 == 0 && update_progress!(progress, k, stats)
            end
            # don't print if it already got printed above
            n % 32 != 0 && update_progress!(progress, n, stats)
        end
    catch e
        if isa(e, InterruptException)
            return results, ntracked
        else
            rethrow(e)
        end
    end
    results, n
end

function update!(stats::SolveStats, R::PathResult)
    if is_singular(R)
        stats.singular_real += is_real(R)
        stats.singular += 1
    else
        stats.regular_real += is_real(R)
        stats.regular += 1
    end
end

function update_progress!(progress, ntracked, stats::SolveStats; finished::Bool=false)
    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

mutable struct BatchTracker{Tracker<:PathTracker, V, R} <: Function
    trackers::Vector{Tracker}
    results::Vector{Union{Nothing, R}}
    ranges::Vector{UnitRange{Int}}
    batch::Vector{V}
    details::Symbol
    all_paths::Bool
    interrupted::Bool
    ntracked::Int
end

function BatchTracker(tracker::PathTracker, batches::BatchIterator, path_result_details::Symbol, save_all_paths::Bool)
    trackers = Threads.resize_nthreads!([tracker])
    batch_size = length(batches.batch)
    batch_results = Vector{Union{Nothing, result_type(tracker)}}(undef, batch_size)
    batch_results .= nothing
    ranges = partition_work(1:batch_size, length(trackers))
    BatchTracker(trackers, batch_results, ranges, batches.batch,
                 path_result_details, save_all_paths, false, 0)
end

function prepare_batch!(batch_tracker::BatchTracker, batch, ntracked)
    n = length(batch)
    if length(batch_tracker.results) ≠ n
        resize!(batch_tracker.results, n)
    end
    batch_tracker.results .= nothing
    batch_tracker.ntracked = ntracked
    partition_work!(batch_tracker.ranges, 1:n, length(batch_tracker.trackers))
end

function (batch::BatchTracker)()
    tid = Threads.threadid()
    try
        track_batch!(batch.results, batch.trackers[tid], batch.ranges[tid],
                 batch.batch, batch.details, batch.all_paths, batch.ntracked)
     catch e
         if isa(e, InterruptException)
             batch.interrupted = true
         else
             rethrow(e)
         end
     end
end
function track_batch!(results, pathtracker, range, starts, details, all_paths, ntracked)
    for k in range
        return_code = track!(pathtracker, starts[k])
        if all_paths ||
           return_code == PathTrackerStatus.success ||
           return_code == PathTrackerStatus.terminated_invalid_startvalue

            results[k] = PathResult(pathtracker, starts[k], ntracked+k; details=details)
        else
            results[k] = nothing
        end
    end
    nothing
end

"""
    path_jumping_check!(results, tracker, details)

Try to detect path jumping by comparing the winding numbers of finite results.
"""
function path_jumping_check!(results::Vector{<:PathResult}, tracker, details::Symbol; finite_results_only=false)
    if finite_results_only
        finite_results_indices = collect(1:length(results))
        finite_results = results
    else
        finite_results_indices = Int[]
        finite_results = Vector{eltype(results)}()
        for (i, r) in enumerate(results)
            if isfinite(r)
                push!(finite_results, r)
                push!(finite_results_indices, i)
            end
        end
    end
    !isempty(finite_results) || return results

    tol = tracker.core_tracker.options.refinement_accuracy
    # find cluster of multiple solutions
    clusters = multiplicities(solution, finite_results; tol=tol)
    while true
        rerun_paths = false
        for cluster in clusters
            m = length(cluster)
            # We have to consider a couple of cases
            # 1) All m solutions in the cluster also have winding number m
            #    -> Everything OK
            # 2) We find two non-singular solutions (winding_number === nothing)
            #    in the same cluster
            #    -> Path jumping happened
            # 3) There can be multiple paths going to the same solution, i.e.,
            #    winding number | m
            #    -> OK
            # 4) It can happen that paths going to a positive dimensional component
            #    have the same endpoint and still winding number = 1.
            #
            # This is quite complicated to tear apart. Therefore we focus for
            # now only on case 2)

            jumping = false
            for i in cluster
                cond = unpack(finite_results[i].condition_jacobian, 1.0)
                w = unpack(finite_results[i].winding_number, 0)
                if w ≤ 1 && cond < 1e10
                    jumping = true
                    break
                end
            end

            if jumping
                rerun_paths = true
                # rerun
                for i in cluster
                    rᵢ = finite_results[i]
                    new_rᵢ = track(tracker, start_solution(rᵢ);
                                            path_number=rᵢ.path_number,
                                            details=details,
                                            accuracy=min(1e-7, accuracy(tracker.core_tracker)),
                                            max_corrector_iters=1)
                    finite_results[i] = new_rᵢ
                    results[finite_results_indices[i]] = new_rᵢ
                end
            end
        end

        rerun_paths || break

        prev_clusters = clusters
        clusters = multiplicities(solution, finite_results; tol=tol)
        if clusters == prev_clusters
            break
        end
    end

    results
end
function path_jumping_check!(results::Vector{<:PathResult}, tracker::PolyhedralTracker, details::Symbol; kwargs...)
    path_jumping_check!(results, tracker.generic_tracker, details; kwargs...)
end

######################
## MultiplicityInfo ##
######################

"""
    MultiplicityInfo

This contains informations about the multiplicities of the solutions.
"""
struct MultiplicityInfo
    multiplicities::Dict{Int, Vector{Vector{Int}}}
    # This set stores all path numbers which we want to ignore if we don't show multiple results
    multiple_indicator::Set{Int32}
end

function MultiplicityInfo(pathresults::Vector{<:PathResult}; tol=1e-6)
    multiple_indicator = Set{Int32}()
    multiplicities = compute_multiplicities(pathresults; tol=float(tol))
    for clusters in values(multiplicities), cluster in clusters
        for i in 2:length(cluster)
            push!(multiple_indicator, pathresults[cluster[i]].path_number)
        end
    end
    MultiplicityInfo(multiplicities, multiple_indicator)
end

function compute_multiplicities(result::Vector{<:PathResult}; kwargs...) where T
    D = Dict{Int, Vector{Vector{Int}}}()
    compute_multiplicities!(D, result; kwargs...)
end
function compute_multiplicities!(D::Dict, result::Vector{<:PathResult}; tol::Float64 = 1e-6) where T
    for m in multiplicities(solution, result, tol = tol)
        if haskey(D, length(m))
            push!(D[length(m)], m)
        else
            D[length(m)] = [m]
        end
    end
    D
end

"""
    assign_multiplicities!(path_results::Vector{<:PathResult}, I::MultiplicityInfo)

Assign the path results their multiplicity according to the multiplicity info `I`.
"""
function assign_multiplicities!(path_results::Vector{<:PathResult}, I::MultiplicityInfo)
    #assign multiplicities
    for r in path_results
        r.multiplicity = 1
    end
    for (k, clusters) in I.multiplicities, cluster in clusters, vᵢ in cluster
        path_results[vᵢ].multiplicity = k
    end
    path_results
end

is_multiple_result(r::PathResult, I::MultiplicityInfo) = r.path_number ∈ I.multiple_indicator

"""
    Result{V<:AbstractVector}

The result of `solve`. This is a wrapper around the results of each single path ([`PathResult`](@ref)) and it contains some additional informations like
a random seed to replicate the result.
"""
mutable struct Result{V}
    pathresults::Vector{PathResult{V}}
    tracked_paths::Int
    seed::Int
    multiplicity_info::MultiplicityInfo
end

function Result(pathresults::Vector{PathResult{V}}, tracked_paths::Int, seed::Int;
                    multiplicity_tol=1e-6) where {V}
    multiplicity_info = MultiplicityInfo(pathresults; tol=multiplicity_tol)
    assign_multiplicities!(pathresults, multiplicity_info)
    Result(pathresults, tracked_paths, seed, multiplicity_info)
end

Base.length(r::Result) = length(r.pathresults)
Base.getindex(r::Result, I) = getindex(r.pathresults, I)

Base.iterate(r::Result) = iterate(r.pathresults)
Base.iterate(r::Result, state) = iterate(r.pathresults, state)
Base.lastindex(r::Result) = lastindex(r.pathresults)
Base.eltype(r::Type{Result{V}}) where {V} = PathResult{V}

is_multiple_result(r::PathResult, R::Result) = is_multiple_result(r, R.multiplicity_info)
is_multiple_result(r::PathResult, R::Vector{<:PathResult}) = false

"""
    multiplicities!(result::Result; tol=1e-6)

Compute the multiplicities of the solutions in `result` with respect to the given tolerance.
"""
function multiplicities!(result::Result; tol=1e-6)
    multiplicity_info = MultiplicityInfo(result.pathresults; tol=tol)
    result.multiplicity_info = multiplicity_info
    assign_multiplicities!(result.pathresults, multiplicity_info)
    result
end

const Results = Union{Result, Vector{<:PathResult}}
const ProjectiveResult = Result{<:PVector}

"""
    nresults(result; only_real=false, real_tol=1e-6, only_nonsingular=false, singular_tol=1e10, onlyfinite=true)

The number of solutions which satisfy the corresponding predicates.

## Example
```julia
result = solve(F)
# Get all non-singular results where all imaginary parts are smaller than 1e-8
nresults(result, only_real=true, real_tol=1e-8, only_nonsingular=true)
```
"""
function nresults(R::Results; only_real=false, real_tol=1e-6,
    only_nonsingular=false, only_singular=false, singular_tol=1e10, onlyfinite=true, multiple_results=false)
    count(R) do r
        (!only_real || is_real(r, real_tol)) &&
        (!only_nonsingular || is_nonsingular(r, singular_tol)) &&
        (!only_singular || is_singular(r, singular_tol)) &&
        (!onlyfinite || isfinite(r) || is_projective(r)) &&
        (multiple_results || !is_multiple_result(r, R))
    end
end

"""
    statistics(R::Result; only_real=false, real_tol=1e-6,
        only_nonsingular=false, only_singular=false, singular_tol=1e10)

Statistic about the number of (real) singular and non-singular solutions etc. Returns a named tuple with the statistics.

## Example
```julia
julia> statistics(solve([x^2+y^2-2, 2x+3y-1]))
(nonsingular = 2, singular = 0, real_nonsingular = 2, real_singular = 0, real = 2, at_infinity = 0, failed = 0, total = 2)
"""
function statistics(R::Results, only_real=false, real_tol=1e-6,
    only_nonsingular=false, only_singular=false, singular_tol=1e10)

    failed = at_infinity = nonsingular = singular = real_nonsingular = real_singular = 0
    singular_with_multiplicity = real_singular_with_multiplicity = 0
    for r in R
        is_multiple_result(r, R) && continue

        if is_failed(r)
            failed += 1
        elseif is_singular(r, singular_tol)
            if is_real(r, real_tol)
                real_singular += 1
                real_singular_with_multiplicity += unpack(multiplicity(r), 1)
            end
            singular += 1
            singular_with_multiplicity += unpack(multiplicity(r), 1)
        elseif !is_projective(r) && !isfinite(r)
            at_infinity += 1
        else # finite, nonsingular
            if is_real(r, real_tol)
                real_nonsingular += 1
            end
            nonsingular += 1
        end
    end
    (nonsingular = nonsingular,
    singular = singular,
    singular_with_multiplicity=singular_with_multiplicity,
    real_nonsingular = real_nonsingular,
    real_singular = real_singular,
    real_singular_with_multiplicity=real_singular_with_multiplicity,
    real = real_nonsingular + real_singular,
    at_infinity = at_infinity,
    failed = failed,
    total = R.tracked_paths)
end

"""
    nfinite(result)

The number of finite solutions.
"""
nfinite(R::Results) = count(isfinite, R)

"""
    nsolutions(result)

The number of solutions.
"""
nsolutions(R::Results) = nresults(R)

"""
    nsingular(result; singular_tol=1e10, multiplicitytol=1e-5, counting_multiplicities=false, kwargs...)

The number of singular solutions. A solution is considered singular if its windingnumber is larger than 1 or the condition number is larger than `tol`.
If `counting_multiplicities=true` the number of singular solutions times their multiplicities is returned.
"""
function nsingular(R::Results; singular_tol=1e10, counting_multiplicities=false, kwargs...)
    S = results(R; only_singular=true, multiple_results=false, singular_tol=singular_tol, kwargs...)
    isempty(S) && return 0
    counting_multiplicities && return sum(multiplicity, S)
    length(S)
end

"""
    nat_infinity(result)

The number of solutions at infinity.
"""
nat_infinity(R::Results) = count(is_at_infinity, R)

"""
    nafailed(result)

The number of failed paths.
"""
nfailed(R::Results) = count(is_failed, R)

"""
    nnonsingular(result; tol=1e-10)

The number of non-singular solutions.
"""
nnonsingular(R::Result; tol = 1e10) = count(r -> is_nonsingular(r, tol), R)

"""
    nreal(result; tol=1e-6)

The number of real solutions where all imaginary parts of each solution
are smaller than `tol`.
"""
nreal(R::Results; tol = 1e-6) = count(r -> is_real(r, tol), R)

"""
    ntracked(R::Result)

Returns the total number of paths tracked.
"""
ntracked(R::Result) = R.tracked_paths

"""
    seed(result)

The random seed used in the computation.
"""
seed(result::Result) = result.seed

# Filtering
"""
    results(result; only_real=false, real_tol=1e-6, only_nonsingular=false,
                onlysigular=false, singular_tol=1e10, onlyfinite=true, multiple_results=false)

Return all `PathResult`s for which the given conditions apply.

## Example

```julia
R = solve(F)

# This gives us all PathResults considered non-singular and real (but still as a complex vector).
real_solutions = results(R, only_real=true, only_nonsingular=true)
```
"""
results(R::Results; kwargs...) = mapresults(identity, R; kwargs...)
results(f::Function, R::Results; kwargs...) = mapresults(f, R; kwargs...)

"""
    mapresults(f::Function, result; conditions...)

Apply the function `f` to all `PathResult`s for which the given conditions apply. For the possible
conditions see [`results`](@ref).

## Example
```julia
# This gives us all solutions considered real (but still as a complex vector).
real_solutions = mapresults(solution, R, only_real=true)
```
"""
function mapresults(f::Function, R::Results;
    only_real=false, real_tol=1e-6, only_nonsingular=false, only_singular=false, singular_tol=1e10,
    onlyfinite=true, multiple_results=false)
    [f(r) for r in R if
        (!only_real || is_real(r, real_tol)) &&
        (!only_nonsingular || is_nonsingular(r, singular_tol)) &&
        (!only_singular || is_singular(r, singular_tol)) &&
        (!onlyfinite || isfinite(r) || is_projective(r)) &&
        (multiple_results || !is_multiple_result(r,R))]
end

"""
    solutions(result; conditions...)

Return all solution (as `Vector`s) for which the given conditions apply.
For the possible `conditions` see [`results`](@ref).

## Example
```julia
julia> @polyvar x y
julia> result = solve([(x-2)y, y+x+3]);
julia> solutions(result)
[[2.0+0.0im, -5.0+0.0im], [-3.0+0.0im, 0.0+0.0im]]
```
"""
function solutions(result::Results; kwargs...)
    mapresults(solution, result; kwargs...)
end

"""
    real_solutions(result; tol=1e-6, conditions...)

Return all real solution (as `Vector`s of reals) for which the given conditions apply.
For the possible `conditions` see [`results`](@ref). Note that `only_real` is always `true`
and `real_tol` is now `tol`.

## Example
```julia
julia> @polyvar x y
julia> result = solve([(x-2)y, y+x+3]);
julia> real_solutions(result)
[[2.0, -5.0], [-3.0, 0.0]]
```
"""
function real_solutions(result::Results; only_real=true, tol=1e-6, kwargs...)
    mapresults(r -> real.(solution(r)), result; only_real=true, real_tol=tol, kwargs...)
end
@deprecate realsolutions(result; kwargs...) real_solutions(result; kwargs...)

"""
    nonsingular(result::Results; conditions...)

Return all `PathResult`s for which the solution is non-singular. This is just a shorthand
for `results(R; only_nonsingular=true, conditions...)`. For the possible `conditions` see [`results`](@ref).
"""
nonsingular(R::Results; kwargs...) = results(R; only_nonsingular=true, kwargs...)

"""
    singular(R::Results; tol=1e10, multiple_results=false, kwargs...)

Return all `PathResult`s for which the solution is singular. A solution is labeled singular if the condition number is greater than `singular_tol`, or if the winding number is > 1.
If `multiple_results=false` only one point from each cluster of multiple solutions is returned.
If If `multiple_results=true` all singular solutions in `R` are returned.
For the possible `kwargs` see [`results`](@ref).
"""
function singular(R::Results; tol=1e10, kwargs...)
    results(R; only_singular = true, singular_tol=tol, kwargs...)
end

"""
    finite(result::AffineResults; conditions...)

Return all `PathResult`s for which the solution is finite. This is just a shorthand
for `results(R; onlyfinite=true, conditions...)`. For the possible `conditions` see [`results`](@ref).
"""
finite(R::Results; kwargs...) = results(R; onlyfinite=true, kwargs...)

"""
    real(result, tol=1e-6)

Get all results where the solutions are real with the given tolerance `tol`.
See [`is_real`](@ref) for details regarding the determination of 'realness'.
"""
Base.real(R::Results; tol=1e-6) = [r for r in R if is_real(r, tol)]

"""
    failed(result)

Get all results where the path tracking failed.
"""
failed(R::Results) = [r for r in R if is_failed(r)]

"""
    at_infinity(result::AffineResult)

Get all results where the solutions is at infinity.
"""
at_infinity(R::Results) = [r for r in R if is_at_infinity(r)]

"""
    multiplicities(V::Results; tol=1e-6)

Returns a `Vector` of `Vector{PathResult}`s grouping the `PathResult`s whose solutions appear with multiplicities *greater* 1 in 'V'.
Two solutions are regarded as equal, when their pairwise distance is less than 'tol'.
"""
function multiplicities(results::Results; tol=1e-6)
    map(i -> results[i], multiplicities(solution, results; tol = tol))
end

####Show function

function Base.show(io::IO, x::Result)
    s = statistics(x)
    println(io, "Result with $(s.nonsingular + s.singular) solutions")
    println(io, "==================================")
    println(io, "• $(s.nonsingular) non-singular $(plural("solution", s.nonsingular)) ($(s.real_nonsingular) real)")
    println(io, "• $(s.singular) singular $(plural("solution", s.singular)) ($(s.real_singular) real)")
    s.at_infinity > 0 &&
        println(io, "• $(s.at_infinity) $(plural("solution", s.at_infinity)) at infinity")
    s.failed > 0 &&
        println(io, "• $(s.failed) failed $(plural("path", s.failed))")
    println(io, "• $(ntracked(x)) paths tracked")
    println(io, "• random seed: $(seed(x))")
    if s.singular > 0
        println(io, "• multiplicity table of singular solutions:")
        singular_multiplicities_table(io, x, s)
    end
end

function Base.show(io::IO, x::ProjectiveResult)
    s = statistics(x)
    println(io, "Result with $(s.nonsingular + s.singular) solutions")
    println(io, "==================================")
    println(io, "• $(s.nonsingular) non-singular $(plural("solution", s.nonsingular)) ($(s.real_nonsingular) real)")
    println(io, "• $(s.singular) singular $(plural("solution", s.singular)) ($(s.real_singular) real)")
    s.failed > 0 &&
        println(io, "• $(s.failed) failed $(plural("path", s.failed))")
    println(io, "• $(ntracked(x)) paths tracked")
    println(io, "• random seed: $(seed(x))")
    if s.singular > 0
        println(io, "• multiplicity table of singular solutions:")
        singular_multiplicities_table(io, x, s)
    end
end

TreeViews.hastreeview(::Result) = true
TreeViews.hastreeview(::ProjectiveResult) = true
TreeViews.numberofnodes(::Result) = 8
TreeViews.numberofnodes(::ProjectiveResult) = 7
TreeViews.treelabel(io::IO, x::Result, ::MIME"application/prs.juno.inline") =
    print(io, "<span class=\"syntax--support syntax--type syntax--julia\">Result</span>")

function TreeViews.nodelabel(io::IO, x::Result, i::Int, ::MIME"application/prs.juno.inline")
    s = statistics(x)
    if i == 1
        print(io, "Paths tracked")
    elseif i == 2 && s.nonsingular > 0
        print(io, "$(s.nonsingular) non-singular ($(s.real_nonsingular) real)")
    elseif i == 3 && s.singular > 0
        print(io, "$(s.singular) singular ($(s.real_singular) real)")
    elseif i == 4 && (s.real_nonsingular+s.real_singular) > 0
        print(io, "$(s.real_nonsingular+s.real_singular) real")
    elseif i == 5 && s.at_infinity > 0
        print(io, "$(s.at_infinity) at_infinity")
    elseif i == 6 && s.failed > 0
        print(io, "$(s.failed) failed")
    elseif i == 7
        print(io, "Random seed used")
    elseif i == 8 && s.singular > 0
        print(io, "  multiplicity table of singular solutions: \n")
        singular_multiplicities_table(io, x, s)
    end
end

function TreeViews.treenode(r::Result, i::Integer)
    s = statistics(r)
    if i == 1
        return ntracked(r)
    elseif i == 2 && s.nonsingular > 0
        return finite(r, only_nonsingular=true)
    elseif i == 3 && s.singular > 0
        return finite(r, only_singular=true)
    elseif i == 4 && (s.real_nonsingular+s.real_singular) > 0
        return finite(r, only_real = true)
    elseif i == 5 && s.at_infinity > 0
        return at_infinity(r)
    elseif i == 6 && s.failed > 0
        return failed(r)
    elseif i == 7
        return seed(r)
    elseif i == 8 && s.singular > 0
        return missing
    end
    missing
end

function TreeViews.nodelabel(io::IO, x::ProjectiveResult, i::Int, ::MIME"application/prs.juno.inline")
    s = statistics(x)
    if i == 1
        print(io, "Paths tracked")
    elseif i == 2 && s.nonsingular > 0
        print(io, "$(s.nonsingular) non-singular ($(s.real_nonsingular) real)")
    elseif i == 3 && s.singular > 0
        print(io, "$(s.singular) singular ($(s.real_singular) real)")
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        print(io, "$(s.real_nonsingular + s.real_singular) real solutions")
    elseif i == 5 && s.failed > 0
        print(io, "$(s.failed) failed")
    elseif i == 6
        print(io, "Random seed used")
    elseif i == 7 && s.singular > 0
        print(io, "  multiplicity table of singular solutions: \n")
        singular_multiplicities_table(io, x, s)
    end
end

function TreeViews.treenode(r::ProjectiveResult, i::Integer)
    s = statistics(r)
    if i == 1
        return length(r)
    elseif i == 2 && s.nonsingular > 0
        return finite(r, only_nonsingular=true)
    elseif i == 3 && s.singular > 0
        return finite(r, only_singular=true)
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        return finite(r, only_real=true)
    elseif i == 5 && s.failed > 0
        return failed(r)
    elseif i == 6
        return seed(r)
    elseif i == 7 && s.singular > 0
        return missing
    end
    missing
end

plural(singularstr, n) = n == 1 ? singularstr : singularstr * "s"

function singular_multiplicities_table(io, result::Result, stats = statistics(result))
    M = result.multiplicity_info.multiplicities
    if isempty(M)
        n_higher_mults_total = 0
    else
        n_higher_mults_total = sum(length, values(M))
    end

    headers = ["mult.", "total", "# real", "# non-real"]
    mult_1_exists = n_higher_mults_total < stats.singular
    data = Matrix{Int}(undef, mult_1_exists+length(M), 4)

    n_higher_real = 0
    i = mult_1_exists + 1
    for k in sort(collect(keys(M)))
    	data[i, 1] = k
    	n_real_solsᵢ = count(Mᵢ -> is_real(result.pathresults[Mᵢ[1]]), M[k])
    	n_solsᵢ = length(M[k])
        data[i, 2] = n_solsᵢ
    	data[i, 3] = n_real_solsᵢ
    	data[i, 4] = n_solsᵢ - n_real_solsᵢ
    	n_higher_real += n_real_solsᵢ
        i += 1
    end

    if mult_1_exists
        data[1,1] = 1
        data[1,2] = stats.singular - n_higher_mults_total
    	data[1,3] = stats.real_singular - n_higher_real
    	data[1,4] = data[1,2] - data[1,3]
    end

    PrettyTables.pretty_table(io, data, headers,
            PrettyTables.unicode;
            alignment=:c,
            header_crayon=PrettyTables.Crayon(bold=false),
            border_crayon=PrettyTables.Crayon(faint=true))
end


text:
function switch_pkgbutler_channel(package_name::AbstractString, channel::Symbol)
    ctx = Pkg.Types.Context()
    haskey(ctx.env.project.deps, package_name) || error("Unkonwn package $package_name.")
    pkg_uuid = ctx.env.project.deps[package_name]
    pkg_path = ctx.env.manifest[pkg_uuid].path
    pkg_path===nothing && error("Package must be deved to enable the Julia Package Butler.")

    PkgButlerEngine.configure_pkg(pkg_path; channel=channel)

    return
end

function switch_pkgbutler_template(package_name::AbstractString, template::Symbol)
    ctx = Pkg.Types.Context()
    haskey(ctx.env.project.deps, package_name) || error("Unkonwn package $package_name.")
    pkg_uuid = ctx.env.project.deps[package_name]
    pkg_path = ctx.env.manifest[pkg_uuid].path
    pkg_path===nothing && error("Package must be deved to enable the Julia Package Butler.")

    PkgButlerEngine.configure_pkg(pkg_path; template=template)

    return
end

function enable_pkgbutler(package_name::AbstractString; channel=:auto, template=:auto, force::Bool=false)
    ctx = Pkg.Types.Context()
    haskey(ctx.env.project.deps, package_name) || error("Unkonwn package $package_name.")
    pkg_uuid = ctx.env.project.deps[package_name]
    pkg_path = ctx.env.manifest[pkg_uuid].path
    pkg_path===nothing && error("Package must be deved to enable the Julia Package Butler.")

    if !isabspath(pkg_path)
        pkg_path = abspath(joinpath(dirname(ctx.env.manifest_file), pkg_path))
    end

    pkg_owner_repo_name = nothing

    pkg_repo = GitRepo(pkg_path)
    try
        pkg_remote = LibGit2.lookup_remote(pkg_repo, "origin")

        pkg_remote===nothing && error("The package must have a remote called origin.")

        pkg_url = LibGit2.url(pkg_remote)

        pkg_owner_repo_name = get_repo_onwer_from_url(pkg_url)
    finally
        close(pkg_repo)
    end

    configure_deploykey(pkg_owner_repo_name, force)

    PkgButlerEngine.configure_pkg(pkg_path; channel=channel, template=template)

    return
end

import GitHub.DEFAULT_API

mutable struct DeployKey <: GitHub.GitHubType
    id::Union{Int, Nothing}
    key::Union{String, Nothing}
    url::Union{GitHub.HTTP.URI, Nothing}
    title::Union{String, Nothing}
    verified::Union{Bool, Nothing}
    created_at::Union{GitHub.Dates.DateTime, Nothing}
    read_only::Union{Bool, Nothing}
end

DeployKey(data::Dict) = GitHub.json2github(DeployKey, data)

GitHub.namefield(key::DeployKey) = key.id

GitHub.@api_default function deploykey(api::GitHub.GitHubAPI, repo, deploykey_obj; options...)
    result = GitHub.gh_get_json(api, "/repos/$(GitHub.name(repo))/keys/$(GitHub.name(deploykey_obj))"; options...)
    return DeployKey(result)
end

GitHub.@api_default function deploykeys(api::GitHub.GitHubAPI, repo; options...)
    results, page_data = GitHub.gh_get_paged_json(api, "/repos/$(GitHub.name(repo))/keys"; options...)
    return map(DeployKey, results), page_data
end

GitHub.@api_default function create_deploykey(api::GitHub.GitHubAPI, repo; options...)
    result = GitHub.gh_post_json(api, "/repos/$(GitHub.name(repo))/keys"; options...)
    return DeployKey(result)
end

GitHub.@api_default function delete_deploykey(api::GitHub.GitHubAPI, repo, item; options...)
    return GitHub.gh_delete(api, "/repos/$(GitHub.name(repo))/keys/$(GitHub.name(item))"; options...)
end

function configure_deploykey(repo, force)
    pub_deploykey, private_deploykey = mktempdir() do f
        p = joinpath(f, "deploykey")
        run(`ssh-keygen -q -t rsa -f $p -N ''`)

        pub_deploykey = read("$p.pub", String)
        private_deploykey = read(p, String)

        rm(p, force=true)
        rm("$p.pub", force=true)

        return pub_deploykey, private_deploykey
    end

    creds = LibGit2.GitCredential(GitConfig(), "https://github.com")

    # TODO Check for creds===nothing if there are no credentials stored
    myauth = GitHub.authenticate(read(creds.password, String))
    Base.shred!(creds.password)

    existing_deploykeys = PkgDev.deploykeys(repo, auth=myauth)

    existing_pkgbutler_deploykeys_indices = findall(i->i.title=="Julia Package Butler", existing_deploykeys[1])

    if !force && length(existing_pkgbutler_deploykeys_indices) > 0
        println("Found one or more existing deployment keys on Github for $repo that are named 'Julia Package Butler'. Should these be deleted? (yes/no)")
        answer = readline()
        if !(answer=="" || lowercase(answer)=="yes" || lowercase(answer)=="y")
            empty!(existing_pkgbutler_deploykeys_indices)
        end
    end

    for dk in existing_pkgbutler_deploykeys_indices
        delete_deploykey(repo, existing_deploykeys[1][dk], auth=myauth)
    end

    create_deploykey(repo, params=Dict("title"=>"Julia Package Butler", "key"=>pub_deploykey, "read_only"=>false), auth=myauth)

    println()
    println("Add the following text as a GitHub Actions Secret with the name JLPKGBUTLER_TOKEN:")
    println()
    println(Base64.base64encode(private_deploykey))
    println()
    println("Add the repository upload token you find at the following URL as a GitHub Actions Secret with the name CODECOV_TOKEN:")
    println()
    println("https://codecov.io/gh/$repo/settings")
end

ext: jl


text:
using CuArrays: @cuindex, cudims

function mapreducedim_kernel_serial(f, op, R, A, range)
    I = @cuindex R
    newrange = map((r, i) -> r === nothing ? i : r, range, I)
    for I′ in CartesianIndices(newrange)
        @inbounds R[I...] = op(R[I...], f(A[I′]))
    end
    return
end

@inline function reduce_block(arr::CuDeviceArray, op)
    sync_threads()
    len = blockDim().x
    while len != 1
        sync_threads()
        skip = (len + 1) >> 1
        reduce_to = threadIdx().x - skip
        if 0 < reduce_to <= (len >> 1)
            arr[reduce_to] = op(arr[reduce_to], arr[threadIdx().x])
        end
        len = skip
    end
    sync_threads()
end

function mapreducedim_kernel_parallel(f, op, R::CuDeviceArray{T}, A::CuDeviceArray{T},
                             CIS, Rlength, Slength) where {T}
    for Ri_base in 0:(gridDim().x * blockDim().y):(Rlength-1)
        Ri = Ri_base + (blockIdx().x - 1) * blockDim().y + threadIdx().y
        Ri > Rlength && return
        RI = Tuple(CartesianIndices(R)[Ri])
        S = @cuStaticSharedMem(T, 512)
        Si_folded_base = (threadIdx().y - 1) * blockDim().x
        Si_folded = Si_folded_base + threadIdx().x
        # serial reduction of A into S by Slength ÷ xthreads
        for Si_base in 0:blockDim().x:(Slength-1)
            Si = Si_base + threadIdx().x
            Si > Slength && break
            SI = Tuple(CIS[Si])
            AI = ifelse.(size(R) .== 1, SI, RI)
            if Si_base == 0
                S[Si_folded] = f(A[AI...])
            else
                S[Si_folded] = op(S[Si_folded], f(A[AI...]))
            end
        end
        # block-parallel reduction of S to S[1] by xthreads
        reduce_block(view(S, (Si_folded_base + 1):512), op)
        # reduce S[1] into R
        threadIdx().x == 1 && (R[Ri] = op(R[Ri], S[Si_folded]))
    end
    return
end

function Base._mapreducedim!(f, op, R::CuArray{T}, A::CuArray{T}) where {T}
    # the kernel as generated from `f` and `op` can require lots of registers (eg. #160),
    # so we need to be careful about how many threads we launch not to run out of them.
    Rlength = length(R)
    Ssize = ifelse.(size(R) .== 1, size(A), 1)
    Slength = prod(Ssize)
    CIS = CartesianIndices(Ssize)

    parallel_args = (f, op, R, A, CIS, Rlength, Slength)
    GC.@preserve parallel_args begin
        parallel_kargs = cudaconvert.(parallel_args)
        parallel_tt = Tuple{Core.Typeof.(parallel_kargs)...}
        parallel_kernel = cufunction(mapreducedim_kernel_parallel, parallel_tt)

        # we are limited in how many threads we can launch...
        ## by the kernel
        kernel_threads = CUDAnative.maxthreads(parallel_kernel)
        ## by the device
        dev = CUDAdrv.device()
        block_threads = (x=attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_MAX_BLOCK_DIM_X),
                         y=attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_MAX_BLOCK_DIM_Y),
                         total=attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK))

        # figure out a legal launch configuration
        y_thr = min(nextpow(2, Rlength ÷ 512 + 1), 512, block_threads.y, kernel_threads)
        x_thr = min(512 ÷ y_thr, Slength, block_threads.x,
                    ceil(Int, block_threads.total/y_thr),
                    ceil(Int, kernel_threads/y_thr))

        if x_thr >= 8
            blk, thr = (Rlength - 1) ÷ y_thr + 1, (x_thr, y_thr, 1)
            parallel_kernel(parallel_kargs...; threads=thr, blocks=blk)
        else
            # not enough work, fall back to serial reduction
            range = ifelse.(length.(axes(R)) .== 1, axes(A), nothing)
            blk, thr = cudims(R)
            @cuda(blocks=blk, threads=thr, mapreducedim_kernel_serial(f, op, R, A, range))
        end
    end

    return R
end

import Base.minimum, Base.maximum, Base.reduce

_reduced_dims(x::CuArray, ::Colon) = Tuple(Base.ones(Int, ndims(x)))
_reduced_dims(x::CuArray, dims) = Base.reduced_indices(x, dims)

_initarray(x::CuArray{T}, dims, init) where {T} = fill!(similar(x, T,     _reduced_dims(x, dims)), init)

function _reduce(op, x::CuArray, init, ::Colon)
    mx = _initarray(x, :, init)
    Base._mapreducedim!(identity, op, mx, x)
    collect(mx)[1]
end

function _reduce(op, x::CuArray, init, dims)
    mx = _initarray(x, dims, init)
    Base._mapreducedim!(identity, op, mx, x)
end

"""
    reduce(op, x::CuArray; dims=:, init)

The initial value `init` is mandatory for `reduce` on `CuArray`'s. It must be a neutral element for `op`.
"""
reduce(op, x::CuArray; dims=:, init) = _reduce(op, x, init, dims)

minimum(x::CuArray{T}; dims=:) where {T} = _reduce(min, x, typemax(T), dims)
maximum(x::CuArray{T}; dims=:) where {T} = _reduce(max, x, typemin(T), dims)


text:
function write_fixed_width(io::IO, minSet::Vector{Tuple{T, T, T}}; sorted=false, header_add=repeat(" ", 10)) where T
    length(header_add) != 10 && throw(ArgumentError("header string must be 10 characters long"))
    ll = sorted ? sort(minSet) : collect(minSet)
    @printf(io, "# === Header Start ===        \n")
    @printf(io, "    # of Elements             \n")
    @printf(io, "  %12d      %s\n", length(minSet), header_add)
    @printf(io, "# === Header End ===          \n")
    for k in ll
        @printf(io, " %9d %9d %9d\n", k...)
    end
end

function write_fixed_width(io::IO, minSet::Vector; sorted=false, header_add=repeat(" ", 10))
    length(header_add) != 10 && throw(ArgumentError("header string must be 10 characters long"))
    ll = sorted ? sort(minSet) : collect(minSet)
    write(io, "# === Header Start ===        \n")
    write(io, "    # of Elements             \n")
    write(io, rpad(" "*string(length(minSet))*" "*header_add, 30, " ")*"\n")
    write(io, "# === Header End ===          \n")
    for k in ll
        for ki in k
            write(io, rpad(k,10))
        end
        write(io, "\n")
    end
end

io_functions = (:write_fixed_width, )

# For all write functions f(io::IO, args...) generate wrappers f(args...) which writes
# directly to stdout and f(s::String, args...) which opens/closes a file and writes to it.
for ff in io_functions
    @eval $ff(s::String, arg, args...; kwargs...) = (open(s,"w") do io; $ff(io,arg,args...; kwargs...); end)
    @eval $ff(arg, args...; kwargs...) = $ff(stdout, arg, args...; kwargs...)
end


text:
@connector function Pin(;name)
    sts = @variables begin
        v(t)                    # Potential at the pin [V]
        i(t), [connect=Flow]    # Current flowing into the pin [A]
    end
    ODESystem(Equation[], t, sts, [], name=name, defaults=Dict(v=>1.0, i=>1.0))
end
Base.@doc "Port for an electrical system." Pin

"""
    OnePort(;name, v_start=0.0, i_start=0.0)

Component with two electrical pins p and n and current i from p to n.

# Parameters:
- `v_start`: [V] Initial voltage across the component
- `i_start`: [A] Initial current through the component
"""
function OnePort(;name, v_start=0.0, i_start=0.0)
    @named p = Pin()
    @named n = Pin()
    sts = @variables begin
        v(t)=v_start
        i(t)=i_start
    end
    eqs = [
        v ~ p.v - n.v
        0 ~ p.i + n.i
        i ~ p.i
    ]
    return compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

@connector function DigitalPin(; name)
    @variables val(t) v(t) i(t)
    eqs = [
        val ~ IfElse.ifelse((0.0 <= v) & (v <= 0.8) | (2.0 <= v) & (v <= 5.0),
                                IfElse.ifelse(v > 2.0, 1, 0), X)
    ]
    ODESystem(Equation[], t, [val, v, i], [], defaults=Dict(val=>0, i=>0), name=name)
end


text:
function e(x)
    n = 1
    res = 1.0
    for i in 1:30
        n *= i
        res += 1 / n;
    end

    res
end

println("e=", e(1))


text:
using Documenter, PersistentCohomology

makedocs(
    modules = [PersistentCohomology],
    format = :html,
    sitename = "PersistentCohomology.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/piever/PersistentCohomology.jl.git",
    target = "build",
    julia = "1.0",
    deps = nothing,
    make = nothing,
)


text:
@testset "Server/Paths                " begin
    # requested path --> a filesystem path
    bk = pwd()
    cd(dirname(dirname(pathof(LiveServer))))
    req = "tmp"
    @test LS.get_fs_path(req) == ""
    req = "/test/dummies/index.html"
    @test LS.get_fs_path(req) == "test/dummies/index.html"
    req = "/test/dummies/r%C3%A9sum%C3%A9/"
    @test LS.get_fs_path(req) == "test/dummies/résumé/index.html"
    cd(bk)
end

#=
NOTE: if extending these tests, please be careful. As they involve @async tasks which,
themselves, spawn async tasks, if your tests fail for some reason you will HAVE to kill
the current Julia session and restart one otherwise the tasks that haven't been killed
due to the tests not being finished properly will keep running and may clash with new
tasks that you will try to start.
=#
@testset "Server/Step-by-step testing " begin
    #
    # STEP 0: cd to dummies
    #
    bk = pwd()
    cd(dirname(dirname(pathof(LiveServer))))
    cd(joinpath("test", "dummies"))
    port = 8123

    isdir("css") && rm("css", recursive=true)
    isfile("tmp.html") && rm("tmp.html")

    write("tmp.html", "blah")
    mkdir("css")
    write("css/foo.css", "body { color: pink; }")

    #
    # STEP 1: launching the listener
    #
    # assert 8000 ≤ port ≦ 9000
    @test_throws ArgumentError serve(port=7000)
    @test_throws ArgumentError serve(port=10000)

    # define filewatcher outside so that can follow it
    fw = LS.SimpleWatcher(LS.file_changed_callback)
    task = @async serve(fw; port=port)
    sleep(0.1) # give it time to get started

    # there should be a callback associated with fw now
    @test fw.status == :runnable
    # the filewatcher should be running
    @test LS.is_running(fw)
    # it also should be empty thus far
    @test isempty(fw.watchedfiles)

    #
    # STEP 2: triggering a request
    #
    response = HTTP.get("http://localhost:$port/")
    @test response.status == 200
    # the browser script should be appended
    @test String(response.body) == replace(read("index.html", String),
                            ""=>"$(LS.BROWSER_RELOAD_SCRIPT)")
    # if one asks for something incorrect, a 404 should be returned
    # XXX ok so actually an ERROR is thrown, that's not good?
    @test_throws HTTP.ExceptionRequest.StatusError HTTP.get("http://localhost:$port/no.html")
    # if one asks for something without a , it should just be appended
    response = HTTP.get("http://localhost:$port/tmp.html")
    @test response.status == 200
    @test String(response.body) == read("tmp.html", String) * LS.BROWSER_RELOAD_SCRIPT
    response = HTTP.get("http://localhost:$port/css/foo.css")
    @test String(response.body) == read("css/foo.css", String)
    @test response.status == 200

    # we asked earlier for index.html therefore that file should be followed
    @test fw.watchedfiles[1].path == "index.html"
    # also tmp
    @test fw.watchedfiles[2].path == "tmp.html"
    # and the css
    @test fw.watchedfiles[3].path == "css/foo.css"

    # if we modify the file, it should trigger the callback function which should open
    # and then subsequently close a websocket. We check this happens properly by adding
    # our own sentinel websocket
    sentinel = HTTP.WebSockets.WebSocket(IOBuffer())
    LS.WS_VIEWERS["tmp.html"] = [sentinel]

    @test sentinel.io.writable
    write("tmp.html", "something new")
    sleep(0.1)
    # the sentinel websocket should be closed
    @test !sentinel.io.writable
    # the websockets should have been flushed
    @test isempty(LS.WS_VIEWERS["tmp.html"])

    # let's do this again with an infra file which will ping all websockets
    sentinel1 = HTTP.WebSockets.WebSocket(IOBuffer())
    sentinel2 = HTTP.WebSockets.WebSocket(IOBuffer())
    push!(LS.WS_VIEWERS["tmp.html"], sentinel1)
    LS.WS_VIEWERS["css/foo.css"] = [sentinel2]
    write("css/foo.css", "body { color:blue; }")
    sleep(0.1)
    # all sentinel websockets should be closed
    @test !sentinel1.io.writable
    @test !sentinel2.io.writable

    # if we remove the files, it shall stop following it
    rm("tmp.html")
    rm("css", recursive=true)
    sleep(0.25)
    # only index.html is still watched
    @test length(fw.watchedfiles) == 1

    #
    # SHUTTING DOWN
    #
    # this should have interrupted the server, so it should be possible
    # to restart one on the same port (otherwise this will throw an error, already in use)
    schedule(task, InterruptException(), error=true)
    sleep(0.25) # give it time to get done
    @test istaskdone(task)
    @test begin
        server = Sockets.listen(port)
        sleep(0.1)
        close(server)
        true
    end == true

    # Check that WS_FILES is properly destroyed
    isempty(LS.WS_VIEWERS)

    cd(bk)
end

@testset "Server/ws_upgrade testing   " begin
    io = IOBuffer()
    s = HTTP.Stream(HTTP.Request("GET", "http://httpbin.org/ip"), io)

    # ws_upgrade
    ws = LS.ws_upgrade(s)

    @test ws.server

    @test occursin("HTTP/1.1 101 Switching Protocols\r\nUpgrade: websocket\r\nConnection: Upgrade\r\nSec-WebSocket-Accept:", String(take!(ws.io)))

    test_string = "blah bLah"
    write(io, test_string)
    @test String(take!(ws.io)) == test_string
    @test isempty(ws.io.data)
end

@testset "Server/ws_tracker testing   " begin
    bk = pwd()
    cd(mktempdir())
    write("test_file.html", "Hello!")

    io = IOBuffer()
    s = HTTP.Stream(HTTP.Request("GET", "http://localhost:8562/test_file.html"), io)

    fs_path = LS.get_fs_path(s.message.target)
    @test fs_path == "test_file.html"

    ws = LS.ws_upgrade(s)
    write(io, "some stuff on the websocket")

    tsk = @async LS.ws_tracker(ws, s.message.target)
    sleep(0.1)

    # simulate a "good" closure (an event caused a write on the websocket and then closes it)
    close(ws.io)
    sleep(0.2)
    @test istaskdone(tsk)
    @test !LS.WS_INTERRUPT[]

    # the websocket should have been added to the list
    @test LS.WS_VIEWERS[fs_path] isa Vector{HTTP.WebSockets.WebSocket}
    @test length(LS.WS_VIEWERS[fs_path]) == 1
    @test LS.WS_VIEWERS[fs_path][1] == ws

    io = IOBuffer()
    s = HTTP.Stream(HTTP.Request("GET", "http://localhost:8562/test_file.html"), io)
    ws = LS.ws_upgrade(s)
    write(io, "la di da, dimension C-137")
    tsk = @async LS.ws_tracker(ws, s.message.target)
    sleep(0.2)

    # simulate a "bad" closure
    schedule(tsk, InterruptException(), error=true)
    sleep(0.1)
    @test istaskdone(tsk)
    @test LS.WS_INTERRUPT[]

    @test length(LS.WS_VIEWERS[fs_path]) == 2
    @test LS.WS_VIEWERS[fs_path][2] == ws

    # cleanup
    empty!(LS.WS_VIEWERS)

    cd(bk)
end


text:
module QuadGG

using Printf

export adaptlob, adaptsim

mutable struct WarnFlag{b<:Bool}
    nowarn::b
end

include("adaptlob.jl")
include("adaptsim.jl")

end # module


text:
function reallocate_ipfs(
    existing_ipfs::Vector{T}, proposed_ipfs::Vector{T}
) where {T<:AbstractString}
    return existing_ipfs ∩ proposed_ipfs
end
function open_ipfs(
    proposed_ipfs::Vector{T}, reallocate_ipfs::Vector{T}
) where {T<:AbstractString}
    return setdiff(proposed_ipfs, reallocate_ipfs)
end
function close_ipfs(
    existing_ipfs::Vector{T}, reallocate_ipfs::Vector{T}, frozenlist::Vector{T}
) where {T<:AbstractString}
    return setdiff(setdiff(existing_ipfs, reallocate_ipfs), frozenlist)
end


text:
using PyPlot
"Rounding for legends"
round_sig(x; sig = 3) = round(x, digits = sig - Int(floor(log10(abs(x))))-1)

"Pyplot 3d surface and countour plot"
function wrap_surface_and_countour(np::Npos; title = "3d surface", levels = 100)
    wrap_surface_and_countour(slice_1(np), slice_2(np), slice_3(np); title = title, levels = levels)
end
function wrap_surface_and_countour(mat_x, mat_y, mat_z; title = "3d surface", levels = 100)
    clf()
    using3D()
    fig = figure("pyplot_surfaceplot", figsize=(10, 10)) # inches
    ax = fig.add_subplot(2,1,1, projection="3d")
    nx, ny = size(mat_x)
    plot_surface(mat_x, mat_y, mat_z, rcount=nx, ccount = ny, edgecolors="k",
                 cmap=ColorMap("Blues"), alpha=0.8, linewidth=0.25)
    xlabel("x")
    ylabel("y")
    PyPlot.title(title)
    #subplot(212)
    ax2 = fig.add_subplot(2,1,2, aspect = "equal")
    ax2.contourf(mat_x, mat_y, mat_z, levels = levels, cmap="Reds")
    cp = ax2.contour(mat_x, mat_y, mat_z, levels = levels, colors="k")
    ax2.clabel(cp, inline=1, fontsize=10, colors = "black")
    xlabel("x")
    ylabel("y")
    PyPlot.title(title * " contour")
    #tight_layout()
end
function wrap_quiver(x::Pos, y::Pos, np::Npos; title = "Quiver")
    if (size(x)[1], size(y)[1]) != size(np)
        error("incorrect dims")
    end
    clf()
    fig = figure("pyplot_quiverplot", figsize = (10, 10))
    q = quiver(x', y', slice_1(np)', slice_2(np)')
    maxlen = maximum(len.(replace(np, [NaN, NaN]=>[0.0]))) |> round_sig
    ax = gca()
    ax.quiverkey(q, X=0.8, Y = 0.05, U = maxlen, coordinates="figure", label="Quiver key, length = $maxlen",labelpos = "E")
    PyPlot.title(title)
end
function wrap_stream_a(x::Matrix, y::Matrix, np::Npos; density = 1., title = "Stream a")
    clf()
    fig = figure("pyplot_streamplot", figsize = (10,10))
    streamplot(x = x, y = y, u = slice_1(np), v = slice_2(np),
            density = density)
    PyPlot.title(title)
end
function wrap_stream_b(x::Matrix, y::Matrix, np::Npos; density = 1., title = "Stream b")
    clf()
    fig = figure("pyplot_streamplot", figsize = (10,10))
    vl = len.(np)
    streamplot(x = x, y = y, u = slice_1(np), v = slice_2(np),
             color = vl,
             linewidth = 2,
             cmap = PyPlot.cm.autumn,
             density = density)
    PyPlot.title(title)
end
function wrap_stream_c(x::Matrix, y::Matrix, np::Npos; density = 1., title = "Stream c")
    clf()
    fig = figure("pyplot_streamplot", figsize = (10,10))
    vl = len.(np)
    maxlen = maximum(len.(replace(np, [NaN, NaN]=>[0.0])))
    lw = 5 .* vl ./ maxlen
    streamplot(x = x, y = y, u = slice_1(np), v = slice_2(np),
        color = "k",
        density = density,
        linewidth = lw)
    PyPlot.title(title)
end
function wrap_stream_d(x::Matrix, y::Matrix, np::Npos; density = 1., title = "Stream d")
    clf()
    fig = figure("pyplot_streamplot", figsize = (10,10))
    vl = len.(np)
    maxlen = maximum(len.(replace(np, [NaN, NaN]=>[0.0])))
    lw = 5 .* vl ./ maxlen
    streamplot(x = x, y = y, u = slice_1(np), v = slice_2(np),
        color = vl,
        density = density,
        cmap = PyPlot.cm.YlGnBu,
        linewidth = lw)
    PyPlot.title(title)
end

max_stars_repo_name: hustf/Gravityfield

max_issues_repo_name: hustf/Gravityfield

max_forks_repo_name: hustf/Gravityfield


text:
using WCharUTF8
using Base.Test

include("wcharutf8_tests.jl")


text:
using Gadfly
using DataFrames, RDatasets

function linescatter()
    x = [1, 2, 3, 4]
    y = [10, 15, 13, 17]
    plot(x=x,y=y, Geom.line)
end
linescatter()

function multiple_scatter_traces()
    x=1:4
    l1 = layer(x = x, y = [10, 15, 13, 17],  Geom.point)
    l2 = layer(x = x, y = [16, 5, 11, 9], Geom.line)
    l3 = layer(x = x, y = [12, 9, 15, 12], Geom.point, Geom.line)
    l4 = layer(x = x, y = [5, 10, 8, 12], Geom.line)
    plot(l1, l2, l3, l4)
end
multiple_scatter_traces()

function data_labels()
    x = 1:5
    y1 = [1, 6, 3, 6, 1]
    y2 = [4, 1, 7, 1, 4]
    y1_labels = ["A-1", "A-2", "A-3", "A-4", "A-5"]
    y2_labels = ["B-a", "B-b", "B-c", "B-d", "B-e"]

    l1 = layer(x = x, y = y1, label=y1_labels, Geom.label, Geom.point)
    l2 = layer(x = x, y = y2, label=y2_labels, Geom.label, Geom.point)
    plot(l1,l2)
end
data_labels()

function area1()
end
area1()

function matrix_subplots()

end

matrix_subplots()

function advanced_layouts()

end
advanced_layouts()

function grouped_bar_example()

end

grouped_bar_example()

function stacked_bar_example()

end
stacked_bar_example()

function two_hists()

end

two_hists()

function box_plot()

end
box_plot()

function data_frame_scatter()
    iris = dataset("datasets", "iris");
    display(head(iris))
end

data_frame_scatter()


text:
using MPI
using MPIAsyncPools
using Test

MPI.Init()
const comm = MPI.COMM_WORLD
const root = 0
const rank = MPI.Comm_rank(comm)
const nworkers = MPI.Comm_size(comm) - 1
isroot() = MPI.Comm_rank(comm) == root
const data_tag = 0
const control_tag = 1

function shutdown(pool::MPIAsyncPool)
    for i in pool.ranks
        MPI.Isend(zeros(1), i, control_tag, comm)
    end
end

function root_main()
    pool = MPIAsyncPool(nworkers)
    @test pool.ranks == collect(1:nworkers)

    sendbuf = Vector{Float64}(undef, 1)
    isendbuf = zeros(eltype(sendbuf), nworkers*length(sendbuf))
    recvbuf = Vector{Float64}(undef, 3*nworkers)
    recvbufs = [view(recvbuf, (i-1)*3+1:i*3) for i in 1:nworkers]
    irecvbuf = copy(recvbuf)
    nwait = 2

    # setup a callback to test that it's called once per receive
    recvs = Set{Int}()
    function recvf(i, epoch, repoch, recvbuf)
        buf = reinterpret(Float64, recvbuf)
        @test length(buf) == 3
        rank, t, sepoch = buf
        @test i == rank
        @test sepoch == repoch
        if repoch != epoch # ignore state results
            return
        end
        if i in recvs
            error("recvf was called twice for worker $i")
        end
        push!(recvs, i)
        return
    end

    # test that we have at least nwait responses from the current epoch in each iteration
    for epoch in 1:100

        # empty the set for each iteration
        empty!(recvs)

        sendbuf[1] = epoch
        delay = @elapsed begin
            repochs = asyncmap!(pool, sendbuf, recvbuf, isendbuf, irecvbuf, comm; nwait, tag=data_tag, recvf)
        end
        from_this_epoch = 0
        for i in 1:nworkers
            wrank, t, wepoch = recvbufs[i]

            # no tests if we've never received from this worker
            if repochs[i] == 0
                continue
            end
            if repochs[i] == epoch
                from_this_epoch += 1
            end

            # test that workers echo what was sent to them
            @test wepoch == repochs[i]
        end
        @test from_this_epoch >= nwait

        # test that the callback was called once per receive
        @test length(recvs) == from_this_epoch
        for i in 1:nworkers
            if repochs[i] == epoch
                @test i in recvs
            end
        end
    end

    # test that all workers have returned after calling waitall!
    for epoch in 1:100
        repochs = asyncmap!(pool, sendbuf, recvbuf, isendbuf, irecvbuf, comm; nwait=1, tag=data_tag)
        repochs = waitall!(pool, recvbuf, irecvbuf)
        @test all(pool.active .== false)
    end

    # test using a custom function to determine when asyncmap! should return
    # (in this case that the first worker has returned)
    f = (epoch, repochs) -> repochs[1] == epoch
    for epoch in 101:200
        delay = @elapsed begin
            repochs = asyncmap!(pool, sendbuf, recvbuf, isendbuf, irecvbuf, comm; nwait=f, tag=data_tag)
        end
        @test repochs[1] == pool.epoch
    end
    shutdown(pool)
end

function worker_main()
    recvbuf = Vector{Float64}(undef, 1)
    sendbuf = Vector{Float64}(undef, 3) # [rank, t, epoch]
    sendbuf[1] = rank
    crreq = MPI.Irecv!(zeros(1), root, control_tag, comm)

    t = 0
    while true
        t += 1

        rreq = MPI.Irecv!(recvbuf, root, data_tag, comm)
        index, _ = MPI.Waitany!([crreq, rreq])
        if index == 1 # exit message on control channel
            break
        end

        epoch = recvbuf[1]
        sendbuf[2] = t
        sendbuf[3] = epoch
        sleep(max(rand()/10, 0.005))

        sreq = MPI.Isend(sendbuf, root, data_tag, comm)
    end
end

if isroot()
    root_main()
else
    worker_main()
end
MPI.Barrier(comm)
if isroot()
    print("All tests pass\n")
end


text:
function threshold(thres::AbstractThreshold, targets, scores)
    return find_threshold(thres, cpu(targets), cpu(scores))[1]
end

@adjoint function threshold(thres::AbstractThreshold, targets, scores)
    t, Δt_s = find_threshold(thres, cpu(targets), cpu(scores))
    t = convert(eltype(scores), t)
    Δt_s = convert(typeof(scores), Δt_s)
    return t, Δ -> (nothing, nothing, Δ .* Δt_s)
end

# methods
quant(surrogate, x, τ) = sum(surrogate, x)/length(x) - τ

struct PatMat <: AbstractThreshold
    τ::Real
    surrogate::Function
end

function find_threshold(thres::PatMat, targets, scores)
    @unpack surrogate, τ = thres

    t = [find_root(t -> quant(surrogate, scores .- t, τ), quantile(vec(scores), 1 - τ))]
    ps = params(scores, t)
    ∇l = gradient(() -> quant(surrogate, scores .- t, τ), ps)

    return t[1], .- ∇l[scores] ./ ∇l[t]
end

struct PatMatNP <: AbstractThreshold
    τ::Real
    surrogate::Function
end

function find_threshold(thres::PatMatNP, targets, scores)
    @unpack surrogate, τ = thres
    inds = find_negatives(targets)
    s = @views scores[inds]

    t = [find_root(t -> quant(surrogate, s .- t, τ), quantile(vec(s), 1 - τ))]
    ps = params(s, t)
    ∇l = gradient(() -> quant(surrogate, s .- t, τ), ps)

    ∇l_s = zero(scores)
    ∇l_s[inds] .= ∇l[s]

    return t[1], .- ∇l_s ./ ∇l[t]
end

struct TopPush <: AbstractThreshold end

function find_threshold(::TopPush, targets, scores)
    inds = find_negatives(targets)
    s = @views scores[inds]

    t, ind = findmax(s)
    Δt_s = zero(scores)
    Δt_s[inds[ind]] = 1

    return t, Δt_s
end

struct TopPushK <: AbstractThreshold
    K::Int
end

function find_threshold(thres::TopPushK, targets, scores)
    @unpack K = thres
    inds = find_negatives(targets)
    s = @views scores[inds]

    ind = partialsortperm(s, 1:K; rev = true)
    Δt_s = zero(scores)
    Δt_s[inds[ind]] .= 1/K

    return sum(s[ind])/K, Δt_s
end

struct τFPL <: AbstractThreshold
    τ::Real
end

function find_threshold(thres::τFPL, targets, scores)
    @unpack τ = thres
    inds = find_negatives(targets)
    s = scores[inds]
    K = floor(Int64, τ*length(inds)) + 1

    ind = partialsortperm(s, 1:K; rev = true)
    Δt_s = zero(scores)
    Δt_s[inds[ind]] .= 1/K

    return sum(s[ind])/K, Δt_s
end

struct TopMean <: AbstractThreshold
    τ::Real
end

function find_threshold(thres::TopMean, targets, scores)
    @unpack τ = thres
    K = floor(Int64, τ*length(scores)) + 1

    ind = partialsortperm(vec(scores), 1:K; rev = true)
    Δt_s = zero(scores)
    Δt_s[ind] .= 1/K

    return sum(scores[ind])/K, Δt_s
end

struct Grill <: AbstractThreshold
    τ::Real
end

function find_threshold(thres::Grill, targets, scores)
    @unpack τ = thres
    t = quantile(vec(scores), 1 - τ)
    Δt_s = zero(scores)

    return t, Δt_s
end

struct GrillNP <: AbstractThreshold
    τ::Real
end

function find_threshold(thres::GrillNP, targets, scores)
    @unpack τ = thres
    inds = find_negatives(targets)
    t = quantile(scores[inds], 1 - τ)
    Δt_s = zero(scores)

    return t, Δt_s
end


text:
module SlurpOptimizations
using  ...PatternStructure.Trees
using  ...PatternStructure.Checks
using  ...PatternStructure.SlurpTypes
using  ...PatternStructure.Special
export optimize_slurps!

#-----------------------------------------------------------------------------------
# dispatch function
#-----------------------------------------------------------------------------------

function optimize_slurps!(node, state)
  slurps = reverse(find(is_slurp, node.children))
  optimize_slurps!(node, slurps, [], state)
end

function optimize_slurps!(node, slurps, after, state)
  isempty(slurps) && return

  i = slurps[1]
  slurp = node.children[i]

  if isempty(after) # last slurp
    if is_simple_slurp(slurp, state)
       node.children[i] = make_simple_last_slurp(slurp, node, i)

    elseif is_simple_alternating_slurp(slurp, state)
      # simple alternating last slurp | (..., *{a,b,c..})
    end
  else
    gates  = find(x->isa(x, PatternGate), node.children[i+1:after[1]-1])
    equals = filter(j->isa(node.children[i+j].check, EqualityCheck), gates)

    if !isempty(equals)
      eqvalues = map(j->(node.children[i+j].check.value, i+j), equals)

      if is_simple_slurp(slurp, state)
        until, index = eqvalues[1]
        node.children[i] = make_simple_slurp_until_one(slurp, node, until, index-i)

      elseif is_simple_alternating_slurp(slurp, state)
      	# alternating slurp until X | (..., *{a,b,c...}, X, ...)
      end
    end
  end
  push!(after, i)
  optimize_slurps!(node, slurps[2:end], after, state)
end

#-----------------------------------------------------------------------------------
# maker functions
#-----------------------------------------------------------------------------------

make_slurp(slurp, head) =
  PatternNode(head, SlurpStep(), slurp.children, slurp.bindings, slurp.depth)

make_simple_last_slurp(slurp, node, i) =
  make_slurp(slurp, SimpleLastSlurp(length(node.children) - i))

make_simple_slurp_until_one(slurp, node, until, index) =
  make_slurp(slurp,
              isa(slurp.head, GreedySlurp) ?
                SimpleGreedySlurpUntil([until], index) :
                SimpleLazySlurpUntil([until], index))

#-----------------------------------------------------------------------------------
# conditions
#-----------------------------------------------------------------------------------

function is_simple_slurp(slurp, state)
  length(slurp.children) == 1 &&
  is_binding(slurp.children[1]) &&
  !is_const_binding(slurp.children[1], state)
end

function is_simple_alternating_slurp(slurp, state)
  length(slurp.children) > 1 &&
  all(is_binding, slurp.children) &&
  !is_const_binding(slurp.children[1], state)
end

function is_binding(node::PatternGate)
  isa(node.check, Binding) &&
  isa(node.child, PatternLeaf)
end
is_binding(x) = false

function is_const_binding(node::PatternGate, state)
  node.check.name in state.consts
end

end


text:
module ContinuousWrapper

using JavaCall

const AbstractRealDistribution = @jimport org.hipparchus.distribution.continuous.AbstractRealDistribution
const BetaDistribution = @jimport org.hipparchus.distribution.continuous.BetaDistribution
const CauchyDistribution = @jimport org.hipparchus.distribution.continuous.CauchyDistribution
const ChiSquaredDistribution = @jimport org.hipparchus.distribution.continuous.ChiSquaredDistribution
const Class = @jimport java.lang.Class
const ConstantRealDistribution = @jimport org.hipparchus.distribution.continuous.ConstantRealDistribution
const EnumeratedRealDistribution = @jimport org.hipparchus.distribution.continuous.EnumeratedRealDistribution
const ExponentialDistribution = @jimport org.hipparchus.distribution.continuous.ExponentialDistribution
const FDistribution = @jimport org.hipparchus.distribution.continuous.FDistribution
const GammaDistribution = @jimport org.hipparchus.distribution.continuous.GammaDistribution
const GumbelDistribution = @jimport org.hipparchus.distribution.continuous.GumbelDistribution
const JString = @jimport java.lang.String
const LaplaceDistribution = @jimport org.hipparchus.distribution.continuous.LaplaceDistribution
const LevyDistribution = @jimport org.hipparchus.distribution.continuous.LevyDistribution
const List = @jimport java.util.List
const LogNormalDistribution = @jimport org.hipparchus.distribution.continuous.LogNormalDistribution
const LogisticDistribution = @jimport org.hipparchus.distribution.continuous.LogisticDistribution
const NakagamiDistribution = @jimport org.hipparchus.distribution.continuous.NakagamiDistribution
const NormalDistribution = @jimport org.hipparchus.distribution.continuous.NormalDistribution
const Object = @jimport java.lang.Object
const ParetoDistribution = @jimport org.hipparchus.distribution.continuous.ParetoDistribution
const RealDistribution = @jimport org.hipparchus.distribution.RealDistribution
const TDistribution = @jimport org.hipparchus.distribution.continuous.TDistribution
const TriangularDistribution = @jimport org.hipparchus.distribution.continuous.TriangularDistribution
const UniformRealDistribution = @jimport org.hipparchus.distribution.continuous.UniformRealDistribution
const WeibullDistribution = @jimport org.hipparchus.distribution.continuous.WeibullDistribution

export AbstractRealDistribution
export BetaDistribution
export CauchyDistribution
export ChiSquaredDistribution
export ConstantRealDistribution
export EnumeratedRealDistribution
export ExponentialDistribution
export FDistribution
export GammaDistribution
export GumbelDistribution
export LaplaceDistribution
export LevyDistribution
export LogNormalDistribution
export LogisticDistribution
export NakagamiDistribution
export NormalDistribution
export ParetoDistribution
export TDistribution
export TriangularDistribution
export UniformRealDistribution
export WeibullDistribution
export cumulative_probability
export density
export equals
export get_alpha
export get_beta
export get_class
export get_degrees_of_freedom
export get_denominator_degrees_of_freedom
export get_location
export get_mean
export get_median
export get_mode
export get_numerator_degrees_of_freedom
export get_numerical_mean
export get_numerical_variance
export get_pmf
export get_scale
export get_shape
export get_standard_deviation
export get_support_lower_bound
export get_support_upper_bound
export hash_code
export inverse_cumulative_probability
export is_support_connected
export log_density
export notify
export notify_all
export probability
export to_string
export wait

include("abstract_real_distribution.jl")
include("beta_distribution.jl")
include("cauchy_distribution.jl")
include("chi_squared_distribution.jl")
include("constant_real_distribution.jl")
include("enumerated_real_distribution.jl")
include("exponential_distribution.jl")
include("f_distribution.jl")
include("gamma_distribution.jl")
include("gumbel_distribution.jl")
include("laplace_distribution.jl")
include("levy_distribution.jl")
include("log_normal_distribution.jl")
include("logistic_distribution.jl")
include("nakagami_distribution.jl")
include("normal_distribution.jl")
include("pareto_distribution.jl")
include("t_distribution.jl")
include("triangular_distribution.jl")
include("uniform_real_distribution.jl")
include("weibull_distribution.jl")

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.dag" begin
    graph = plaindigraph(10, [
        (1, 4), (2, 4), (3, 7), (4, 5), (4, 6),
        (4, 7), (6, 9), (7, 8), (7, 9), (7, 10)
    ])
    @test issetequal(ancestors(graph, 7), [1, 2, 3, 4])
    @test issetequal(descendants(graph, 4), [5, 6, 7, 8, 9, 10])
    nodes = topologicalsort(graph)
    @test nodes[4] == 4
    @test nodes[end] in [5, 8, 9, 10]
end # graph.shortestpath


text:
server = nothing
port = 8000

function start_server(_port::Int64)
	log_message("starting up server...")
	global port = _port
	global server = Sockets.listen(Sockets.InetAddr(Sockets.localhost, port))
    # global server = HTTP.Servers.Server(create_base_router())
    try
        HTTP.serve(create_base_router(), "127.0.0.1", port, verbose = true, server = server)
    finally
        close(server)
	end
end

stop_server() = close(server)


text:
using DelimitedFiles, Statistics

"""

    housing([test]; [url, file])

Return (xtrn,ytrn,xtst,ytst) from the [UCI Boston
Housing](https://archive.ics.uci.edu/ml/machine-learning-databases/housing)
dataset The dataset has housing related information for 506
neighborhoods in Boston from 1978. Each neighborhood has 14
attributes, the goal is to use the first 13, such as average number of
rooms per house, or distance to employment centers, to predict the
14’th attribute: median dollar value of the houses.

`test=0` by default and `xtrn` (13,506) and `ytrn` (1,506) contain the
whole dataset. Otherwise data is shuffled and split into train and
test portions using the ratio given in `test`. xtrn and xtst are
always normalized by subtracting the mean and dividing into standard
deviation.

"""
function housing(test=0.0;
                 file=joinpath(@__DIR__, "housing", "housing.data"),
                 url="https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data")
    if !isfile(file)
        isdir(dirname(file)) || mkpath(dirname(file))
        @info("Downloading $url to $file")
        download(url, file)
    end
    data = readdlm(file)'
    # @show size(data) # (14,506)
    x = data[1:13,:]
    y = data[14:14,:]
    x = (x .- mean(x,dims=2)) ./ std(x,dims=2) # Data normalization
    if test == 0
        xtrn = xtst = x
        ytrn = ytst = y
    else
        r = randperm(size(x,2))          # trn/tst split
        n = round(Int, (1-test) * size(x,2))
        xtrn=x[:,r[1:n]]
        ytrn=y[:,r[1:n]]
        xtst=x[:,r[n+1:end]]
        ytst=y[:,r[n+1:end]]
    end
    return (xtrn, ytrn, xtst, ytst)
end

nothing


text:
type ForwardNet{T<:Real}
    dag::DiGraph # graph over nodes
    nodes::Vector{Node{T}}
    name_to_index::Dict{Symbol,Int} # Symbol → index in dag and nodes (note: not all Nodes have names)
    name_counter::Int

    ForwardNet{T}() where T = new(DiGraph(0), Node[], Dict{Symbol,Int}(), 0)
end

Base.getindex(net::ForwardNet, index::Int) = net.nodes[index]
Base.getindex(net::ForwardNet, name::Symbol) = net.nodes[net.name_to_index[name]]
indexof(net::ForwardNet, name::Symbol) = net.name_to_index[name]
lastindex(net::ForwardNet) = nv(net.dag)
function Base.push!{T}(net::ForwardNet{T}, node::Node{T}, parents::Vector{Int}=Int[])

    add_vertex!(net.dag)
    push!(net.nodes, node)
    i = length(net.nodes)

    nodename = ForwardNets.name(node)
    @assert(!haskey(net.name_to_index, nodename))
    net.name_to_index[nodename] = i

    for parent in parents
        add_edge!(net.dag, parent, i)
    end

    net
end
function Base.push!{T}(net::ForwardNet{T}, node::Node{T}, parents::Union{Vector{Symbol}, Vector{NameOrIndex}})
    int_parents = Array{Int}(length(parents))
    for (i,p) in enumerate(parents)
        if isa(p, Int)
            int_parents[i] = p
        elseif isa(p, Symbol)
            int_parents[i] = net.name_to_index[p]
        end
    end

    push!(net, node, int_parents)
end
Base.push!{T}(net::ForwardNet{T}, node::Node{T}, parent::Int) = push!(net, node, Int[parent])
Base.push!{T}(net::ForwardNet{T}, node::Node{T}, parent::Symbol) = push!(net, node, Int[net.name_to_index[parent]])

function next_name_number!(net::ForwardNet)
    n = net.name_counter += 1
    n
end
create_next_name!(net::ForwardNet, prefix::String) = convert(Symbol, prefix * "_" * string(next_name_number!(net)))


text:
"""
    TimeStepStrategy

Abstract type for strategies for updating the time step with
[`update_dτ()`](@ref). Implemented strategies:

   * [`ConstantTimeStep`](@ref)
"""
abstract type TimeStepStrategy end

"""
    ConstantTimeStep <: TimeStepStrategy

Keep `dτ` constant.
"""
struct ConstantTimeStep <: TimeStepStrategy end

"""
    update_dτ(s<:TimeStepStrategy, dτ, tnorm) -> new dτ
Update the time step according to the strategy `s`.
"""
update_dτ(::ConstantTimeStep, dτ, args...) = dτ


text:
function open_imagefile(filename::AbstractString)
    if is_apple()
        run(`open $(filename)`)
    elseif is_linux() || is_bsd()
        run(`xdg-open $(filename)`)
    elseif is_windows()
        run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    else
        warn("Showing plots is not supported on OS $(string(Compat.KERNEL))")
    end
end

max_stars_repo_name: wakakusa/ScatterMatrixPlots

max_issues_repo_name: wakakusa/ScatterMatrixPlots

max_forks_repo_name: wakakusa/ScatterMatrixPlots