text:
# ] Gadfly@1.2.1 TreeParzen MLJ DataFrames XGBoost MLJModels MLJTuning CategoricalArrays ComputationalResources

import Gadfly, MLJ, DataFrames, MLJTuning, MLJModels, CategoricalArrays, ComputationalResources

using TreeParzen

import Random: seed!

Gadfly.set_default_plot_size(20Gadfly.cm, 16Gadfly.cm)

# this isn't great but it will do
seed!(1999)

# Helper function for turning strings into not strings,
# seems to be required for use with XGBoostRegressor,
# Ints and other things make it complain
conv(x::CategoricalArrays.CategoricalArray) = Float64.(MLJ.int(x))
conv(x) = Float64.(x)

function opt_hist_plot(tuning_history, title, path)

    metric_name = String(Symbol(first(first(tuning_history).measure)))
    metric = first.(getfield.(tuning_history, :measurement))
    cummetric = accumulate(min, metric)

    upper_bound = min(minimum(metric) * 3, maximum(metric))

    plotobj = Gadfly.plot(
        Gadfly.layer(x=1:length(metric), y=metric, Gadfly.Geom.line, Gadfly.Theme(default_color=Gadfly.colorant"orange")),
        Gadfly.layer(x=1:length(metric), y=cummetric, Gadfly.Geom.step),
        Gadfly.Guide.ylabel("MLJ Tuning optimisation measurement: $metric_name"; orientation=:vertical),
        Gadfly.Guide.xlabel("Sequence"),
        Gadfly.Coord.Cartesian(ymax=upper_bound),
        Gadfly.Guide.Title(title),
    )

    plotobj |> Gadfly.SVGJS(path)
    Gadfly.display(plotobj)
end

#### constants ####
NUM_CV_FOLDS = 4
PCT_TRAIN_DATA = 0.75
NUM_TP_ITER_SMALL = 25
NUM_TP_ITER_LARGE = 250
#### end constants ####

MLJ.@load XGBoostRegressor

Features, targets = MLJ.@load_reduced_ames
# This one turns strings into not strings, and everything into Float64s
IntCat_Features = NamedTuple{keys(Features)}(conv.(values(deepcopy(Features))))
Features = DataFrames.DataFrame(Features)
IntCat_Features = DataFrames.DataFrame(IntCat_Features)

# Do hold-out partitioning. If you want same results each time use shuffle=false or set RNG seed
train, test = MLJ.partition(eachindex(targets), PCT_TRAIN_DATA, shuffle=true)

train_features = IntCat_Features[train, :]
train_targets = targets[train]

test_features = IntCat_Features[test, :]
test_targets = targets[test]

num_train_data = length(train_targets)
training_data_per_fold = (num_train_data / NUM_CV_FOLDS) * (NUM_CV_FOLDS - 1)

# Some prior decisions (almost arbitrary, there aren't strong reasons to make these decisions)
#
# - Test holdout : 25%
# - 4 fold cross validation -- this leaves each individual training set with ~820 data points and ~200 for evaluation on each fold
# - Optimising using MAE as the metric because we will assess final quality on RMSL (and because for this example we're going to optimise directly in logspace)

# We identified 7 potentially interesting optimisation parameters for gradient boosted trees.
# One of the first things we note is that even if we just selected 2 points on the relevant ranges,
# a cartesian product (grid search) over this space would be 128 evaluations
# Search over 3 points on each axis (say, min, median, max) would take the search to well over 2000 evaluations.
space = Dict(
    :num_round => HP.QuantUniform(:num_round, 1., 500., 1.),
    :eta => HP.LogUniform(:eta, -3., 0.),
    :gamma => HP.LogUniform(:gamma, -3., 3.),
    :max_depth => HP.QuantUniform(:max_depth, 1., ceil(log2(training_data_per_fold)), 1.0),
    :min_child_weight => HP.LogUniform(:min_child_weight, -5., 2.),
    :lambda => HP.LogUniform(:lambda, -5., 2.),
    :alpha => HP.LogUniform(:alpha, -5., 2.),
)

model_tpl = XGBoostRegressor()

# Here we are tuning the model by MAE in logspace, which isn't necessarily
# right but we also need to ensure only positive valued outputs are produced.
# And this is easier for TreeParzen demonstration purposes than constructing a full
# learning network which log targets and exp predictions
tuning = MLJTuning.TunedModel(
    model=model_tpl,
    ranges=space,
    tuning=TreeParzen.MLJTreeParzen.MLJTreeParzenTuning(),
    n=NUM_TP_ITER_SMALL,
    resampling=MLJ.CV(nfolds=NUM_CV_FOLDS),
    measure=MLJ.mav,
)

mach = MLJ.machine(tuning, train_features, log.(train_targets))
MLJ.fit!(mach)

# julia> fit!(mach)
# [ Info: Training Machine{DeterministicTunedModel{MLJTreeParzenTuning,…}} @ 5…58.
# [ Info: Attempting to evaluate 25 models.
# Evaluating over 20 metamodels: 100%[=========================] Time: 0:03:21
# Evaluating over 1 metamodels: 100%[=========================] Time: 0:00:07
# Evaluating over 1 metamodels: 100%[=========================] Time: 0:00:02
# Evaluating over 1 metamodels: 100%[=========================] Time: 0:00:14
# Evaluating over 1 metamodels: 100%[=========================] Time: 0:00:05
# Evaluating over 1 metamodels: 100%[=========================] Time: 0:00:08
# Machine{DeterministicTunedModel{MLJTreeParzenTuning,…}} @ 5…58
#
# julia>

# perform the evaluation(s) -- predict for a TunedModel will use best one (or best parameters trained on whole data, depending on settings)
pred = exp.(MLJ.predict(mach, test_features))
@show MLJ.rmsl(test_targets, pred)

best_model = MLJ.fitted_params(mach).best_model

for x in keys(space) println("$x = $(getproperty(best_model, x))") end

opt_hist_plot(mach.report.history, "Simple optimisation of tree boosting", joinpath(@__DIR__, "../../../examples/simple_mlj_demo/images/simple_tree_tuning.svg"))

# To demonstrate use of suggestions, we can take the best result from last tuning run.
# BEAR IN MIND that this is cheating from a DS perspective, this is just to demonstrate the functionality.
suggestion = Dict(key => getproperty(best_model, key) for key in keys(space))

search = TreeParzen.MLJTreeParzen.MLJTreeParzenSpace(space, suggestion)

tuning = MLJTuning.TunedModel(
    model=model_tpl,
    ranges=search,
    tuning=TreeParzen.MLJTreeParzen.MLJTreeParzenTuning(;random_trials=3),
    n=NUM_TP_ITER_SMALL,
    resampling=MLJ.CV(nfolds=NUM_CV_FOLDS),
    measure=MLJ.mav,
)

mach = MLJ.machine(tuning, train_features, log.(train_targets))
MLJ.fit!(mach)

# perform the evaluation(s) -- predict for a TunedModel will use best one (or best parameters trained on whole data, depending on settings)
pred = exp.(MLJ.predict(mach, test_features))
@show MLJ.rmsl(test_targets, pred)

best_model = MLJ.fitted_params(mach).best_model

for x in keys(space) println("$x = $(getproperty(best_model, x))") end

# we can also accelerate learning by using parallelism
# notice how above we had 20 metamodels followed by sequences of 1?
# What is happening here is that there is an initial stage (by default 20)
# of models drawn entirely at random (using specified prior distributions)
# without probabilistic modelling. Once probabilistic modelling kicks in,
# we draw a suggestion and then update with result. This doesn't seem amenable
# to parallelism, and it isn't. However there is an additional parameter
# `max_simultaneous_draws` which allows the system to draw `n` samples before
# `updating the distribution. Whilst intuitively this allows parallelism, it
# also enables more exploration before drawing again according to updated distribution.
# If using this parameter, consider increasing the value of `linear_forgetting`
# from its default of 25 to a higher number -- a good place is probably at least `n*25`. The
# `linear_forgetting` parameter keeps most n recently observed results equally weighted
# and older observations are  weighted by a linear ramp according to their age.

# grab previous best result again, hey we started cheating, might as well continue.
suggestion = Dict(key => getproperty(best_model, key) for key in keys(space))
search = TreeParzen.MLJTreeParzen.MLJTreeParzenSpace(space, suggestion)

tuning = MLJTuning.TunedModel(
    model=model_tpl,
    ranges=space,
    tuning=TreeParzen.MLJTreeParzen.MLJTreeParzenTuning(;random_trials=3, max_simultaneous_draws=2, linear_forgetting=50),
    n=NUM_TP_ITER_SMALL,
    resampling=MLJ.CV(nfolds=NUM_CV_FOLDS),
    measure=MLJ.mav,
)

mach = MLJ.machine(tuning, train_features, log.(train_targets))
MLJ.fit!(mach)

# julia> fit!(mach)
# [ Info: Training Machine{DeterministicTunedModel{MLJTreeParzenTuning,…}} @ 4…10.
# [ Info: Attempting to evaluate 25 models.
# Evaluating over 3 metamodels: 100%[=========================] Time: 0:00:35
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:07
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:00
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:50
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:12
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:24
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:23
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:03
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:16
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:03
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:23
# Evaluating over 2 metamodels: 100%[=========================] Time: 0:00:08
# Machine{DeterministicTunedModel{MLJTreeParzenTuning,…}} @ 4…10

# It's worth noting, that this changed behaviour of optimiser but didn't
# introduce parallelism directly: That's up to MLJ. So lets parallelise.
# It wouldn't really accelerate computations in this case possibly due to
# overheads, but is an option for larger tasks. Also note that we can instead
# parallelise via `resampling_acceleration` if we don't want to alter behaviour of
# optimiser, but this might be more limited.

suggestion = Dict(key => getproperty(best_model, key) for key in keys(space))
search = TreeParzen.MLJTreeParzen.MLJTreeParzenSpace(space, suggestion)

tuning = MLJTuning.TunedModel(
    model=model_tpl,
    ranges=space,
    tuning=TreeParzen.MLJTreeParzen.MLJTreeParzenTuning(;random_trials=3, max_simultaneous_draws=2, linear_forgetting=50),
    n=NUM_TP_ITER_SMALL,
    resampling=MLJ.CV(nfolds=NUM_CV_FOLDS),
    measure=MLJ.mav,
    acceleration=ComputationalResources.CPUProcesses(),
)

mach = MLJ.machine(tuning, train_features, log.(train_targets))
MLJ.fit!(mach)

# An interesting feature of note is that TreeParzen supports "tree-structured"
# parameter spaces; hence the name. Originally conceived for optimising
# Deep Belief Networks (DBNs), where parameters take form such as number of layers,
# and how many nodes are within each layer. Here, we can see the idea that if
# we have 2 layers, then number of nodes in layer 3 as a parameter is not relevant.
# We can construct a similar example using XGBoost - either we can use boosted trees
# or boosted linears. Some parameters are relevant in both cases, but even so we
# would want to model them differently, e.g. `num_iterations` might need to be lower for
# boosting trees than for boosting linear functions.

# First, create a container model (because of space nesting, we need to do this)
# plus its constructor:

mutable struct tuned_xgb <: MLJ.Deterministic
    xgb::XGBoostRegressor
end
tuned_xgb(;xgb=Dict{Symbol, Any}()) = tuned_xgb(XGBoostRegressor(;xgb...))

# quick fit and predict methods
MLJ.fit(t::tuned_xgb, verbosity::Int, X, y, w=nothing) = MLJ.fit(t.xgb, verbosity, X, y, w)
MLJ.predict(t::tuned_xgb, fitted, X) = MLJ.predict(t.xgb, fitted, X)

# Define search parameters for the tree search space
tree_space = Dict(
    :booster => "gbtree",
    :num_round => HP.QuantUniform(:num_round_tree, 50., 750., 1.),
    :eta => HP.LogUniform(:eta_tree, -3., 0.),
    :gamma => HP.LogUniform(:gamma_tree, -3., 3.),
    :max_depth => HP.QuantUniform(:max_depth_tree, 1., ceil(log2(training_data_per_fold)), 1.0),
    :min_child_weight => HP.LogUniform(:min_child_weight_tree, -5., 1.),
    :lambda => HP.LogUniform(:lambda_tree, -5., 1.),
    :alpha => HP.LogUniform(:alpha_tree, -5., 1.),
)

# Define search parameters for the linear search space
linear_space = Dict(
    :booster => "gblinear",
    :updater => "coord_descent",
    :num_round => HP.QuantUniform(:num_round_linear, 500., 1000., 1.),
    :lambda => HP.LogUniform(:lambda_linear, -10., 0.),
    :alpha => HP.LogUniform(:alpha_linear, -10., 0.),
    :feature_selector => HP.Choice(:feature_selector_linear, ["cyclic", "greedy"]),
)

# Now we combine them so that it chooses either one search space or another (and then select parameters for each)
joint_space = Dict(
    :xgb => HP.Choice(:xgb, [linear_space, tree_space])
)

# Because we have a top level conditional, crank up the number of random trials to start off with
tuning = MLJTuning.TunedModel(
    model=tuned_xgb(),
    ranges=joint_space,
    tuning=TreeParzen.MLJTreeParzen.MLJTreeParzenTuning(;random_trials=50),
    n=NUM_TP_ITER_LARGE,
    resampling=MLJ.CV(nfolds=NUM_CV_FOLDS),
    measures=MLJ.mav,
)

# Do the usual fitting
mach = MLJ.machine(tuning, train_features, log.(train_targets))
MLJ.fit!(mach)

# Print out final evaluation, but we dont know what the best model is yet
pred = exp.(MLJ.predict(mach, test_features))
@show MLJ.rmsl(test_targets, pred)

# Take a look at the best performing model
best_model = MLJ.fitted_params(mach).best_model.xgb
if best_model.booster == "gbtree"
    println("Tree params")
    for x in keys(tree_space) println("$x = $(getproperty(best_model, x))") end
else
    println("Linear params")
    for x in keys(linear_space) println("$x = $(getproperty(best_model, x))") end
end

# Look at learning history
opt_hist_plot(mach.report.history, "Conditional optimisation of tree/linear boosting", joinpath(@__DIR__, "../../../examples/simple_mlj_demo/images/joint_linear_tree_tuning.svg"))

# Find times we selected linear:
linear_selected = getfield.(getfield.(first.(mach.report.history), :xgb), :booster) .== "gblinear"
linear_hist = mach.report.history[linear_selected]

opt_hist_plot(linear_hist, "History (partial) of evaluations conditioned on linear boosts", joinpath(@__DIR__, "../../../examples/simple_mlj_demo/images/joint_linear_only_iterations.svg"))

tree_hist = mach.report.history[.!linear_selected]

opt_hist_plot(tree_hist, "History (partial) of evaluations conditioned on tree boosts", joinpath(@__DIR__, "../../../examples/simple_mlj_demo/images/joint_tree_only_iterations.svg"))


text:
# This file was generated by JuDoc, do not modify it. # hide
vv = inverse_transform(stand, w)
@show sum(abs.(vv .- v))

ext: jl


text:
"""
    TypeV()

Lorem ipsum dolor sit amet.
"""
struct TypeV

end


text:
using SupervisedLearning
using UnicodePlots

import EmpiricalRisks
import Regression

d = 10
n = 100
p = 3
w = randn(d+1)
Xtr = randn(d, n)
Xte = randn(d, n)
for i = 2:d
  Xtr[i,:] =  Xtr[1,:] .^ i
  Xtr[i,:] = (Xtr[i,:] - mean(Xtr[i,:])) ./ std(Xtr[i,:])
  Xte[i,:] =  Xte[1,:] .^ i
  Xte[i,:] = (Xte[i,:] - mean(Xte[i,:])) ./ std(Xte[i,:])
end
tr = sign(Xtr[1:p,:]'w[1:p] + w[d+1] + 1.5 * randn(n))
te = sign(Xte[1:p,:]'w[1:p] + w[d+1])

trainSet = dataSource(Xtr, tr, SignedClassEncoding(["no","yes"]))
testSet = dataSource(Xte, te, SignedClassEncoding(["no","yes"]))

print(barplot(classDistribution(trainSet)..., title = "Trainset", width = 30))
print(barplot(classDistribution(testSet)..., title = "Testset", width = 30))

models = [("without Regularization", Classifier.LogisticRegression()),
          ("with L1 penalty", Classifier.LogisticRegression(l1_coef = .1)),
          ("with L2 penalty", Classifier.LogisticRegression(l2_coef = .1))]

solvers = [Solver.GradientDescent()]

for solver in solvers, (desc, model) in models

  @test state(model) == :uninitialized
  @test iterations(model) == 0
  @test_throws StateError trainingCurve(model)
  @test_throws StateError cost(model)

  #train!(model, data, solver, max_iter = 50, break_every = 5)
  maxiter = 1000
  train!(model, trainSet, solver, max_iter = maxiter, break_every = 5) do
    @test state(model) == :training
    @test iterations(model) % 5 == 0
    remember!(model, :testSet, cost(model, testSet))
  end

  @test state(model) == :trained
  @test 0 < iterations(model) <= maxiter
  @test_approx_eq cost(model) -fitness(model)
  @test_approx_eq cost(model, trainSet) -fitness(model, trainSet)
  @test_approx_eq cost(model) cost(model, trainSet)

  x, y = trainingCurve(model)#history(model, :tst)
  x, y = history(model, :testSet)
  plt = lineplot(x, y, ylim=[floor(minimum(y)), ceil(maximum(y))], xlim=[0, maxiter], width = 30, height = 2)
  annotate!(plt, :r, 1, "$(name(model)) (test acc: $(round(accuracy(model, testSet),2)))")
  annotate!(plt, :r, 2, "$(typeof(solver)) $desc")
  annotate!(plt, :bl, ""); annotate!(plt, :br, "")
  print(plt)

  y = predict(model, trainSet)
  #@test sum(y .== groundtruth(trainSet)) > .8 * nobs(trainSet)
  @test_approx_eq accuracy(model, trainSet) (sum(y .== groundtruth(trainSet)) / length(groundtruth(trainSet)))

  p = predictProb(model, trainSet)
  @test length(y) == length(p)
  @test maximum(p) <= 1.
  @test minimum(p) >= 0.
end


text:
include("dependencies_for_runtests.jl")
using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField, VectorInvariant

@testset "VerticalVorticityField with HydrostaticFreeSurfaceModel" begin

    for arch in archs
        @testset "VerticalVorticityField with HydrostaticFreeSurfaceModel [$arch]" begin
            @info "  Testing VerticalVorticityField with HydrostaticFreeSurfaceModel [$arch]..."

            grid = LatitudeLongitudeGrid(arch, size = (3, 3, 3),
                                         longitude = (0, 60),
                                         latitude = (15, 75),
                                         z = (-1, 0))

            model = HydrostaticFreeSurfaceModel(grid = grid,
                                                momentum_advection = VectorInvariant())

            ψᵢ(λ, φ, z) = rand()
            set!(model, u=ψᵢ, v=ψᵢ)

            ζ = VerticalVorticityField(model)

            compute!(ζ)

            @test all(isfinite.(ζ.data))
        end
    end
end


text:
#int ldap_is_ldap_url( const char *url )
function ldap_is_ldap_url(url::AbstractString)
    ccall((:ldap_is_ldap_url, libldap), Cint, (Cstring,), url)
end

# int ldap_initialize(LDAP **ldp, char *uri)
function ldap_initialize(ldp_handle_ref::Ref{Ptr{Cvoid}}, uri::AbstractString)
    ccall((:ldap_initialize, libldap), Cint, (Ref{Ptr{Cvoid}}, Cstring), ldp_handle_ref, uri)
end

#char *ldap_err2string( int err );
function ldap_err2string(err::Integer)
    ccall((:ldap_err2string, libldap), Cstring, (Cint,), err)
end

# int ldap_url_parse( const char *url, LDAPURLDesc **ludpp )
function ldap_url_parse(url::AbstractString, ldap_url_desc_handle_ref::Ref{Ptr{LDAPURLDesc}})
    ccall((:ldap_url_parse, libldap), Cint, (Cstring, Ref{Ptr{LDAPURLDesc}}), url, ldap_url_desc_handle_ref)
end

# void ldap_free_urldesc( LDAPURLDesc *ludp );
function ldap_free_urldesc(ldap_url_desc_handle::Ptr{LDAPURLDesc})
    ccall((:ldap_free_urldesc, libldap), Cvoid, (Ptr{LDAPURLDesc},), ldap_url_desc_handle)
end

# int ldap_get_option(LDAP *ld, int option, void *outvalue);
function ldap_get_option(ldp_handle::Ptr{Cvoid}, option::LDAPOption, out_value_ref::Ref{Cint})
    opt_int = Cint(option)
    ccall((:ldap_get_option, libldap), Cint, (Ptr{Cvoid}, Cint, Ref{Cint}), ldp_handle, opt_int, out_value_ref)
end

# int ldap_set_option(LDAP *ld, int option, const void *invalue);
function ldap_set_option(ldp_handle::Ptr{Cvoid}, option::LDAPOption, in_value::T) where {T}
    opt_int = Cint(option)
    in_value_ref = Ref(Cint(in_value))
    ccall((:ldap_set_option, libldap), Cint, (Ptr{Cvoid}, Cint, Ref{Cint}), ldp_handle, opt_int, in_value_ref)
end

# int ldap_simple_bind_s(LDAP *ld, const char *who, const char *passwd);
function ldap_simple_bind_s(ldp_handle::Ptr{Cvoid}, who::AbstractString, password::AbstractString)
    ccall((:ldap_simple_bind_s, libldap), Cint, (Ptr{Cvoid}, Cstring, Cstring), ldp_handle, who, password)
end

# int ldap_unbind_s(LDAP *ld);
function ldap_unbind_s(ldp_handle::Ptr{Cvoid})
    ccall((:ldap_unbind_s, libldap), Cint, (Ptr{Cvoid},), ldp_handle)
end

# int ldap_msgfree( LDAPMessage *msg );
function ldap_msgfree(msg_handle::Ptr{Cvoid})
    ccall((:ldap_msgfree, libldap), Cint, (Ptr{Cvoid},), msg_handle)
end

#=
int ldap_search_ext_s(
      LDAP *ld,
      char *base,
      int scope,
      char *filter,
      char *attrs[],
      int attrsonly,
      LDAPControl **serverctrls,
      LDAPControl **clientctrls,
      struct timeval *timeout,
      int sizelimit,
      LDAPMessage **res );
=#
function ldap_search_ext_s(
            ldp_handle::Ptr{Cvoid},
            base::AbstractString,
            scope::Integer,
            ldap_filter::Union{Nothing, AbstractString},
            attr_desc_only::Bool,
            size_limit::Integer,
            result::Ref{Ptr{Cvoid}}
        )
    ccall((:ldap_search_ext_s, libldap), Cint,
        (Ptr{Cvoid}, Cstring, Cint, Cstring, Ptr{Ptr{UInt8}}, Cint,
            Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ref{Ptr{Cvoid}}),
        ldp_handle, base, scope, ldap_filter == nothing ? C_NULL : ldap_filter, C_NULL, attr_desc_only, C_NULL, C_NULL, C_NULL, size_limit, result)
end

# int ldap_count_messages( LDAP *ld, LDAPMessage *result )
function ldap_count_messages(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_count_messages, libldap), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

# LDAPMessage *ldap_first_message( LDAP *ld, LDAPMessage *result )
function ldap_first_message(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_first_message, libldap), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

# LDAPMessage *ldap_next_message( LDAP *ld, LDAPMessage *message )
function ldap_next_message(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_next_message, libldap), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

# int ldap_msgtype( LDAPMessage *msg );
function ldap_msgtype(msg_handle::Ptr{Cvoid})
    ccall((:ldap_msgtype, libldap), Cint, (Ptr{Cvoid},), msg_handle)
end

# int ldap_msgid( LDAPMessage *msg );
function ldap_msgid(msg_handle::Ptr{Cvoid})
    ccall((:ldap_msgid, libldap), Cint, (Ptr{Cvoid},), msg_handle)
end

#ldap_count_references LDAP_P((
#    LDAP *ld,
#    LDAPMessage *chain ));
function ldap_count_references(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_count_references, libldap), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

#LDAP_F( int )
#ldap_count_entries LDAP_P((
#    LDAP *ld,
#    LDAPMessage *chain ));
function ldap_count_entries(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_count_entries, libldap), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

#char * ldap_get_dn LDAP_P((
#    LDAP *ld,
#    LDAPMessage *entry ));
function ldap_get_dn(ldp_handle::Ptr{Cvoid}, entry_handle::Ptr{Cvoid})
    ccall((:ldap_get_dn, libldap), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, entry_handle)
end

#LDAP_F( void )
#ldap_memfree LDAP_P((
#    void* p ));
function ldap_memfree(handle::Ptr{T}) where {T}
    ccall((:ldap_memfree, libldap), Cvoid, (Ptr{T},), handle)
end

#LDAP_F( LDAPMessage * )
#ldap_first_entry LDAP_P((
#    LDAP *ld,
#    LDAPMessage *chain ));
function ldap_first_entry(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_first_entry, libldap), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

#LDAP_F( LDAPMessage * )
#ldap_next_entry LDAP_P((
#    LDAP *ld,
#    LDAPMessage *entry ));
function ldap_next_entry(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_next_entry, libldap), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

#LDAP_F( int )
#ldap_get_entry_controls LDAP_P((
#  LDAP      *ld,
#  LDAPMessage   *entry,
#  LDAPControl   ***serverctrls));
function ldap_get_entry_controls(ldp_handle::Ptr{Cvoid}, entry_handle::Ptr{Cvoid}, controls_vec_ref::Ref{Ptr{Ptr{LDAPControl}}})
    ccall((:ldap_get_entry_controls, libldap), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{Ptr{Ptr{LDAPControl}}}), ldp_handle, entry_handle, controls_vec_ref)
end

# ldap_controls_free() frees a NULL-terminated array of controls.
#=
LDAP_F( void )
ldap_controls_free LDAP_P((
    LDAPControl **ctrls ));
=#
function ldap_controls_free(ctrls_handle::Ptr{Ptr{LDAPControl}})
    ccall((:ldap_controls_free, libldap), Cvoid, (Ptr{Ptr{LDAPControl}},), ctrls_handle)
end

#=
LDAP_F( int )
ldap_parse_result LDAP_P((
    LDAP            *ld,
    LDAPMessage     *res,
    int             *errcodep,
    char            **matcheddnp,
    char            **errmsgp,
    char            ***referralsp,
    LDAPControl     ***serverctrls,
    int             freeit ));
=#
function ldap_parse_result(ldp_handle::Ptr{Cvoid}, result_handle::Ptr{Cvoid},
        err_ref::Ref{Cint}, match_str_ref::Ref{Ptr{UInt8}}, err_msg_ref::Ref{Ptr{UInt8}},
        referralsp_ref::Ref{Ptr{Ptr{UInt8}}}, controls_ref::Ref{Ptr{Ptr{LDAPControl}}}, freeit::Bool)

    ccall((:ldap_parse_result, libldap), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{Cint}, Ref{Ptr{UInt8}}, Ref{Ptr{UInt8}}, Ref{Ptr{Ptr{UInt8}}}, Ref{Ptr{Ptr{LDAPControl}}}, Cint),
        ldp_handle, result_handle, err_ref, match_str_ref, err_msg_ref, referralsp_ref, controls_ref, freeit)
end

#LDAP_F( LDAPMessage * )
#ldap_first_reference LDAP_P((
#    LDAP *ld,
#    LDAPMessage *chain ));
function ldap_first_reference(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_first_reference, libldap), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

#LDAP_F( LDAPMessage * )
#ldap_next_reference LDAP_P((
#    LDAP *ld,
#    LDAPMessage *ref ));
function ldap_next_reference(ldp_handle::Ptr{Cvoid}, msg_handle::Ptr{Cvoid})
    ccall((:ldap_next_reference, libldap), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, msg_handle)
end

#=
LDAP_F( void )
ldap_memvfree LDAP_P((
    void** v ));
=#
function ldap_memvfree(ptr::Ptr{Ptr{T}}) where {T}
    ccall((:ldap_memvfree, libldap), Cvoid, (Ptr{Ptr{T}},), ptr)
end

#LDAP_F( char * )
#ldap_first_attribute LDAP_P((
#  LDAP *ld,
#  LDAPMessage *entry,
#  BerElement **ber ));
function ldap_first_attribute(ldp_handle::Ptr{Cvoid}, entry_handle::Ptr{Cvoid}, ber_element_handle_ref::Ref{Ptr{Cvoid}})
    ccall((:ldap_first_attribute, libldap), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{Ptr{Cvoid}}), ldp_handle, entry_handle, ber_element_handle_ref)
end

#=
LDAP_F( char * )
ldap_next_attribute LDAP_P((
    LDAP *ld,
    LDAPMessage *entry,
    BerElement *ber ));
=#
function ldap_next_attribute(ldp_handle::Ptr{Cvoid}, entry_handle::Ptr{Cvoid}, ber_element_handle::Ptr{Cvoid})
    ccall((:ldap_next_attribute, libldap), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), ldp_handle, entry_handle, ber_element_handle)
end

#LDAP_F( char ** )
#ldap_get_values LDAP_P((
#    LDAP *ld,
#    LDAPMessage *entry,
#    LDAP_CONST char *target ));
function ldap_get_values(ldp_handle::Ptr{Cvoid}, entry_handle::Ptr{Cvoid}, attribute::AbstractString)
    ccall((:ldap_get_values, libldap), Ptr{Ptr{UInt8}}, (Ptr{Cvoid}, Ptr{Cvoid}, Cstring), ldp_handle, entry_handle, attribute)
end

#LDAP_F( void )
#ldap_value_free LDAP_P((
#    char **vals ));
function ldap_value_free(vals::Ptr{Ptr{UInt8}})
    ccall((:ldap_value_free, libldap), Cvoid, (Ptr{Ptr{UInt8}},), vals)
end


text:
# module Jive

"""
    Skipped

Skipped expressions are in `Skipped.expressions`
"""
module Skipped
expressions = []
end # Jive.Skipped

isjiveskip() = get(ENV, "JIVE_SKIP", "1") == "1"

"""
    @skip

skip the expression.
"""
macro skip(expr::Expr)
    if isjiveskip()
        typ = expr.head
        if typ in (:module, :struct)
            name = expr.args[2]
            push!(Skipped.expressions, typ=>name)
        elseif typ in (:function, :macro)
            fexpr = expr.args[1]
            push!(Skipped.expressions, typ=>first(fexpr.args))
        elseif typ === :call
            push!(Skipped.expressions, typ=>first(expr.args))
        else
            # @info :typ (typ, expr.args)
            # do if quote block using return (=) vect
            push!(Skipped.expressions, typ)
        end
        nothing
    else
        esc(expr)
    end
end

macro skip(node::QuoteNode)
    if isjiveskip()
        push!(Skipped.expressions, node.value)
        nothing
    else
        esc(node)
    end
end

macro skip(sym::Symbol)
    if isjiveskip()
        if sym === :nothing
            push!(Skipped.expressions, nothing)
        else
            push!(Skipped.expressions, sym)
        end
        nothing
    else
        esc(sym)
    end
end

macro skip(val::Any)
    if isjiveskip()
        push!(Skipped.expressions, val)
        nothing
    else
        esc(val)
    end
end

# module Jive


text:
module TriangleIntersect

import Base: convert, intersect
using LinearAlgebra: cross, dot, norm, normalize
using GeometryBasics: AbstractPoint, Ngon, Point, Vec, Triangle, decompose

export IntersectableTriangle, Ray, Intersection, no_intersection

const AbstractTriangle{D} = Ngon{D, <:Number, 3, <:AbstractPoint{D}}

# 3 dimensions
const D = 3;

struct IntersectableTriangle
    a::Point{D}
    v1::Point{D}
    v2::Point{D}
    normal::Vec{D}
    v1v1::Float64
    v2v2::Float64
    v1v2::Float64
    denom::Float64
    function IntersectableTriangle(a::Point{D}, b::Point{D}, c::Point{D})
        v1 = b - a
        v2 = c - a
        normal = cross(v1, v2)
        v1v1 = dot(v1, v1)
        v2v2 = dot(v2, v2)
        v1v2 = dot(v1, v2)
        denom = v1v2 * v1v2 - v1v1 * v2v2
        new(a, v1, v2, normal, v1v1, v2v2, v1v2, denom)
    end
end

convert(::Type{IntersectableTriangle}, t::AbstractTriangle{D}) = IntersectableTriangle(decompose(Point, t)...)

struct Ray
    origin::Point{D}
    direction::Vec{D}
    Ray(a, b) = new(a, normalize(b))
end

struct Intersection
    ip::Point{D} # intersecting point
    id::Float64 # intersecting distance
    is_intersection::Bool
end

const no_intersection = Intersection(Point(0, 0, 0), 0.0, false)

function intersect(r::Ray, t::IntersectableTriangle)
    denom = dot(t.normal, r.direction)
    denom == 0 && return no_intersection
    ri = dot(t.normal, (t.a - r.origin)) / denom
    ri <= 0 && return no_intersection
    plane_intersection =  ri * r.direction + r.origin
    w = plane_intersection - t.a
    wv1 = dot(w, t.v1)
    wv2 = dot(w, t.v2)
    s_intersection = (t.v1v2 * wv2 - t.v2v2 * wv1) / t.denom
    s_intersection <= 0 && return no_intersection
    s_intersection >= 1 && return no_intersection
    t_intersection = (t.v1v2 * wv1 - t.v1v1 * wv2) / t.denom
    t_intersection <= 0 && return no_intersection
    t_intersection >= 1 && return no_intersection
    s_intersection + t_intersection >= 1 && return no_intersection
    Intersection(t.a + s_intersection * t.v1 + t_intersection * t.v2, ri, true)
end

intersect(r::Ray, t::AbstractTriangle{D}) = intersect(r, convert(IntersectableTriangle, t))

end # module


text:
using Revise
using Layered
using Animations
using Colors

function test()

    duration = 4

    a_points = Animation(
        0, reshape(P.(rand(25), rand(25)), (5, 5)),
        sineio(n=2, yoyo=true, prewait=0.3),
        duration, P.(grid(LinRange(0, 1, 5), LinRange(0, 1, 5))...))

    function frame(t)
        c, tl = canvas(2, 2)

        l = layer_in_rect!(tl, c.rect, :w, :norm, margin=30)

        circles!(l) do
            Circle.(a_points(t), 0.05)
        end + Fill(:frac => f -> LCHuv(70, 70, f * 360)) + Stroke(nothing)

        c

    end

    filepath = joinpath(dirname(pathof(Animations)), "../misc/example_array.gif")
    record(frame, filepath, 30, duration)

end; test()


text:
# In our fist example, we will build a basic left-right
# visualization using the standard layout options,
# to see how a 'test set' gets transformed by the
# anayltic continuation of `Arctan`.
using ComplexVisual
@ComplexVisual.import_huge

# Test if user has `ComplexVisualGtk` installed. It is
# possible to use `ComplexVisual` as a standalone package,
# though the interactive mode is not possible
const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

# If the goal function has a discontinuity, it is
# recommended to add a `cut_test`. This helps in visualization
cut_test1 = cv_create_angle_cross_test(+π/2, 1.0, Inf; δ=1e-2)
cut_test2 = cv_create_angle_cross_test(-π/2, 1.0, Inf; δ=1e-2)

cut_test = (w,z) -> cut_test1(w, z) || cut_test2(w, z)

# Our transformation is the complex definition of Arctan
trafo = z -> log((1+z*1im)/(1-z*1im))/(2im)

# We now set up our layout
domain   = CV_Math2DCanvas(-2.0 + 2.0im, 2.0 - 2.0im, 100)
codomain = CV_Math2DCanvas(-2.0 + 2.0im, 2.0 - 2.0im, 100)

# The `lr_painters_kwargs` gives options to change the shape,
# style, colors, etc. of the interactive parts.
# Similarly, `lr_start_kwargs` has options to change how the
# interactive parts start. See `lrdomains.jl` for details.
# Change `:z_start` if you want to save different pictures
lr_start_kwargs = Dict(:z_start => 0.2 +0.8im)
scene = cv_scene_lr_std(trafo, domain, codomain; cut_test, lr_start_kwargs)

# `scene` is a CV_SceneSetupChain, which can be used to
# interact.
# `cv_get_can_layout` is an internal getter-function.
# It serves to get `scene.layout`, which is all that
# we need to create a `png` image.
cv_save_image(cv_get_can_layout(scene), "./Arctan.png")

# will not do anything unless ComplexVisualGtk is installed
if CVG !== nothing
    handler = CVG.cvg_visualize(scene)
    CVG.cvg_wait_for_destroy(handler.window)

    CVG.cvg_close(handler); CVG.cv_destroy(scene)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
using Images
using StaticArrays
using VideoIO
using SLAM

function main()
    save_dir = "/home/pxl-th/projects/slam-uni"
    vpath = "/home/pxl-th/Videos/MOV_0014.MOV"

    focal = 910.0
    width, height = 1920.0, 1080.0
    cx, cy = width / 2.0, height / 2.0
    δt = 1.0 / 30.0
    timestamp = 0.0

    camera = Camera(;fx=focal, fy=focal, cx, cy, height, width)
    params = Params(;stereo=false, do_local_bundle_adjustment=true, map_filtering=false)
    saver = ReplaySaver()
    slam_manager = SlamManager(params, camera; slam_io=saver)
    slam_manager_thread = Threads.@spawn run!(slam_manager)

    t1 = time()

    for (i, frame) in enumerate(VideoIO.openvideo(vpath))
        i == 150 && break

        frame = Gray{Float64}.(frame)
        add_image!(slam_manager, frame, timestamp)
        timestamp += δt

        q_size = get_queue_size(slam_manager)
        f_size = length(slam_manager.mapper.estimator.frame_queue)
        m_size = length(slam_manager.mapper.keyframe_queue)
        while q_size > 0 || f_size > 0 || m_size > 0
            sleep(1e-2)
            q_size = get_queue_size(slam_manager)
            f_size = length(slam_manager.mapper.estimator.frame_queue)
            m_size = length(slam_manager.mapper.keyframe_queue)
        end

        sleep(1e-2)
    end

    slam_manager.exit_required = true
    wait(slam_manager_thread)

    t2 = time()
    @info "SLAM took: $(t2 - t1) seconds."

    SLAM.save(saver, save_dir)
    @info "Saved ReplaySaver."
end


text:
export solve, Result, nresults, nsolutions, nfinite, nsingular, natinfinity,
    nfailed, nnonsingular, nreal, ntracked, finite, results, mapresults,
    failed, atinfinity, singular, nonsingular, seed,
    solutions, realsolutions, multiplicities, uniquesolutions, statistics

"""
    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 rather than projective space. Currently this is only supported for parameter homotopies.
* `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...)
    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!(results, tracker, path_result_details; finite_results_only=!save_all_paths)
    end
    Result(results, ntracked, seed(tracker))
end

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, 0.1, "Tracking $n paths... ")
    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)
            k = 0
            for batch in batches
                prepare_batch!(batch_tracker, batch)
                ccall(:jl_threading_run, Ref{Cvoid}, (Any,), batch_tracker)

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

                update_progress!(progress, k, 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
            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 issingular(R)
        stats.singular_real += isreal(R)
        stats.singular += 1
    else
        stats.regular_real += isreal(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)
    if finished
        ProgressMeter.finish!(progress; showvalues=showvalues)
    end

    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
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)
end

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

    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)
     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)
    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], 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 == 0 && cond < 1e12
                    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

"""
    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.
"""
struct Result{V}
    pathresults::Vector{PathResult{V}}
    tracked_paths::Int
    seed::Int
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}

const Results = Union{Result, Vector{<:PathResult}}
const ProjectiveResult = Result{<:PVector}
"""
    nresults(result; onlyreal=false, realtol=1e-6, onlynonsingular=false, singulartol=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, onlyreal=true, realtol=1e-8, onlynonsingular=true)
```
"""
function nresults(R::Results; onlyreal=false, realtol=1e-6,
    onlynonsingular=false, onlysingular=false, singulartol=1e10, onlyfinite=true)
    count(R) do r
        (!onlyreal || isreal(r, realtol)) &&
        (!onlynonsingular || isnonsingular(r, singulartol)) &&
        (!onlysingular || issingular(r, singulartol)) &&
        (!onlyfinite || isfinite(r) || isprojective(r))
    end
end

"""
    statistics(R::Result; onlyreal=false, realtol=1e-6,
        onlynonsingular=false, onlysingular=false, singulartol=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, atinfinity = 0, failed = 0, total = 2)
"""
function statistics(R::Results, onlyreal=false, realtol=1e-6,
    onlynonsingular=false, onlysingular=false, singulartol=1e10)

    failed = atinfinity = nonsingular = singular = real_nonsingular = real_singular = 0

    for r in R
        if isfailed(r)
            failed += 1
        elseif issingular(r, singulartol)
            if isreal(r, realtol)
                real_singular += 1
            end
            singular += 1
        elseif !isprojective(r) && !isfinite(r)
            atinfinity += 1
        else # finite, nonsingular
            if isreal(r, realtol)
                real_nonsingular += 1
            end
            nonsingular += 1
        end
    end
    (nonsingular = nonsingular,
    singular = singular,
    real_nonsingular = real_nonsingular,
    real_singular = real_singular,
    real = real_nonsingular + real_singular,
    atinfinity = atinfinity,
    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) = count(issuccess, R)

"""
    nsingular(result; tol=1e10)

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`.
"""
nsingular(R::Results; tol=1e10) = count(r -> issingular(r, tol), R)

"""
    natinfinity(result)

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

"""
    nafailed(result)

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

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

The number of non-singular solutions.
"""
nnonsingular(R::Result; tol = 1e10) = count(r -> isnonsingular(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 -> isreal(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; onlyreal=false, realtol=1e-6, onlynonsingular=false, onlysigular=false, singulartol=1e10, onlyfinite=true)

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).
realsolutions = results(R, onlyreal=true, onlynonsingular=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).
realsolutions = mapresults(solution, R, onlyreal=true)
```
"""
function mapresults(f::Function, R::Results;
    onlyreal=false, realtol=1e-6, onlynonsingular=false, onlysingular=false, singulartol=1e10,
    onlyfinite=true)
    [f(r) for r in R if
        (!onlyreal || isreal(r, realtol)) &&
        (!onlynonsingular || isnonsingular(r, singulartol)) &&
        (!onlysingular || issingular(r, singulartol)) &&
        (!onlyfinite || isfinite(r) || isprojective(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

"""
    realsolutions(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 `onlyreal` is always `true`
and `realtol` is now `tol`.

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

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

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

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

Return all `PathResult`s for which the solution is singular. This is just a shorthand
for `results(R; onlysingular=true, conditions...)`. For the possible `conditions` see [`results`](@ref).
"""
function singular(R::Results; singulartol=1e10, tol=singulartol, kwargs...)
    results(R; onlysingular=true, singulartol=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 [`isreal`](@ref) for details regarding the determination of 'realness'.
"""
Base.real(R::Results; tol=1e-6) = [r for r in R if isreal(r, tol)]

"""
    failed(result)

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

"""
    atinfinity(result::AffineResult)

Get all results where the solutions is at infinity.
"""
atinfinity(R::Results) = [r for r in R if isatinfinity(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

"""
    uniquesolutions(R::Result; tol=1e-6, multiplicities=false, conditions...)

Return all *unique* solutions. If `multiplicities` is `true`, then
all *unique* solutions with their correspnding multiplicities as pairs `(s, m)`
where `s` is the solution and `m` the multiplicity are returned.
For the possible `conditions` see [`results`](@ref).

## Example
```julia-repl
julia> @polyvar x;
julia> uniquesolutions([(x-3)^3*(x+2)], multiplicities=true)
[([3.0+0.0im], 3), ([-2.0+0.0im], 1)]
julia> uniquesolutions([(x-3)^3*(x+2)])
[[3.0+0.0im], [-2.0+0.0im]]
```
"""
function uniquesolutions(R::Results; tol=1e-6, multiplicities=false, conditions...)
    uniquesolutions(R, Val(multiplicities); tol=tol, conditions...)
end

function uniquesolutions(R::Results, ::Val{Multiplicities}; tol=1e-6, conditions...)  where {Multiplicities}
    sols = solutions(R; conditions...)
    M = multiplicities(sols; tol=tol)
    indicator = trues(length(sols))
    uniques = map(M) do m
        for k in m
            indicator[k] = false
        end
        if Multiplicities
            (sols[m[1]], length(m))
        else
            sols[m[1]]
        end
    end
    for (k, s) in enumerate(sols)
        if indicator[k]
            if Multiplicities
                push!(uniques, (s, 1))
            else
                push!(uniques, s)
            end
        end
    end
    uniques
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.atinfinity > 0 &&
        println(io, "• $(s.atinfinity) $(plural("solution", s.atinfinity)) 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 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 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.atinfinity > 0
        print(io, "$(s.atinfinity) atinfinity")
    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 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, onlynonsingular=true)
    elseif i == 3 && s.singular > 0
        return finite(r, onlysingular=true)
    elseif i == 4 && (s.real_nonsingular+s.real_singular) > 0
        return finite(r, onlyreal = true)
    elseif i == 5 && s.atinfinity > 0
        return atinfinity(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 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, onlynonsingular=true)
    elseif i == 3 && s.singular > 0
        return finite(r, onlysingular=true)
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        return finite(r, onlyreal=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))
    multiplicities_dict = Dict{Int, Vector{Vector{Int}}}()
    for m in multiplicities(solution, singular(result))
    	if haskey(multiplicities_dict, length(m))
    		push!(multiplicities_dict[length(m)], m)
    	else
    		multiplicities_dict[length(m)] = [m]
    	end
    end
    if isempty(multiplicities_dict)
        n_mults_total = 0
    else
        n_mults_total = sum(M -> sum(length, M), values(multiplicities_dict))
    end

    curr_n_real_sols = 0
    curr_n_sols = 0

    off = n_mults_total < stats.singular
    data = Matrix{Int}(undef, off+length(multiplicities_dict), 4)
    for (i, k) in enumerate(sort(collect(keys(multiplicities_dict))))
    	data[i+off, 1] = k
    	n_real_solsᵢ = count(M -> isreal(result.pathresults[first(M)]), multiplicities_dict[k])
    	n_solsᵢ = length(multiplicities_dict[k])
    	data[i+off, 2] = n_real_solsᵢ
    	data[i+off, 3] = n_solsᵢ - n_real_solsᵢ
    	data[i+off, 4] = n_solsᵢ

    	curr_n_real_sols += k * n_real_solsᵢ
    	curr_n_sols += k * n_solsᵢ
    end

    if off
    	data[1,1] = 1
    	data[1,2] = stats.real_singular - curr_n_real_sols
    	data[1,3] = (stats.singular - curr_n_sols) - (stats.real_singular - curr_n_real_sols)
    	data[1,4] = stats.singular - curr_n_sols
    end

    headers = ["multiplicity", "# real", "# non-real", "# total"]
    PrettyTables.pretty_table(io, data, headers;
            alignment=:c,
            header_crayon=PrettyTables.Crayon(bold=false))
end


text:
using LinearAlgebra

abstract type DescentDirectionMethod end

basis(i, n) = [k == i ? 1.0 : 0.0 for k in 1:n]

diff_complex(f, x; h = 1e-20) = imag(f(x .+ h * im)) / h

function backtracking_line_search(f, ∇f, x, d, α; p=0.5, β=1e-4)
    y, g = f(x), ∇f(x)
    while f(x + α * d) > y + β * α * (g ⋅ d)
        α *= p
    end
    return α
end


text:
# MDP parameters, ValueIteration minimizes the cost => cost is positive, reward is negative
no_states = 10
horizon = 5
actions = [:l, :r]
actionCost = 1.
actionsImpact = Base.ImmutableDict(:l => -1, :r => 1)
reward_states = [1, no_states]
reward = -10.
discount_factor = 1.
noise = .6

# MDPs initialization
mdp = FHExample(no_states, horizon, actions, actionCost, actionsImpact, reward_states, reward, discount_factor, noise)

# initialize the solver
# max_iterations: maximum number of iterations value iteration runs for (default is 100)
# belres: the value of Bellman residual used in the solver
solver = ValueIterationSolver(max_iterations=10, belres=1e-3, include_Q=true);

# Solve Value Iteration
VIPolicy = DiscreteValueIteration.solve(solver, mdp);

# Solve Finite Horizon by Value Iteration
FHPolicy = FiniteHorizonPOMDPs.solve(mdp);

# Compare resulting policies
@test all((FiniteHorizonPOMDPs.action(FHPolicy, s) == action(VIPolicy, s) for s in states(mdp)))

# # Evaluates to false because of different elements Types (ExampleState vs FHExampleState)
fh_states = collect(Iterators.flatten([FiniteHorizonPOMDPs.stage_states(mdp, i) for i=1:mdp.horizon]))
ih_states = states(mdp)

z = zip(fh_states, ih_states)

@test all((fh == ih for (fh, ih) in z))


text:
# # Cyclic symmetry for Sums of Hermitian Squares

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Symmetry/cyclic.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Symmetry/cyclic.ipynb)

using Test #src

# We start by defining the Cyclic group.

using GroupsCore
import PermutationGroups

struct CyclicElem <: GroupElement
    n::Int
    id::Int
end
Base.:(==)(a::CyclicElem, b::CyclicElem) = a.n == b.n && a.id == b.id
Base.inv(el::CyclicElem) = CyclicElem(el.n, (el.n - el.id) % el.n)

function Base.:*(a::CyclicElem, b::CyclicElem)
    return CyclicElem(a.n, (a.id + b.id) % a.n)
end
Base.:^(el::CyclicElem, k::Integer) = CyclicElem(el.n, (el.id * k) % el.n)

Base.conj(a::CyclicElem, b::CyclicElem) = inv(b) * a * b
Base.:^(a::CyclicElem, b::CyclicElem) = conj(a, b)

function PermutationGroups.order(el::CyclicElem)
    return div(el.n, gcd(el.n, el.id))
end

struct CyclicGroup <: Group
    n::Int
end
Base.eltype(::CyclicGroup) = CyclicElem
Base.one(c::Union{CyclicGroup, CyclicElem}) = CyclicElem(c.n, 0)
PermutationGroups.gens(c::CyclicGroup) = [CyclicElem(c.n, 1)]
PermutationGroups.order(::Type{T}, c::CyclicGroup) where {T} = convert(T, c.n)
function Base.iterate(c::CyclicGroup, prev::CyclicElem=CyclicElem(c.n, -1))
    id = prev.id + 1
    if id >= c.n
        return nothing
    else
        next = CyclicElem(c.n, id)
        return next, next
    end
end

# Now we define that the cyclic group acts on monomial by permuting variables
# cyclically. So for instance, `CyclicElem(3, 1)` would transform
# `x_1^3*x_2*x_3^4` into `x_1^4*x_2^3*x_3`.

import MultivariatePolynomials
const MP = MultivariatePolynomials

using SumOfSquares

struct Action{V<:MP.AbstractVariable} <: Symmetry.OnMonomials
    variables::Vector{V}
end
Symmetry.SymbolicWedderburn.coeff_type(::Action) = Float64
function Symmetry.SymbolicWedderburn.action(a::Action, el::CyclicElem, mono::MP.AbstractMonomial)
    return prod(MP.powers(mono), init=MP.constantmonomial(mono)) do (var, exp)
        index = findfirst(isequal(var), a.variables)
        new_index = mod1(index + el.id, el.n)
        return a.variables[new_index]^exp
    end
end

using DynamicPolynomials
@polyvar x[1:3]
action = Action(x)
g = CyclicElem(3, 1)
@test Symmetry.SymbolicWedderburn.action(action, g, x[1]^3 * x[2] * x[3]^4) == x[1]^4 * x[2]^3 * x[3] #src
Symmetry.SymbolicWedderburn.action(action, g, x[1]^3 * x[2] * x[3]^4)

# The following polynomial `poly` is invariant under the action of the group `G`.

N = 3
G = CyclicGroup(N)
poly = sum(x[i] * x[mod1(i + 1, N)] for i in 1:N) + sum(x.^2)
@test Symmetry.SymbolicWedderburn.action(action, g, poly) == poly #src
Symmetry.SymbolicWedderburn.action(action, g, poly)

# Let's now find the minimum of `p` by exploiting this symmetry.

import CSDP
solver = CSDP.Optimizer
model = Model(solver)
@variable(model, t)
@objective(model, Max, t)
pattern = Symmetry.Pattern(G, action)
con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
optimize!(model)
@test termination_status(model) == MOI.OPTIMAL #src
@test objective_value(model) ≈ 0.0 atol=1e-4 #src
solution_summary(model)

# Let's look at the symmetry adapted basis used.

gram = gram_matrix(con_ref).sub_gram_matrices #src
@test length(gram) == 3                       #src
@test gram[1].Q ≈ [0 0; 0 2] #src
@test length(gram[1].basis.polynomials) == 2 #src
@test gram[1].basis.polynomials[1] == 1 #src
@test gram[1].basis.polynomials[2] ≈ -sum(x)/√3 #src
@test gram[2].Q ≈ [0.5;;] #src
@test length(gram[2].basis.polynomials) == 1 #src
@test gram[2].basis.polynomials[1] ≈ (-2x[1] + x[2] + x[3])/√6 #src
@test gram[3].Q == gram[2].Q #src
@test length(gram[3].basis.polynomials) == 1 #src
@test gram[3].basis.polynomials[1] ≈ (-x[2] + x[3])/√2 #src
for gram in gram_matrix(con_ref).sub_gram_matrices
    println(gram.basis.polynomials)
    display(gram.Q)
end

# Let's look into more details at the last two elements of the basis.

basis = [(-2x[1] + x[2] + x[3])/√6, (-x[2] + x[3])/√2]

# This actually constitutes the basis for an invariant subspace corresponding
# to a group character of degree 2 and multiplicity 1.
# This means that it decomposes the semidefinite matrix into 2 blocks of size
# 1-by-1 that are equal. Indeed, we see above that `gram.Q` is identically equal for both.
# As the group is generated by one element `g`, we can just verify it by verifying
# its invariance under `g`.
# The image of each element under the basis is:

image = [Symmetry.SymbolicWedderburn.action(action, g, p) for p in basis]

# We can see that they are both still in the same 2-dimensional subspace.

a = -1/2
b = -√3/2
@test all(image .≈ [a -b; b a] * basis) #src
[a -b; b a] * basis

# In fact, these last two basis comes from the real decomposition of a complex one.

import CSDP
solver = CSDP.Optimizer
model = Model(solver)
add_bridge(model, SumOfSquares.COI.Bridges.Variable.HermitianToSymmetricPSDBridge)
add_bridge(model, SumOfSquares.COI.Bridges.Constraint.SplitZeroBridge)
MOI.Bridges.add_bridge(backend(model).optimizer, PolyJuMP.ZeroPolynomialBridge{Complex{Float64}})
MOI.Bridges.add_bridge(backend(model).optimizer, SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Complex{Float64}})
@variable(model, t)
@objective(model, Max, t)
pattern = Symmetry.Pattern(G, action)
cone = SumOfSquares.NonnegPolyInnerCone{SumOfSquares.COI.HermitianPositiveSemidefiniteConeTriangle}()
pp = (1.0 + 0.0im) * poly - (1.0 + 0.0im) * t
con_ref = @constraint(model, pp in cone, symmetry = pattern)
optimize!(model)
@test termination_status(model) == MOI.OPTIMAL #src
@test objective_value(model) ≈ 0.0 atol=1e-4 #src
solution_summary(model)

gram = gram_matrix(con_ref).sub_gram_matrices #src
@test length(gram) == 3                       #src
@test gram[1].Q ≈ [0 0; 0 2] #src
@test length(gram[1].basis.polynomials) == 2 #src
@test gram[1].basis.polynomials[1] == 1 #src
@test gram[1].basis.polynomials[2] ≈ -sum(x)/√3 #src
@test gram[2].Q ≈ [0.5;;] #src
@test length(gram[2].basis.polynomials) == 1 #src
@test gram[2].basis.polynomials[1] ≈ (basis[1] + basis[2] * im) / √2  #src
@test gram[3].Q ≈ [0.5;;] #src
@test length(gram[3].basis.polynomials) == 1 #src
@test gram[3].basis.polynomials[1] ≈ (basis[1] - basis[2] * im) / √2 #src
for gram in gram_matrix(con_ref).sub_gram_matrices
    println(gram.basis.polynomials)
    display(gram.Q)
end

# We can see that the real invariant subspace was in fact coming from two complex conjugate complex invariant subspaces:

complex_basis = basis[1] + im * basis[2]
image = Symmetry.SymbolicWedderburn.action(action, g, complex_basis)

# And there is a direct correspondance between the representation of the real and complex versions:

@test all(image .≈ (a + b * im) * complex_basis) #src
(a + b * im) * complex_basis

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
## ------------------------------------------------------
function plot_time_series(simid, D, ϵs, cg)

   Dstr = string(round(D; sigdigits = 3))
   tslim = [Inf, -Inf]

   Xp = plot()
   Xlim = [0.0, -Inf]

   sgp = plot()
   sglim = [0.0, -Inf]

   max_ϵs_ = maximum(ϵs)
   ls0 = 2.5

   for (ϵi, ϵ) in enumerate(ϵs)

      status = get_status(simid, (;D, ϵ, cg))
      (status == UNDONE_SIM_STATUS) && continue

      @info("Doing", D, ϵ, cg, status)

      (ts, Xts) = collect_ts(:Xts, simid, (;D, ϵ, cg))
      (ts, sgts) = collect_ts(:sgts, simid, (;D, ϵ, cg))

      ts = ts ./ maximum(ts) .* 100.0
      tslim = (-5.0, 105.0)

      for (p, panel, ys, ylim, ylabel) in [
         [Xp, "A", Xts, Xlim, _textbf("X")],
         [sgp, "B", sgts, sglim, _textbf("s_g")]
      ]
         ylim .= lims(ys, ylim)
         plot!(p, ts, ys;
            title = _textbf(panel),
            xlim = tslim, ylim,
            label = "",
            ylabel,
            color = :black,
            xlabel = _textbf("simulation progress (\\%)"),
            alpha = 0.7,
            lw = ls0 + (2.0 * ls0) * (ϵ / max_ϵs_)
         )
      end
   end

   # annotation
   # fontsize = 16
   # text = _textbf("D=", Dstr)
   # annotate!(Xp, [(70.0, last(Xlim)* 0.4, (text, fontsize, :left, :top, :black))])
   # annotate!(sgp, [(70.0, last(sglim)* 0.4, (text, fontsize, :left, :top, :black))])

   # ffile = sfig(ChE, [Xp, sgp],
   #    @fileid, "dyn_time_series", (;D, cg), ".png";
   #    layout = (1, 2)
   # )
   # @info("Done", D, ffile)

   return [Xp, sgp]
end