text:
const OCR_DATA_DIR = abspath(joinpath(pathof(@__MODULE__), "../../tessdata"))

@enum PageSegMode begin
    PSM_OSD_ONLY               = 0
    PSM_AUTO_OSD               = 1
    PSM_AUTO_ONLY              = 2
    PSM_AUTO                   = 3
    PSM_SINGLE_COLUMN          = 4
    PSM_SINGLE_BLOCK_VERT_TEXT = 5
    PSM_SINGLE_BLOCK           = 6
    PSM_SINGLE_LINE            = 7
    PSM_SINGLE_WORD            = 8
    PSM_CIRCLE_WORD            = 9
    PSM_SINGLE_CHAR            = 10
    PSM_SPARSE_TEXT            = 11
    PSM_SPARSE_TEXT_OSD        = 12
    PSM_RAW_LINE               = 13
end

download_ocr_language(lang, target = OCR_DATA_DIR) =
    Tesseract.download_languages(lang, target = target)

default_ocr_instance(ctx::VsContext) = ctx.ocr_instances[ctx.config.ocr_language]

function get_pix(img::AbstractMatrix)
    f = IOBuffer()
    s = Stream{format"PNG"}(f)
    save(s, img)
    buffer = take!(f)
    pix = pix_read(buffer)
end

create_ocr_instance(
    language::AbstractString;
    tess_datapath = OCR_DATA_DIR
) = TessInst(language, tess_datapath)

function _set_page_seg_mode(inst::TessInst, mode::PageSegMode)
    ccall(
        (:TessBaseAPISetPageSegMode, Tesseract.TESSERACT),
        Cvoid,
        (Ptr{Cvoid}, PageSegMode),
        inst, mode
    )
end

function _ocr(img::AbstractMatrix, instance::TessInst; resolution = 72, page_seg_mode = PSM_SINGLE_LINE)
    pix = get_pix(img)
    tess_image(instance, pix)
    tess_resolution(instance, resolution)
    _set_page_seg_mode(instance, page_seg_mode)
end

"""
    ocr([Type=String], image, language/instance/context; resolution = 72, strip = true)

Use Tesseract to OCR an image.
"""
function ocr(
    ::Type{String},
    img::AbstractMatrix,
    instance::TessInst;
    resolution = 72,
    strip = true,
    page_seg_mode::PageSegMode = PSM_SINGLE_LINE
)
    length(img) == 0 && return ""
    _ocr(img, instance; resolution = resolution, page_seg_mode = page_seg_mode)
    text = tess_text(instance)
    if strip
        text = Base.strip(text)
    end
    text
end

function ocr(
    ::Type{String},
    img::AbstractMatrix, language::AbstractString;
    tess_datapath = OCR_DATA_DIR,
    kwargs...
)
    instance = create_ocr_instance(language, tess_datapath = tess_datapath)
    ocr(String, img, instance; kwargs...)
end

ocr(
    ::Type{String},
    img::AbstractMatrix,
    ctx::VsContext;
    language = ctx.config.ocr_language,
    kwargs...
) = let text = ocr(String, img, ctx.ocr_instances[language]; kwargs...)
    is_cjk(language) ? remove_spaces(text) : text
end

function ocr(::Type{Int}, img::AbstractMatrix, args...; default::Int = 0, kwargs...)
    text = ocr(String, img, args...; kwargs...)
    parse_int(text, default)
end

ocr(img::AbstractMatrix, args...; kwargs...) = ocr(String, img, args...; kwargs...)

function ocr_tsv(img::AbstractMatrix, instance; kwargs...)
    length(img) == 0 && return ""
    _ocr(img, instance; kwargs...)
    tess_parsed_tsv(instance)
end

is_cjk(lang) = split(lang, '+') ∩ ("chi_sim", "chi_tra", "jpn", "kor") |> !isempty
function init_multiple_ocr!(ctx::VsContext, languages::AbstractVector{<:AbstractString}, tess_datapath = OCR_DATA_DIR)
    for lang in languages
        lang in keys(ctx.ocr_instances) && continue
        ctx.ocr_instances[lang] = create_ocr_instance(lang, tess_datapath = tess_datapath)
    end
    ctx
end

function ocr_multiple_lang(img::AbstractMatrix, ctx::VsContext)
    _, best_lang, best_words = find_closest(
        ctx.ocr_instances,
        should_break = (x) -> x[1] < 20  # Set a threshold to break early
    ) do (lang, inst)
        tsvs = [t for t in ocr_tsv(img, inst) if t.level == 5]  # filter the words
        average_conf = if length(tsvs) == 0
            -1.0
        else
            mean(t.conf for t in tsvs)
        end
        100 - average_conf, lang, tsvs
    end
    if length(best_words) == 0
        "", best_lang
    else
        join((t.text for t in best_words), is_cjk(best_lang) ? "" : " "), best_lang
    end
end


text:
# Latent ODE model
#
# Based on
# https://arxiv.org/abs/1806.07366
# https://arxiv.org/abs/2003.10775

struct LatentODE <: LatentDE end

@doc raw"""
    apply_feature_extractor(encoder::Encoder{GOKU}, x::Vector{Array})
    apply_feature_extractor(encoder::Encoder{LatentODE}, x::Vector{Array})
Applies the feature extractor layer contained in the `encoder` to the batch of input data x, usually reducing its dimensionality (i.e. extracting features).
# Arguments
`encoder`: Structure of type `Encoder` containing all the encoders layers.\
`x`: Input data, consists of matrices corresponding to the time frames and having size (`input data dimension` x `batch size`)
"""
apply_feature_extractor(encoder::Encoder{LatentODE}, x) = encoder.feature_extractor.(x)

@doc raw"""
    apply_pattern_extractor(encoder::Encoder{LatentODE}, fe_out)
Passes extracted features through the pattern_extractor layer contained in the `encoder`, returning a matrix containing patterns for the initial conditions.
# Arguments
`fe_out`: Output of feature extractor layer
"""
function apply_pattern_extractor(encoder::Encoder{LatentODE}, fe_out)
    pe_z₀ = encoder.pattern_extractor

    # reverse sequence
    fe_out_rev = reverse(fe_out)

    # pass it through the recurrent layers
    pe_out = map(pe_z₀, fe_out_rev)[end]

    # reset hidden states
    Flux.reset!(pe_z₀)

    return pe_out
end

@doc raw"""
    apply_latent_in(encoder::Encoder{LatentODE}, pe_out)
    apply_latent_in(encoder::Encoder{LatentODE}, pe_out)
Applies the `encoder`'s `latent_in` layer to `pe_out`, obtaining representations of the mean and log-variance of the initial conditions (and parameters in the case of GOKU) to use for sampling.
# Arguments
`pe_out`: output of pattern extractor layer
"""
function apply_latent_in(encoder::Encoder{LatentODE}, pe_out)
    li_μ_z₀, li_logσ²_z₀ = encoder.latent_in

    z₀_μ = li_μ_z₀(pe_out)
    z₀_logσ² = li_logσ²_z₀(pe_out)

    return z₀_μ, z₀_logσ²
end

@doc raw"""
    apply_latent_out(decoder::Decoder{LatentODE}, z̃₀)
Applies the `decoder`'s `latent_out` layer to `z̃₀`, obtaining the inferred initial conditions for use in the differential equation layer.
# Arguments
`z̃₀`: Sampled abstract representations of the initial conditions.
"""
apply_latent_out(decoder::Decoder{LatentODE}, z̃₀) = decoder.latent_out(z̃₀)

@doc raw"""
    diffeq_layer(decoder::Decoder{LatentODE}, ẑ₀, t)
Uses decoder.diffeq's DE solver to extrapolate the latent states (from the initial states ẑ₀) to time t.
"""
function diffeq_layer(decoder::Decoder{LatentODE}, ẑ₀, t)
    dudt = decoder.diffeq.dudt
    solver = decoder.diffeq.solver
    neural_model = decoder.diffeq.neural_model
    augment_dim = decoder.diffeq.augment_dim
    kwargs = decoder.diffeq.kwargs
    # sensealg = decoder.diffeq.sensealg

    # nODE = neural_model(dudt, (t[1], t[end]), solver, sensealg = sensealg, saveat = t)
    nODE = neural_model(dudt, (t[1], t[end]), solver; saveat = t, kwargs...)
    nODE = augment_dim == 0 ? nODE : AugmentedNDELayer(nODE, augment_dim)
    ẑ = Array(nODE(ẑ₀))

    # Transform the resulting output (mainly used for Kuramoto-like systems)
    ẑ = transform_after_diffeq(ẑ, decoder.diffeq)
    ẑ = Flux.unstack(ẑ, 3)

    return ẑ
end

@doc raw"""
    apply_reconstructor(decoder::Decoder{GOKU}, ẑ)
    apply_reconstructor(decoder::Decoder{LatentODE}, ẑ)
Reconstruct the initial data from the extrapolated latent states `ẑ`.
# Arguments
`decoder`: Structure of type Decoder containing all of the decoder layers.
`ẑ`: Latent states, consists of matrices corresponding to the time frames and having size (`latent data dimension` x `batch size`)
"""
apply_reconstructor(decoder::Decoder{LatentODE}, ẑ) = decoder.reconstructor.(ẑ)

@doc raw"""
    sample(μ::T, logσ²::T, model::LatentDiffEqModel{LatentODE}) where T <: Array
Samples abstract representations of the initial state from the normal distribution with mean μ and variance exp(logσ²).
"""
function sample(μ::T, logσ²::T, model::LatentDiffEqModel{LatentODE}) where T <: Array
    z₀_μ = μ
    z₀_logσ² = logσ²

    ẑ₀ = z₀_μ + randn(Float32, size(z₀_logσ²)) .* exp.(z₀_logσ²/2f0)

    return ẑ₀
end

function sample(μ::T, logσ²::T, model::LatentDiffEqModel{LatentODE}) where T <: Flux.CUDA.CuArray
    z₀_μ = μ
    z₀_logσ² = logσ²

    ẑ₀ = z₀_μ + gpu(randn(Float32, size(z₀_logσ²))) .* exp.(z₀_logσ²/2f0)

    return ẑ₀
end

@doc raw"""
    default_layers(model_type, input_dim::Int, diffeq;
        device = cpu,
        hidden_dim_resnet = 200, rnn_input_dim = 32,
        rnn_output_dim = 16, latent_dim = 16,
        latent_to_diffeq_dim = 200, θ_activation = softplus,
        output_activation = σ, init = Flux.kaiming_uniform(gain = 1/sqrt(3)))
Generates default encoder and decoder layers that are to be fed into the LatentDiffEqModel.
# Arguments
`model_type`: GOKU() or LatentODE()\
`input_dim`: Dimension of input\
`diffeq`: Differential equations structure, containing fields `prob`, `solver` and `sensealg`, which correspond to DifferentialEquations.jl's problem, solver and sensitivity algorithm, respectively. For an example, see https://github.com/gabrevaya/LatentDiffEq.jl/blob/master/examples/pendulum_friction-less/pendulum.jl
"""
function default_layers(model_type::LatentODE, input_dim, diffeq; device = cpu,
                            hidden_dim_resnet = 200, rnn_input_dim = 32,
                            rnn_output_dim = 32, latent_to_diffeq_dim = 200,
                            output_activation = σ, init = Flux.kaiming_uniform(gain = 1/sqrt(3)))

    latent_dim_in = diffeq.latent_dim_in
    latent_dim_out = diffeq.latent_dim_out

    ######################
    ### Encoder layers ###
    ######################

    # Resnet
    l1 = Dense(input_dim, hidden_dim_resnet, relu, init = init)
    l2 = Dense(hidden_dim_resnet, hidden_dim_resnet, relu, init = init)
    l3 = Dense(hidden_dim_resnet, hidden_dim_resnet, relu, init = init)
    l4 = Dense(hidden_dim_resnet, rnn_input_dim, relu, init = init)
    feature_extractor = Chain(l1,
                                SkipConnection(l2, +),
                                SkipConnection(l3, +),
                                l4) |> device

    # RNN
    pattern_extractor = Chain(RNN(rnn_input_dim, rnn_output_dim, relu, init = init),
                                RNN(rnn_output_dim, rnn_output_dim, relu, init = init)) |> device

    # final fully connected layers before sampling
    li_μ_z₀ = Dense(rnn_output_dim, latent_dim_in, init = init) |> device
    li_logσ²_z₀ = Dense(rnn_output_dim, latent_dim_in, init = init) |> device

    latent_in = (li_μ_z₀, li_logσ²_z₀)

    encoder_layers = (feature_extractor, pattern_extractor, latent_in)

    ######################
    ### Decoder layers ###
    ######################

    # going back to the input dimensions
    # Resnet
    l1 = Dense(latent_dim_out, hidden_dim_resnet, relu, init = init)
    l2 = Dense(hidden_dim_resnet, hidden_dim_resnet, relu, init = init)
    l3 = Dense(hidden_dim_resnet, hidden_dim_resnet, relu, init = init)
    l4 = Dense(hidden_dim_resnet, input_dim, output_activation, init = init)
    reconstructor = Chain(l1,
                            SkipConnection(l2, +),
                            SkipConnection(l3, +),
                            l4)  |> device

    decoder_layers = (x -> x, diffeq, reconstructor)

    return encoder_layers, decoder_layers
end

ext: jl


text:
""" Comissão de vendas
# Dado o salário fixo,
# o valor das vendas efetuadas pelo vendedor de uma empresa e
# sabendo-se que ele recebe uma comissão de 3% sobre o total das vendas até R$ 1.500,00
# e 5% sobre o que ultrapassar este valor, calcular e escrever o seu salário total.
"""

println("informe o salario")
salary = parse(Float16, readline())

println("informe o valor das vendas")
commission = parse(Float16, readline())

if commission <= 1500
    salary += (commission*0.03)
else
    salary += (commission*0.05)
end
println(salary)

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
runall = false;
@testset "Testing traitLikDiscrete" begin
@testset "Testing Substitution Models, P and Q matrices" begin

m1 = BinaryTraitSubstitutionModel(1.0, 2.0);
@test_nowarn show(DevNull, m1)
m1 = BinaryTraitSubstitutionModel(1.0,2.0, ["carnivory", "non-carnivory"]);
@test nStates(m1)==2
@test PhyloNetworks.nparams(m1)==2
@test_nowarn show(DevNull, m1)
@test_throws ErrorException PhyloNetworks.BinaryTraitSubstitutionModel(-1.0,2.0)
m2 = EqualRatesSubstitutionModel(4, 3.0);
@test nStates(m2)==4
@test PhyloNetworks.nparams(m2)==1
m2 = EqualRatesSubstitutionModel(4, 3.0, ["S1","S2","S3","S4"]);
@test_nowarn show(DevNull, m2)
m3 = TwoBinaryTraitSubstitutionModel([2.0,1.2,1.1,2.2,1.0,3.1,2.0,1.1],
["carnivory", "noncarnivory", "wet", "dry"]);
@test_nowarn show(DevNull, m3)
@test nStates(m3)==4
@test PhyloNetworks.nparams(m3)==8

@test Q(m1) == SMatrix{2,2}(-1.0, 2.0, 1.0, -2.0)
@test Q(m2) == SMatrix{4,4}(-9.0, 3, 3, 3, 3, -9, 3, 3, 3, 3, -9, 3, 3, 3, 3, -9)
@test Q(m3) ≈ SMatrix{4,4}(-3.0, 3.1, 1.2, 0.0, 1.0, -4.2, 0.0, 2.2, 2.0, 0.0, -3.2, 1.1, 0.0, 1.1, 2.0, -3.3) atol=1e-4
@test P(m1, 0.5) ≈ SMatrix{2,2}(0.7410433867161432,0.5179132265677134,0.2589566132838567,0.4820867734322865) atol=1e-15
@test P(m2, 0.1) ≈ SMatrix{4,4}(0.4758956589341516,0.17470144702194956,0.17470144702194962,0.17470144702194945,0.17470144702194956,0.47589565893415153,0.17470144702194967,0.17470144702194945,0.17470144702194962,0.17470144702194967,0.4758956589341516,0.17470144702194945,0.17470144702194945,0.17470144702194945,0.17470144702194945,0.4758956589341518) atol=1e-15
@test P(m3, 0.5) ≈ SMatrix{4,4}(0.39839916380463375,0.36847565707393248,0.23055614536582461,0.22576141081414305,0.1545971371424259,0.25768553619230444,0.14816051303688715,0.24300762748855972,0.29194735222005136,0.20198250750421617,0.35349416558860958,0.20267178083856716,0.15505634683288913,0.17185629922954704,0.26778917600867863,0.32855918085873009) atol=1e-15
@test P(m1, [0.02,0.01]) ≈ StaticArrays.SArray{Tuple{2,2},Float64,2,4}[[0.980588 0.0194118; 0.0388236 0.961176], [0.990149 0.00985149; 0.019703 0.980297]] atol=1e-6
@test P(m2, [0.02,0.01]) ≈ Array{Float64,2}[[0.839971 0.053343 0.053343 0.053343; 0.053343 0.839971 0.053343 0.053343; 0.053343 0.053343 0.839971 0.053343; 0.053343 0.053343 0.053343 0.839971], [0.91519 0.0282699 0.0282699 0.0282699; 0.0282699 0.91519 0.0282699 0.0282699; 0.0282699 0.0282699 0.91519 0.0282699; 0.0282699 0.0282699 0.0282699 0.91519]] atol=1e-6

end

@testset "Testing random discrete trait simulation" begin

m1 = BinaryTraitSubstitutionModel(1.0,2.0, ["carnivory", "non-carnivory"]);
m2 = EqualRatesSubstitutionModel(4, 3.0, ["S1","S2","S3","S4"]);
# on a single branch
srand(12345);
@test randomTrait(m1, 0.2, [1,2,1,2,2]) == [1,2,1,1,2]
srand(12345);
@test randomTrait(m2, 0.05, [1,3,4,2,1]) == [1,3,4,2,1]
# on a network
net = readTopology("(A:1.0,(B:1.0,(C:1.0,D:1.0):1.0):1.0);")
srand(12345);
a,b = randomTrait(m1, net)
@test a == [1 2 1 1 1 1 2]
@test b == ["-2", "-3", "-4", "D", "C", "B", "A"]
if runall
    for e in net.edge e.length = 10.0; end
    @time a,b = randomTrait(m1, net; ntraits=100000) # ~ 0.014 seconds
    mean(a[:,1]) # expect 1.5 at root
    mean(a[:,2]) # expect 1.333 at other nodes
    @time a,b = randomTrait(m2, net; ntraits=100000) # ~ 0.02 seconds
    length([x for x in a[:,1] if x==4])/length(a[:,1]) # expect 0.25
    length([x for x in a[:,2] if x==4])/length(a[:,2])
    length([x for x in a[:,3] if x==4])/length(a[:,3])
    length([x for x in a[:,4] if x==4])/length(a[:,4])
    length([x for x in a[:,5] if x==4])/length(a[:,5])
    length([x for x in a[:,6] if x==4])/length(a[:,6])
    length([x for x in a[:,7] if x==4])/length(a[:,7]) # expect 0.25
end

net2 = readTopology("(((A:4.0,(B:1.0)#H1:1.1::0.9):0.5,(C:0.6,#H1:1.0):1.0):3.0,D:5.0);")
srand(12345);
a,b = randomTrait(m1, net2; keepInternal=false)
@test a == [1  1  1  2]
@test b == ["D", "C", "B", "A"]
srand(12345);
a,b = randomTrait(m1, net2; keepInternal=true)
@test a == [1  2  1  1  1  1  1  1  1]
@test b == ["-2", "D", "-3", "-6", "C", "-4", "#H1", "B", "A"]
if runall
    for e in net2.edge
        if e.hybrid
            e.length = 0.0
        end
    end
    a,b = randomTrait(m1, net2; ntraits=100000)
    # plot(net2, showNodeNumber=true) shows: H1 listed 7th, parents listed 4th and 6th
    c = map( != , a[:, 4],a[:, 6] ); # traits when parents have different traits
    n1 = sum(map( ==, a[c,7],a[c,6] )) # 39644 traits: hybrid ≠ major parent
    n2 = sum(map( ==, a[c,7],a[c,4] )) #  4401 traits: hybrid ≠ minor parent
    n1/sum(c) # expected 0.9
    n2/sum(c) # expected 0.1
    for e in net2.edge
        e.length = 0.0
    end
    net2.edge[4].length = 10.0
    a,b = randomTrait(m1, net2; ntraits=100000);
    a[:, 1] == a[:, 2]  # true: root = leaf D, as expected
    a[:, 1] == a[:, 5]  # true: root = leaf C
    mean(a[:, 6]) # expected 1.3333
    a[:, 6] == a[:, 9] # true: major hybrid parent node = leaf A
end

end

@testset "Test discrete likelihood, fixed topology" begin

# test on a tree
#=
likelihood calculated in R using a fixed Q matrix, first with ace() then
with fitDiscrete(), then with fitMK(). problem: they give different results,
see http://blog.phytools.org/2015/09/the-difference-between-different.html
- ace: misses log(#states) in its log-likelihood
- fitDiscrete in geiger: uses empirical prior at root, not stationary dist,
  but "lik" object is very flexible
- fitMk is correct. also great for 2 correlated binary traits

library(ape)
mytree = read.tree(text = "(A:3.0,(B:2.0,(C:1.0,D:1.0):1.0):1.0);")
states = c(1,1,2,2)
names(states)  = mytree$tip.label
fitER = ace(states, mytree, model="ER", type="discrete")
print(fitER$loglik, digits=17) # log-likelihood = -1.9706530878326345
print(fitER$loglik - log(2), digits=17) #         -2.6638002683925799
print(fitER$rates, digits=17)  # rates = 0.3743971742794559
print(fitER$lik.anc, digits=17)# posterior probs of states at nodes: 3x2 matrix (3 internal nodes, 2 states)

library(geiger)
fitER = fitDiscrete(mytree, states, model="ER")
print(fitER$opt$q12, digits=17) # rates = 0.36836216513047726
print(fitER$opt$lnL, digits=17) # log-likelihood = -2.6626566310743804
lik = fitER$lik
lik(0.3743971742794559, root="given",root.p=c(.5,.5)) # -2.6638002630818232: same as ace + log(2)

library(phytools)
Q2 = matrix(c(-1,1,1,-1),2,2)*fitER$opt$q12
fit2 = fitMk(mytree, states, model="ER", fixedQ=Q2)
print(fit2$logLik, digits=17) # log-likelihood = -2.6638637960257574

fitER = fitDiscrete(mytree, states, model="ARD")
lik = fitER$lik
Q = c(0.29885191850718751, 0.38944304456937912) # q12, q21
lik(Q, root="given", root.p=Q[2:1]/sum(Q)) # -2.6457428692377234
lik(Q, root="flat") # -2.6447321523303113
Q = c(0.2, 0.3) # q12, q21
lik(Q, root="flat") # -2.6754091090953693 .1,.7: -3.3291679800706073
optim(Q, lik, lower=1e-8, control=list(fnscale=-1), root="flat")
# rates = 0.29993140042699212 0.38882902905265493 loglik=-2.6447247349802496

states=c(1,2,1); names(states)=c("A","B","D")
fitER = fitDiscrete(mytree, states, model="ARD"); lik = fitER$lik
lik(Q, root="flat") # -2.1207856874033491
=#

net = readTopology("(A:3.0,(B:2.0,(C:1.0,D:1.0):1.0):1.0);");
tips = Dict("A" => "lo", "B" => "lo", "C" => "hi", "D" => "hi");
m1 = EqualRatesSubstitutionModel(2,0.36836216513047726, ["lo", "hi"]);
fit1 = (@test_nowarn fitDiscrete(net, m1, tips; fixedparam=true));
@test_nowarn show(DevNull, fit1)
@test loglikelihood(fit1) ≈ -2.6638637960257574
species = ["G","C","A","B","D"]
dat1 = DataFrame(trait = ["hi","hi","lo","lo","hi"], species = species)
m2 = BinaryTraitSubstitutionModel(0.2, 0.3, ["lo", "hi"])
fit2 = (@test_nowarn fitDiscrete(net, m2, dat1; fixedparam=true))
@test fit2.trait == [[1],[1],[2],[2]]
@test loglikelihood(fit2) ≈ -2.6754091090953693
originalSTDOUT = STDOUT
redirect_stdout(open("/dev/null", "w"))
fit2 = @test_nowarn fitDiscrete(net, m2, dat1; verbose=true) # 65 iterations
redirect_stdout(originalSTDOUT)
@test fit2.model.rate ≈ [0.29993140042699212, 0.38882902905265493] atol=2e-4
@test loglikelihood(fit2) ≈ -2.6447247349802496 atol=2e-4
m2.rate = [0.2, 0.3];
dat2 = DataFrame(trait1= ["hi","hi","lo","lo","hi"], trait2=["hi",missing,"lo","hi","lo"]);
fit3 = (@test_nowarn fitDiscrete(net, m2, species, dat2; fixedparam=true))
@test fit3.loglik ≈ (-2.6754091090953693 - 2.1207856874033491)
PhyloNetworks.fit!(fit3; fixedparam=false)
@test fit3.model.rate ≈ [0.3245640354187991, 0.5079501745877728]
fit3.net = readTopology("(A,(B,(C,D):1.0):1.0);"); # no branch lengths
@test_throws ErrorException PhyloNetworks.fit!(fit3; fixedparam=false)

# test on a network, 1 hybridization
net = readTopology("(((A:4.0,(B:1.0)#H1:1.1::0.9):0.5,(C:0.6,#H1:1.0::0.1):1.0):3.0,D:5.0);")
# function below used to check that simulation proportions == likelihood
m1 = BinaryTraitSubstitutionModel([1.0, 2.0], [1,2]) # model.label = model.index
function traitprobabilities(model, net, ntraits=10)
    res, lab = randomTrait(model, net; ntraits=ntraits)
    tips = findin(lab, tipLabels(net)) # indices of tips: columns in res
    dat = DataFrame(species = lab[tips])
    tmp = StatsBase.countmap([res[i,tips] for i in 1:ntraits])
    i = 0
    prop = Float64[]
    for (k,v) in tmp
        i += 1
        dat[Symbol("x",i)] = k
        push!(prop, v/ntraits)
    end
    npatterns = i
    lik = Float64[]
    for i in 1:npatterns
        fit = fitDiscrete(net, model, dat[[:species, Symbol("x",i)]]; fixedparam=true)
        push!(lik, fit.loglik)
    end
    return dat, prop, lik
end
#=
using PhyloNetworks, StatsBase, DataFrames
d, p, ll = traitprobabilities(m1, net, 100000000);
all(isapprox.(log.(p), ll, atol=1e-3)) # true
hcat(log.(p), ll)
 -1.62173  -1.62184
 -3.00805  -3.00807
 -4.39506  -4.39436
 -3.00747  -3.0082
 -3.70119  -3.70121
 -3.00759  -3.0082
 -2.31516  -2.31505
 -2.31554  -2.31499
 -3.0083   -3.0082
 -3.008    -3.0082
 -2.31475  -2.31505
 -3.702    -3.70135
 -3.00836  -3.00813
 -3.70033  -3.70121
 -2.31546  -2.31499
 -3.70124  -3.70135
=#
d = DataFrame(species=["D","C","B","A"], x1=[1,1,1,1], x2=[1,2,2,1], x3=[2,2,2,2], x4=[1,1,2,2],
    x5=[2,2,2,1], x6=[2,2,1,1], x7=[1,1,2,1], x8=[2,1,1,1], x9=[2,1,2,1], x10=[1,2,1,2],
    x11=[1,2,1,1], x12=[2,2,1,2], x13=[2,1,1,2], x14=[1,2,2,2], x15=[1,1,1,2], x16=[2,1,2,2])
lik = Float64[]
for i in 1:16
    fit = fitDiscrete(net, m1, d[[:species, Symbol("x",i)]]; fixedparam=true)
    push!(lik, fit.loglik)
end
@test lik ≈ [-1.6218387598967712, -3.008066347196894, -4.3943604143403245, -3.008199100743402,
    -3.70121329832901, -3.0081981601869483, -2.315051933868397, -2.314985711030534,
    -3.0081988850020873, -3.0081983709272504, -2.3150512090547584, -3.70134532205944,
    -3.008132923628349, -3.7012134632082083, -2.3149859724945876, -3.7013460518770915]
fit1 = fitDiscrete(net, m1, d[[:species, :x6]])

# with parameter estimation
net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);")
m1 = BinaryTraitSubstitutionModel([1.0, 1.0], ["lo", "hi"])
dat = DataFrame(species=["C","A","B","D"], trait=["hi","lo","lo","hi"])
fit1 = fitDiscrete(net, m1, dat; fixedparam=true)
@test fit1.loglik ≈ -2.77132013004859
PhyloNetworks.fit!(fit1; fixedparam=false)
@test fit1.model.rate ≈ [0.2722263130324768, 0.34981109618902395] atol=1e-4
@test fit1.loglik ≈ -2.727701700695741
# for information only: function used locally to check for correct parameter estimation
function simulateManyTraits_estimate(ntraits)
    m1 = BinaryTraitSubstitutionModel([1.0, 0.5], [1,2])
    res, lab = randomTrait(m1, net; ntraits=ntraits)
    tips = findin(lab, tipLabels(net)) # indices of tips: columns in res
    dat = DataFrame(transpose(res[:,tips])); species = lab[tips]
    return fitDiscrete(net, m1, species, dat)
end
# simulateManyTraits_estimate(100000)
# α=1.1124637623451075, β=0.5604529225895175, loglik=-25587.1  with ntraits=10000
# α=0.9801472136310236, β=0.4891696992781437, loglik=-255755.6 with ntraits=100000
# time with ntraits=100000: 907.2s = 15min 7s (one single processor, no binning of traits with same pattern)

# ancestral state reconstruction - fixit!!
fit1.model.rate[1] = 0.2722263130324768;
fit1.model.rate[2] = 0.34981109618902395;
@test_throws ErrorException ancestralStateReconstruction(fit1, 4) # 1 trait, not 4: error
asr = ancestralStateReconstruction(fit1)
@test names(asr) == [:nodenumber, :nodelabel, :lo, :hi]
@test asr[:nodenumber] == collect(1:9)
@test asr[:nodelabel] == ["A","B","C","D","5","6","7","8","#H1"]
@test asr[:lo] ≈ [1.,1.,0.,0., 0.28602239466671175, 0.31945742289603263,
    0.16855042517785512, 0.7673588716207436, 0.7827758475866091]
@test asr[:hi] ≈ [0.,0.,1.,1.,0.713977605333288, 0.6805425771039674,
    0.8314495748221447, 0.23264112837925616, 0.21722415241339132]
@test fit1.postltw ≈ [-0.08356534477069566, -2.5236181051014333]
end # end of testset, fixed topology

end # of nested testsets


text:
# EXCLUDE FROM TESTING
using VectorEngine, Test

# pass a simple but mutable struct to VE device

mutable struct xm
      x::Int32
      m::Int64
end

# VE side function that modifies the struct
function pass_struct!(r::Ref{xm})
    r[].m += 1
    return
end

a = xm(1, 100)

@synchronize @veda pass_struct!(Ref(a))
@test (a.x == 1) && (a.m == 101)


text:
export monodromy_solve,
    find_start_pair,
    MonodromyOptions,
    MonodromyResult,
    is_success,
    is_heuristic_stop,
    parameters,
    permutations,
    trace,
    verify_solution_completeness

#####################
# Monodromy Options #
#####################

## Options
"""
    MonodromyOptions(; options...)

Options for [`monodromy_solve`](@ref).
"""
Base.@kwdef struct MonodromyOptions{D,GA<:Union{Nothing,GroupActions}}
    check_startsolutions::Bool = true
    group_actions::GA = nothing
    loop_finished_callback = always_false
    parameter_sampler = independent_normal
    equivalence_classes::Bool = !isnothing(group_actions)
    # stopping heuristic
    trace_test_tol::Float64 = 1e-6
    trace_test::Bool = true
    target_solutions_count::Union{Nothing,Int} = nothing
    timeout::Union{Nothing,Float64} = nothing
    min_solutions::Union{Nothing,Int} = nothing
    max_loops_no_progress::Int = 5
    reuse_loops::Symbol = :all
    permutations::Bool = false
    # unique points options
    distance::D = EuclideanNorm()
    unique_points_atol::Union{Nothing,Float64} = nothing
    unique_points_rtol::Union{Nothing,Float64} = nothing
end

"""
    independent_normal(p::AbstractVector)

Sample a vector where each entry is drawn independently from the complex Normal distribution
by calling [`randn(ComplexF64)`](@ref).

    independent_normal(L::LinearSubspace)

Creates a random linear subspace by calling [`rand_subspace`](@ref).
"""
independent_normal(p::AbstractVector{T}) where {T} = randn(ComplexF64, length(p))
independent_normal(L::LinearSubspace) = rand_subspace(ambient_dim(L); dim = dim(L))

#############################
# Loops and Data Structures #
#############################

struct LoopTrackingJob
    id::Int
    loop_id::Int
end

#######################
# Loop data structure
#######################

struct MonodromyLoop{P<:Union{LinearSubspace{ComplexF64},Vector{ComplexF64}}}
    # p -> p₁ -> p₂ -> p
    p::P
    p₀₁::P # halfway
    p₁::P
    p₂::P
end

function MonodromyLoop(base::AbstractVector, options::MonodromyOptions)
    p = convert(Vector{ComplexF64}, base)
    p₁ = convert(Vector{ComplexF64}, options.parameter_sampler(p))
    p₂ = convert(Vector{ComplexF64}, options.parameter_sampler(p))

    MonodromyLoop(p, 0.5 .* (p₁ .- p), p₁, p₂)
end

function MonodromyLoop(base::LinearSubspace, options::MonodromyOptions)
    L = convert(LinearSubspace{ComplexF64}, base)
    # the second linear space is just a translation in order to perform a trace test
    # To still find new solutions quickly we translate the linear space
    # by a larger distance
    v = LA.rmul!(LA.normalize!(randn(ComplexF64, codim(L))), 5)
    L₀₁ = translate(L, v, Extrinsic)
    L₁ = translate(L₀₁, v, Extrinsic)
    L₂ = convert(LinearSubspace{ComplexF64}, options.parameter_sampler(L))

    MonodromyLoop(L, L₀₁, L₁, L₂)
end

##########################
## Monodromy Statistics ##
##########################

Base.@kwdef mutable struct MonodromyStatistics
    tracked_loops::Threads.Atomic{Int} = Threads.Atomic{Int}(0)
    tracking_failures::Threads.Atomic{Int} = Threads.Atomic{Int}(0)
    generated_loops::Threads.Atomic{Int} = Threads.Atomic{Int}(0)
    solutions::Vector{Int} = Int[]
    permutations::Vector{Vector{Int}} = Vector{Int}[]
end

MonodromyStatistics(nsolutions::Int) =
    MonodromyStatistics(solutions_development = [nsolutions])

function Base.show(io::IO, S::MonodromyStatistics)
    println(io, "MonodromyStatistics")
    print_fieldname(io, S, :tracked_loops)
    print_fieldname(io, S, :tracking_failures)
    print_fieldname(io, S, :solutions)
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", S::MonodromyStatistics) = S

# update routines
function loop_tracked!(stats::MonodromyStatistics)
    Threads.atomic_add!(stats.tracked_loops, 1)
    stats
end
function loop_failed!(stats::MonodromyStatistics)
    Threads.atomic_add!(stats.tracking_failures, 1)
    stats
end
function loop_finished!(stats, nsolutions)
    push!(stats.solutions, nsolutions)
end

function add_permutation!(stats::MonodromyStatistics, loop_id, start_id, end_id)
    perms = stats.permutations[loop_id]
    # ensure that perms has correct length
    while length(perms) < start_id
        push!(perms, 0)
    end
    perms[start_id] = end_id
    stats
end

function solutions_current_loop(stats::MonodromyStatistics, nsolutions)
    nsolutions - stats.solutions[end]
end
function solutions_last_loop(stats::MonodromyStatistics)
    length(stats.solutions) > 1 ? stats.solutions[end] - stats.solutions[end-1] : 0
end

function loops_no_change(stats::MonodromyStatistics, nsolutions)
    k = 0
    for i = length(stats.solutions):-1:1
        stats.solutions[i] == nsolutions || break
        k += 1
    end
    return max(k - 1, 0)
end

@noinline function make_showvalues(stats::MonodromyStatistics; queued::Int, solutions::Int)
    [
        ("tracked loops (queued)", "$(stats.tracked_loops[]) ($queued)"),
        (
            "solutions in current (last) loop",
            "$(solutions_current_loop(stats, solutions)) ($(solutions_last_loop(stats)))",
        ),
        (
            "generated loops (no change)",
            "$(stats.generated_loops[]) ($(loops_no_change(stats, solutions)))",
        ),
    ]
end

############
## Result ##
############

"""
    MonodromyResult

The monodromy result contains the result of the [`monodromy_solve`](@ref) computation.
"""
struct MonodromyResult{P,LP}
    returncode::Symbol
    results::Vector{PathResult}
    parameters::P
    loops::Vector{MonodromyLoop{LP}}
    statistics::MonodromyStatistics
    equivalence_classes::Bool
    seed::UInt32
    trace::Union{Nothing,Float64}
end

Base.show(io::IO, ::MIME"application/prs.juno.inline", x::MonodromyResult) = x
function Base.show(io::IO, result::MonodromyResult{T}) where {T}
    println(io, "MonodromyResult")
    println(io, "="^length("MonodromyResult"))
    println(io, "• return_code → :$(result.returncode)")
    if result.equivalence_classes
        println(io, "• $(nsolutions(result)) classes of solutions (modulo group action)")
    else
        println(io, "• $(nsolutions(result)) solutions")
    end
    println(io, "• $(result.statistics.tracked_loops[]) tracked loops")
    print(io, "• random_seed → ", sprint(show, result.seed))
    if !isnothing(result.trace)
        print(io, "\n• trace → ", sprint(show, result.trace))
    end
end

TreeViews.hastreeview(::MonodromyResult) = true
TreeViews.numberofnodes(::MonodromyResult) = 7
TreeViews.treelabel(io::IO, x::MonodromyResult, ::MIME"application/prs.juno.inline") =
    print(
        io,
        "MonodromyResult",
    )
#
function TreeViews.nodelabel(
    io::IO,
    x::MonodromyResult,
    i::Int,
    ::MIME"application/prs.juno.inline",
)
    if i == 1
        if x.equivalence_classes
            print(io, "$(nsolutions(x)) solutions (modulo group action)")
        else
            print(io, "$(nsolutions(x)) solutions")
        end
    elseif i == 2
        print(io, "Return code")
    elseif i == 3
        print(io, "Tracked loops")
    elseif i == 4
        print(io, "Statistics")
    elseif i == 5
        print(io, "Parameters")
    elseif i == 6
        print(io, "Seed")
    elseif i == 7
        print(io, "Trace")
    end
end
function TreeViews.treenode(r::MonodromyResult, i::Integer)
    if i == 1
        return r.results
    elseif i == 2
        return r.returncode
    elseif i == 3
        return r.statistics.tracked_loops[]
    elseif i == 4
        return r.statistics
    elseif i == 5
        return r.parameters
    elseif i == 6
        return r.seed
    elseif i == 7
        return something(r.trace, missing)
    end
    missing
end

"""
    is_success(result::MonodromyResult)

Returns true if the monodromy computation achieved its target solution count.
"""
is_success(result::MonodromyResult) = result.returncode == :success

"""
    is_heuristic_stop(result::MonodromyResult)

Returns true if the monodromy computation stopped due to the heuristic.
"""
is_heuristic_stop(result::MonodromyResult) = result.returncode == :heuristic_stop

"""
    solutions(result::MonodromyResult)

Return all solutions.
"""
solutions(r::MonodromyResult) = solutions(r.results)

"""
    nsolutions(result::MonodromyResult)

Returns the number solutions of the `result`.
"""
nsolutions(r::MonodromyResult) = length(r.results)

"""
    results(result::MonodromyResult)

Returns the computed [`PathResult`](@ref)s.
"""
results(r::MonodromyResult) = r.results

"""
    nresults(result::MonodromyResult)

Returns the number of results computed.
"""
nresults(r::MonodromyResult) = length(r.results)

"""
    parameters(result::MonodromyResult)

Return the parameters corresponding to the given result `r`.
"""
ModelKit.parameters(r::MonodromyResult) = r.parameters

"""
    seed(result::MonodromyResult)

Return the random seed used for the computations.
"""
seed(r::MonodromyResult) = r.seed

"""
    trace(result::MonodromyResult)

Return the result of the trace test computed during the monodromy.
"""
trace(r::MonodromyResult) = r.trace

"""
    permutations(r::MonodromyResult; reduced=true)

Return the permutations of the solutions that are induced by tracking over the loops.
If `reduced=false`, then all permutations are returned.
If `reduced=true` then permutations without repetitions are returned.

If a solution was not tracked in the loop, then the corresponding entry is 0.

Example: monodromy loop for a varying line that intersects two circles.
```julia
using LinearAlgebra
@var x[1:2] a b c
c1 = (x - [2, 0]) ⋅ (x - [2, 0]) - 1
c2 = (x - [-2, 0]) ⋅ (x - [-2, 0]) - 1
F = [c1 * c2; a * x[1] + b * x[2] - c]
S = monodromy_solve(F, [[1, 0]], [1, 1, 1], parameters = [a, b, c], permutations = true)

permutations(S)
```

will return

```julia
2×2 Array{Int64,2}:
 1  2
 2  1
```

and `permutations(S, reduced = false)` returns

```julia
2×12 Array{Int64,2}:
 1  2  2  1  1  …  1  2  1  1  1
 2  1  1  2  2     2  1  2  2  2
```

"""
function permutations(r::MonodromyResult; reduced::Bool = true)
    π = r.statistics.permutations
    if reduced
        π = unique(π)
    end
    N = nresults(r)
    A = zeros(Int, N, length(π))
    for (j, πⱼ) in enumerate(π), i = 1:N
        A[i, j] = πⱼ[i]
    end

    A
end

#####################
## monodromy solve ##
#####################

mutable struct MonodromySolver{
    Tracker<:EndgameTracker,
    P,
    UP<:UniquePoints,
    MO<:MonodromyOptions,
}
    trackers::Vector{Tracker}
    loops::Vector{MonodromyLoop{P}}
    unique_points::UP
    unique_points_lock::ReentrantLock
    options::MO
    statistics::MonodromyStatistics
    # We save the sums of the solutions for three
    # values to check if the trace is colinear
    # We augment the sums by 1 in each coordinate to
    # make this work for the two variables case
    # (n + 1) × 3 matrix
    trace::Matrix{ComplexF64}
    trace_lock::ReentrantLock
end

function MonodromySolver(
    F::Union{System,AbstractSystem},
    parameters;
    intrinsic = nothing,
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    tracker_options = TrackerOptions(),
    options = MonodromyOptions(),
)
    if !isa(options, MonodromyOptions)
        options = MonodromyOptions(; options...)
    end

    if parameters isa LinearSubspace
        H = linear_subspace_homotopy(
            F,
            parameters,
            parameters;
            compile = compile,
            intrinsic = intrinsic,
        )
        P = LinearSubspace{ComplexF64}
    else
        H = parameter_homotopy(F; generic_parameters = parameters, compile = compile)
        P = Vector{ComplexF64}
    end
    egtracker = EndgameTracker(
        H;
        tracker_options = tracker_options,
        options = EndgameOptions(endgame_start = 0.0),
    )
    trackers = [egtracker]
    group_actions = options.equivalence_classes ? options.group_actions : nothing
    x₀ = zeros(ComplexF64, size(F, 2))
    #TODO: This is wrong for IntrinsicSubspaceHomotopy of projective varieties
    if !isnothing(group_actions) && (egtracker.tracker.homotopy isa AffineChartHomotopy)
        group_actions = let H = egtracker.tracker.homotopy, group_actions = group_actions
            (cb, s) -> apply_actions(v -> cb(set_solution!(v, H, v, 0)), group_actions, s)
        end
    end

    unique_points =
        UniquePoints(x₀, 1; metric = options.distance, group_actions = group_actions)

    trace = zeros(ComplexF64, length(x₀) + 1, 3)
    MonodromySolver(
        trackers,
        MonodromyLoop{P}[],
        unique_points,
        ReentrantLock(),
        options,
        MonodromyStatistics(),
        trace,
        ReentrantLock(),
    )
end

function add_loop!(MS::MonodromySolver, p)
    push!(MS.loops, MonodromyLoop(p, MS.options))
    Threads.atomic_add!(MS.statistics.generated_loops, 1)
    if MS.options.permutations
        push!(MS.statistics.permutations, zeros(Int, length(MS.unique_points)))
    end
    MS
end
loop(MS::MonodromySolver, i) = MS.loops[i]
nloops(MS::MonodromySolver) = length(MS.loops)

function reset_trace!(MS::MonodromySolver)
    MS.trace .= 0
    MS.trace[end, :] .= 1
    MS
end

# Check if trace is colinear
function trace_colinearity(MS)
    σ = LA.svdvals(MS.trace)
    σ[3] / σ[1]
end

"""
    monodromy_solve(F, [sols, p]; options..., tracker_options = TrackerOptions())

Solve a polynomial system `F(x;p)` with specified parameters and initial solutions `sols`
by monodromy techniques. This makes loops in the parameter space of `F` to find new solutions.
If `F` the parameters `p` only occur *linearly* in `F` it is eventually possible to compute
a *start pair* ``(x₀, p₀)`` automatically. In this case `sols` and `p` can be omitted and
the automatically generated parameters can be obtained with the [`parameters`](@ref) function
from the [`MonodromyResult`](@ref).

    monodromy_solve(F, [sols, L]; dim, codim, intrinsic = nothing, options...,
                                  tracker_options = TrackerOptions())

Solve the polynomial system `[F(x); L(x)] = 0` where `L` is a `[LinearSubspace`](@ref).
If `sols` and `L` are not provided it is necesary to provide `dim` or `codim` where `(co)dim`
is the expected (co)dimension of a component of `V(F)`.
See also [`linear_subspace_homotopy`](@ref) for the `intrinsic` option.

## Options

* `catch_interrupt = true`: If true catches interruptions (e.g. issued by pressing Ctrl-C)
  and returns the partial result.
* `check_startsolutions = true`: If `true`, we do a Newton step for each entry of `sols`for
  checking if it is a valid startsolutions. Solutions which are not valid are sorted out.
* `compile = $(COMPILE_DEFAULT[])`: If `true` then a `System` (resp. `Homotopy`) is compiled
  to a straight line program ([`CompiledSystem`](@ref) resp. [`CompiledHomotopy`](@ref))
  for evaluation. This induces a compilation overhead. If `false` then the generated program
  is only interpreted ([`InterpretedSystem`](@ref) resp. [`InterpretedHomotopy`](@ref)).
  This is slower than the compiled version, but does not introduce compilation overhead.
* `distance = EuclideanNorm()`: The distance function used for [`UniquePoints`](@ref).
* `loop_finished_callback = always_false`: A callback to end the computation. This function is
  called with all current [`PathResult`](@ref)s after a loop is exhausted.
  2 arguments. Return `true` if the compuation should be stopped.
* `equivalence_classes=true`: This only applies if there is at least one group action
  supplied. We then consider two solutions in the same equivalence class if we can transform
  one to the other by the supplied group actions. We only track one solution per equivalence
  class.
* `group_action = nothing`: A function taking one solution and returning other solutions if
  there is a constructive way to obtain them, e.g. by symmetry.
* `group_actions = nothing`: If there is more than one group action you can use this to chain
  the application of them. For example if you have two group actions `foo` and `bar` you can
  set `group_actions=[foo, bar]`. See [`GroupActions`](@ref) for details regarding the
  application rules.
* `max_loops_no_progress = 5`: The maximal number of iterations (i.e. loops generated)
  without any progress.
* `min_solutions`: The minimal number of solutions before a stopping heuristic is applied.
  By default no lower limit is enforced.
* `parameter_sampler = independent_normal`: A function taking the parameter `p` and
  returning a new random parameter `q`. By default each entry of the parameter vector
  is drawn independently from Normal distribution.
* `permutations = false`: Whether to keep track of the permutations induced by the loops.
* `resuse_loops = :all`: Strategy to reuse other loops for new found solutions. `:all`
  propagates a new solution through all other loops, `:random` picks a random loop, `:none`
  doesn't reuse a loop.
* `target_solutions_count`: The computation is stopped if this number of solutions is
  reached.
* `threading = true`: Enable multithreading of the path tracking.
* `timeout`: The maximal number of *seconds* the computation is allowed to run.
* `trace_test = true`: If `true` a trace test is performed to check whether all solutions
  are found. This is only applicable if monodromy is performed with a linear subspace.
  See also [`trace_test`](@ref).
* `trace_test_tol = 1e-10`: The tolerance for the trace test to be successfull.
  The trace is divided by the number of solutions before compared to the trace_test_tol.
* `unique_points_rtol`: the relative tolerance for [`unique_points`](@ref).
* `unique_points_atol`: the absolute tolerance for [`unique_points`](@ref).
"""
function monodromy_solve(F::Vector{<:ModelKit.Expression}, args...; parameters, kwargs...)
    monodromy_solve(System(F; parameters = parameters), args...; kwargs...)
end
function monodromy_solve(
    F::AbstractVector{<:MP.AbstractPolynomial},
    args...;
    parameters,
    variables = setdiff(MP.variables(F), parameters),
    variable_ordering = variables,
    kwargs...,
)
    sys = System(F, variables = variable_ordering, parameters = parameters)
    monodromy_solve(sys, args...; kwargs...)
end
function monodromy_solve(
    F::Union{System,AbstractSystem},
    args...;
    seed = rand(UInt32),
    tracker_options = TrackerOptions(),
    show_progress::Bool = true,
    threading::Bool = Threads.nthreads() > 1,
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    catch_interrupt::Bool = true,
    dim = nothing,
    codim = nothing,
    intrinsic = nothing,
    # monodromy options
    options = nothing,
    group_action = nothing,
    group_actions = isnothing(group_action) ? nothing : GroupActions(group_action),
    kwargs...,
)
    if isnothing(options)
        if !(group_actions isa GroupActions) && !isnothing(group_actions)
            group_actions = GroupActions(group_actions)
        end
        options = MonodromyOptions(; group_actions = group_actions, kwargs...)
    end

    if !isnothing(seed)
        Random.seed!(seed)
    end
    if length(args) == 0
        start_pair = find_start_pair(fixed(F; compile = compile))
        if isnothing(start_pair)
            error(
                "Cannot compute a start pair (x, p) using `find_start_pair(F)`." *
                " You need to explicitly pass a start pair.",
            )
        end
        x, p = start_pair
        S = [x]
        # if we have no parameters, then we intersect with a linear subspace.
        # For this we need to know the intended dimension or codimension.
        # We don't guess dim to catch the case that the user just forgot to pass
        # parameters
        if isnothing(p)
            if isnothing(dim) && isnothing(codim)
                error(
                    "Given system doesn't have any parameters. If you intended to intersect " *
                    "with a linear subspace it is necessary to provide a " *
                    "dimension (`dim`) or codimension (`codim`) of the component of interest.",
                )
            else
                projective = is_homogeneous(System(F))
                if !isnothing(codim)
                    codim += projective
                end
                p = rand_subspace(
                    x;
                    dim = codim,
                    codim = dim,
                    affine = !is_homogeneous(System(F)),
                )
            end
        end
    else
        S, p = args
    end
    MS = MonodromySolver(
        F,
        p;
        compile = compile,
        options = options,
        tracker_options = tracker_options,
        intrinsic = intrinsic,
    )

    monodromy_solve(
        MS,
        S,
        p,
        seed;
        show_progress = show_progress,
        threading = threading,
        catch_interrupt = catch_interrupt,
    )
end

"""
    find_start_pair(F; max_tries = 100, atol = 0.0, rtol = 1e-12)

Try to find a pair `(x,p)` for the system `F` such that `F(x,p) = 0` by randomly sampling
a pair `(x₀, p₀)` and performing Newton's method in variable *and* parameter space.
"""
function find_start_pair(
    F::System;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    max_tries = 1_000,
    kwargs...,
)
    if nparameters(F) != 0
        start_pair = find_start_pair_linear_in_params(F; kwargs...)
        if !isnothing(start_pair)
            return start_pair
        end
    end
    find_start_pair(fixed(F; compile = compile); max_tries = max_tries, kwargs...)
end
function find_start_pair(
    F::AbstractSystem;
    max_tries::Int = 1_000,
    atol::Float64 = 0.0,
    rtol::Float64 = 1e-12,
)
    if nparameters(F) == 0
        SF = F
    else
        SF = StartPairSystem(F)
    end
    m = nparameters(F)
    n = nvariables(F)
    c = NewtonCache(SF)
    for i = 1:max_tries
        xp₀ = randn(ComplexF64, size(SF, 2))
        res = newton(SF, xp₀, nothing, InfNorm(), c, atol = atol, rtol = rtol)
        if is_success(res)
            x, p = res.x[1:n], res.x[n+1:end]
            res2 = newton(F, x, p, InfNorm())
            if is_success(res2)
                return res2.x, isempty(p) ? nothing : p
            end
        end
    end
    nothing
end
function find_start_pair_linear_in_params(F; kwargs...)
    x₀ = randn(ComplexF64, nvariables(F))
    F_x₀ = F(x₀, parameters(F))
    all(!ModelKit.is_number, F_x₀) || return nothing

    Ab = construct_linear_system(F_x₀, parameters(F))
    if !isnothing(Ab)
        A, b = Ab
        m, n = size(A)
        m ≤ n || return nothing
        # if b is 0 then we compute the a basis vector of the nullspace
        if iszero(b)
            # don't want to have 0 as a parameter
            m == n && return nothing
            p₀ = LA.nullspace(A)[:, 1]
        else
            p₀ = LA.qr(A, Val(true)) \ b
        end
        return x₀, p₀
    else
        return nothing
    end
end

function construct_linear_system(F, x)
    A = zeros(ComplexF64, length(F), length(x))
    b = zeros(ComplexF64, length(F))

    try
        for (i, f) in enumerate(F)
            M, coeffs = exponents_coefficients(f, x)
            for (j, exp) in enumerate(eachcol(M))
                var_found = false
                for (k, d) in enumerate(exp)
                    if d < 0 || d > 1 || (d == 1 && var_found)
                        return nothing
                    end
                    if d == 1
                        var_found = true
                        A[i, k] = coeffs[j]
                    end
                end
                if !var_found
                    b[i] = -coeffs[j]
                end
            end
        end
    catch err
        if err isa ModelKit.PolynomialError
            return nothing
        else
            rethrow(err)
        end
    end
    A, b
end

monodromy_solve(MS::MonodromySolver, x::AbstractVector{<:Number}, p, seed; kwargs...) =
    monodromy_solve(MS, [x], p, seed; kwargs...)

function monodromy_solve(
    MS::MonodromySolver,
    X::AbstractVector{<:AbstractVector{<:Number}},
    p,
    seed;
    show_progress::Bool = true,
    threading::Bool = Threads.nthreads() > 1,
    catch_interrupt::Bool = true,
)
    if !show_progress
        progress = nothing
    else
        if MS.options.equivalence_classes
            desc = "Solutions (modulo group action) found:"
        else
            desc = "Solutions found:"
        end
        progress = ProgressMeter.ProgressUnknown(; dt = 0.4, desc = desc, output = stdout)
        progress.tlast += 0.3
    end
    MS.statistics = MonodromyStatistics()
    empty!(MS.unique_points)
    reset_trace!(MS)
    results = check_start_solutions(MS, X, p)
    if isempty(results)
        @warn "None of the provided solutions is a valid start solution."
        retcode = :invalid_startvalue
    else
        retcode = :default
        try
            # convert to complex base
            if p isa LinearSubspace
                cp = convert(LinearSubspace{ComplexF64}, p)
            else
                cp = convert(Vector{ComplexF64}, p)
            end
            if threading
                retcode = threaded_monodromy_solve!(results, MS, cp, seed, progress)
            else
                retcode = serial_monodromy_solve!(results, MS, cp, seed, progress)
            end
        catch e
            if !catch_interrupt ||
               !(
                isa(e, InterruptException) ||
                (isa(e, TaskFailedException) && isa(e.task.exception, InterruptException))
            )
                rethrow(e)
            end
            retcode = :interrupted
        end
    end

    MonodromyResult(
        retcode,
        results,
        p,
        MS.loops,
        MS.statistics,
        MS.options.equivalence_classes,
        seed,
        p isa LinearSubspace ? trace_colinearity(MS) : nothing,
    )
end

function check_start_solutions(MS::MonodromySolver, X, p)
    tracker = MS.trackers[1]
    parameters!(tracker, p, p)
    results = PathResult[]
    for x in X
        res = track(tracker, x)
        if is_success(res)
            _, added = add!(MS, res, length(results) + 1)
            if added
                push!(results, res)
            end
        end
    end

    results
end

function add!(MS::MonodromySolver, res::PathResult, id)
    if isnothing(MS.options.unique_points_rtol)
        rtol = uniqueness_rtol(res)
    else
        rtol = MS.options.unique_points_rtol
    end
    if isnothing(MS.options.unique_points_atol)
        atol = 1e-14
    else
        atol = MS.options.unique_points_atol
    end

    add!(MS.unique_points, solution(res), id; atol = atol, rtol = rtol)
end

function serial_monodromy_solve!(
    results::Vector{PathResult},
    MS::MonodromySolver,
    p,
    seed,
    progress,
)
    queue = LoopTrackingJob[]
    tracker = MS.trackers[1]
    t₀ = time()
    retcode = :in_progress
    stats = MS.statistics
    while true
        loop_finished!(stats, length(results))

        if MS.options.loop_finished_callback(results)
            retcode = :terminated_callback
            @goto _return
        end
        if p isa LinearSubspace &&
           nloops(MS) > 0 &&
           MS.options.trace_test &&
           trace_colinearity(MS) < MS.options.trace_test_tol
            retcode = :success
            @goto _return
        end
        if isnothing(MS.options.target_solutions_count) &&
           loops_no_change(stats, length(results)) ≥ MS.options.max_loops_no_progress
            retcode = :heuristic_stop
            @goto _return
        end
        if length(results) == something(MS.options.target_solutions_count, -1)
            retcode = :success
            @goto _return
        end

        add_loop!(MS, p)
        reset_trace!(MS)
        # schedule all jobs
        new_loop_id = nloops(MS)
        for i = 1:length(results)
            push!(queue, LoopTrackingJob(i, new_loop_id))
        end

        while !isempty(queue)
            job = popfirst!(queue)
            res = track(
                tracker,
                results[job.id],
                loop(MS, job.loop_id),
                MS.trace,
                MS.trace_lock;
                collect_trace = MS.options.trace_test && nloops(MS) == job.loop_id,
            )
            if !isnothing(res)
                loop_tracked!(stats)

                # 1) check whether solutions already exists
                id, got_added = add!(MS, res, length(results) + 1)

                if MS.options.permutations
                    add_permutation!(stats, job.loop_id, job.id, id)
                end

                got_added || @goto _update
                # 2) doesn't exist, so add to results
                push!(results, res)

                # 3) schedule on same loop again
                push!(queue, LoopTrackingJob(id, job.loop_id))

                # 4) schedule on other loops
                if MS.options.reuse_loops == :all
                    for k = 1:nloops(MS)
                        k != job.loop_id || continue
                        push!(queue, LoopTrackingJob(id, k))
                    end
                elseif MS.options.reuse_loops == :random && nloops(MS) ≥ 2
                    k = rand(2:nloops(MS))
                    if k ≤ job.loop_id
                        k -= 1
                    end
                    push!(queue, LoopTrackingJob(id, k))
                end
            else
                loop_failed!(stats)
                if MS.options.permutations
                    add_permutation!(stats, job.loop_id, job.id, 0)
                end
            end
            # Update progress
            @label _update
            update_progress!(
                progress,
                stats;
                solutions = length(results),
                queued = length(queue),
            )

            if length(results) == something(MS.options.target_solutions_count, -1) &&
               # only terminate after a completed loop to ensure that we collect proper
               # permutation informations
               !MS.options.permutations
                retcode = :success
                @goto _return
            elseif !isnothing(MS.options.timeout) &&
                   time() - t₀ > (MS.options.timeout::Float64)
                retcode = :timeout
                @goto _return
            end
        end
    end

    @label _return

    update_progress!(
        progress,
        stats;
        finish = true,
        solutions = length(results),
        queued = length(queue),
    )

    return retcode
end

function threaded_monodromy_solve!(
    results::Vector{PathResult},
    MS::MonodromySolver,
    p,
    seed,
    progress,
)
    queue = Channel{LoopTrackingJob}(Inf)
    Threads.resize_nthreads!(MS.trackers)
    data_lock = ReentrantLock()
    t₀ = time()
    retcode = :in_progress
    stats = MS.statistics
    notify_lock = ReentrantLock()
    cond_queue_emptied = Threads.Condition(notify_lock)
    workers_idle = fill(true, Threads.nthreads())
    interrupted = Ref(false)
    queued = Ref(0)
    try
        workers = map(enumerate(MS.trackers)) do (tid, tracker)
            @tspawnat tid begin
                for job in queue
                    stop_queue = false
                    Base.@lock notify_lock begin
                        if interrupted[]
                            # mark worker as idle
                            workers_idle[tid] = true
                            # if all idle notify that we want a new low
                            if all(workers_idle)
                                notify(cond_queue_emptied)
                            end
                            stop_queue = true
                        else
                            workers_idle[tid] = false
                        end
                    end
                    if stop_queue
                        break
                    end
                    res = track(
                        tracker,
                        results[job.id],
                        loop(MS, job.loop_id),
                        MS.trace,
                        MS.trace_lock;
                        collect_trace = MS.options.trace_test && nloops(MS) == job.loop_id,
                    )
                    if !isnothing(res)
                        loop_tracked!(stats)

                        # 1) check whether solutions already exists
                        lock(data_lock)
                        id, got_added = add!(MS, res, length(results) + 1)

                        if MS.options.permutations
                            add_permutation!(stats, job.loop_id, job.id, id)
                        end

                        if !got_added
                            unlock(data_lock)
                            @goto _update
                        end
                        # 2) doesn't exist, so add to results
                        push!(results, res)
                        unlock(data_lock)

                        # 3) schedule on same loop again
                        push!(queue, LoopTrackingJob(id, job.loop_id))

                        # 4) schedule on other loops
                        if MS.options.reuse_loops == :all
                            for k = 1:nloops(MS)
                                k != job.loop_id || continue
                                push!(queue, LoopTrackingJob(id, k))
                            end
                        elseif MS.options.reuse_loops == :random && nloops(MS) ≥ 2
                            k = rand(2:nloops(MS))
                            if k ≤ job.loop_id
                                k -= 1
                            end
                            push!(queue, LoopTrackingJob(id, k))
                        end
                    else
                        loop_failed!(stats)
                        if MS.options.permutations
                            add_permutation!(stats, job.loop_id, job.id, 0)
                        end
                    end
                    # Update progress
                    @label _update
                    update_progress!(
                        progress,
                        stats;
                        solutions = length(results),
                        queued = Base.n_avail(queue),
                    )

                    # mark worker as idle
                    Base.@lock notify_lock begin
                        workers_idle[tid] = true
                        # if queue is empty, check whether all other are also waiting
                        if !isready(queue) && all(workers_idle)
                            notify(cond_queue_emptied)
                        end
                    end

                    if length(results) ==
                       something(MS.options.target_solutions_count, -1) &&
                       # only terminate after a completed loop to ensure that we collect proper
                       # permutation informations
                       !MS.options.permutations
                        retcode = :success
                        Base.@lock notify_lock begin
                            interrupted[] = true
                        end

                    elseif !isnothing(MS.options.timeout) &&
                           time() - t₀ > (MS.options.timeout::Float64)
                        retcode = :timeout
                        Base.@lock notify_lock begin
                            interrupted[] = true
                        end
                    end
                end
            end
        end

        t = Threads.@spawn begin
            Base.@lock notify_lock begin
                while true
                    if interrupted[]
                        break
                    end
                    loop_finished!(stats, length(results))

                    if MS.options.loop_finished_callback(results)
                        retcode = :terminated_callback
                        break
                    end
                    if loops_no_change(stats, length(results)) ≥
                       MS.options.max_loops_no_progress
                        retcode = :heuristic_stop
                        break
                    end
                    if length(results) == something(MS.options.target_solutions_count, -1)
                        retcode = :success
                        break
                    end
                    if p isa LinearSubspace &&
                       nloops(MS) > 0 &&
                       MS.options.trace_test &&
                       trace_colinearity(MS) < MS.options.trace_test_tol
                        retcode = :success
                        break
                    end

                    add_loop!(MS, p)
                    reset_trace!(MS)
                    # schedule all jobs
                    new_loop_id = nloops(MS)
                    for i = 1:length(results)
                        push!(queue, LoopTrackingJob(i, new_loop_id))
                    end

                    wait(cond_queue_emptied)
                    retcode == :in_progress || break
                end
            end
        end

        wait(t)
    catch e
        close(queue)
        interrupted[] = true
        rethrow(e)
    finally
        queued[] = Base.n_avail(queue)
        close(queue)
    end

    update_progress!(
        progress,
        stats;
        finish = true,
        solutions = length(results),
        queued = queued[],
    )

    return retcode
end

function uniqueness_rtol(res::PathResult)
    # we should have a region of uniqueness of 0.5inv(ω) * norm(x)
    # the norm(x) factors comes from the fact that we computed with a
    # weighted norm.
    # To be more pessimistic we only consider 0.25inv(ω)^2
    # and require at most sqrt(res.accuracy) and at least 1e-14.
    clamp(0.25 * inv(res.ω)^2, 1e-14, sqrt(res.accuracy))
end

"""
    track(tracker, x, edge::LoopEdge, loop::MonodromyLoop, stats::MonodromyStatistics)

Track `x` along the edge `edge` in the loop `loop` using `tracker`. Record statistics
in `stats`.
"""
function track(
    egtracker::EndgameTracker,
    r::PathResult,
    loop::MonodromyLoop,
    trace,
    trace_lock;
    collect_trace::Bool = false,
)
    tracker = egtracker.tracker

    if loop.p isa LinearSubspace && collect_trace
        parameters!(tracker, loop.p, loop.p₀₁)
        x₀ = solution(r)
        retcode =
            track!(tracker, x₀; ω = r.ω, μ = r.μ, extended_precision = r.extended_precision)
        x₀₁ = solution(tracker)
        if !is_success(retcode)
            return nothing
        end

        parameters!(tracker, loop.p₀₁, loop.p₁)
        retcode = track!(
            tracker,
            x₀₁;
            ω = tracker.state.ω,
            μ = tracker.state.μ,
            extended_precision = tracker.state.extended_prec,
        )
        x₁ = solution(tracker)
        if !is_success(retcode)
            return nothing
        end

        Base.@lock trace_lock begin
            for i = 1:length(x₀)
                trace[i, 1] += x₀[i]
                trace[i, 2] += x₀₁[i]
                trace[i, 3] += x₁[i]
            end
        end
    else
        parameters!(tracker, loop.p, loop.p₁)
        retcode = track!(
            tracker,
            solution(r);
            ω = r.ω,
            μ = r.μ,
            extended_precision = r.extended_precision,
        )
        if !is_success(retcode)
            return nothing
        end
        x₁ = solution(tracker)
    end

    parameters!(tracker, loop.p₁, loop.p₂)
    retcode = track!(
        tracker,
        x₁;
        ω = tracker.state.ω,
        μ = tracker.state.μ,
        extended_precision = tracker.state.extended_prec,
    )
    if !is_success(retcode)
        return nothing
    end

    x₂ = solution(tracker)
    parameters!(tracker, loop.p₂, loop.p)
    result = track(
        egtracker,
        x₂;
        ω = tracker.state.ω,
        μ = tracker.state.μ,
        extended_precision = tracker.state.extended_prec,
    )
    if is_success(result)
        result
    else
        nothing
    end
end

update_progress!(progress::Nothing, stats::MonodromyStatistics; kwargs...) = nothing
function update_progress!(
    progress,
    stats::MonodromyStatistics;
    queued::Int,
    solutions::Int,
    finish::Bool = false,
)
    if finish
        ProgressMeter.update!(progress, solutions)
        showvalues = make_showvalues(stats; queued = queued, solutions = solutions)
        ProgressMeter.finish!(progress; showvalues = showvalues)
    elseif time() > progress.tlast + progress.dt
        showvalues = make_showvalues(stats; queued = queued, solutions = solutions)
        ProgressMeter.update!(progress, solutions, showvalues = showvalues)
    end
    # yield()
    nothing
end

## Parameter homotopy constructor
solve(F, R::MonodromyResult; target_parameters, kwargs...) = solve(
    F,
    solutions(R);
    start_parameters = parameters(R),
    target_parameters = target_parameters,
    kwargs...,
)

##################
## VERIFICATION ##
##################
"""
    verify_solution_completeness(F::System, monodromy_result; options...)
    verify_solution_completeness(F::System, solutions, parameters;
        trace_tol = 1e-14,
        show_progress = true,
        compile = COMPILE_DEFAULT[],
        monodromy_options = (compile = compile,),
        parameter_homotopy_options = (compile = compile,),
    )

Verify that a monodromy computation found all solutions by [`monodromy_solve`](@ref).
This uses the trace test described in [^dCR17] and [^LRS18].
The trace is a numerical value which is 0 if all solutions are found, for this the
`trace_tol` keyword argument is used. The function returns `nothing` if some computation
couldn't be carried out. Otherwise returns a boolean. Note that this function requires the
computation of solutions to another polynomial system using monodromy. This routine can
return `false` although all solutions are found if this additional solution set is not
complete.

### Example
```julia
@var x y a b c;
f = x^2+y^2-1;
l = a*x+b*y+c;
sys = System([f, l]; parameters = [a, b, c]);
mres = monodromy_solve(sys, [-0.6-0.8im, -1.2+0.4im], [1,2,3]);
show(mres);
verify_solution_completeness(sys, mres)
```
```
MonodromyResult
==================================
• 2 solutions (0 real)
• return code → heuristic_stop
• 44 tracked paths
• seed → 367230
julia> verify_solution_completeness(sys, mres)
[ Info: Certify provided solutions...
[ Info: Got 2 dinstinct solutions.
[ Info: Compute additional witnesses for completeness...
┌ Info: MonodromyResult
│ ===============
│ • return_code → :heuristic_stop
│ • 4 solutions
│ • 28 tracked loops
└ • random_seed → 0x21e7406a
[ Info: Certify additional witnesses...
[ Info: Computed 2 additional witnesses
[ Info: Compute trace using two parameter homotopies...
[ Info: Norm of trace: 9.33238819760471e-17
true
```

[^dCR17]:
    del Campo, Abraham Martín, and Jose Israel Rodriguez.
    "Critical points via monodromy and local methods."
    Journal of Symbolic Computation 79 (2017): 559-574.

[^LRS18]:
    Leykin, Anton, Jose Israel Rodriguez, and Frank Sottile. "Trace test."
    Arnold Mathematical Journal 4.1 (2018): 113-125.
"""
verify_solution_completeness(F::System, mres::MonodromyResult; kwargs...) =
    verify_solution_completeness(F, solutions(mres), parameters(mres); kwargs...)
function verify_solution_completeness(
    F::System,
    sols::AbstractVector{<:AbstractVector},
    q::AbstractVector;
    show_progress = true,
    compile = COMPILE_DEFAULT[],
    monodromy_options = (compile = compile,),
    parameter_homotopy_options = (compile = compile,),
    trace_tol = 1e-14,
)
    n = nvariables(F)
    m = nparameters(F)
    @unique_var t v[1:m] a[1:n] λ

    x = variables(F)
    p = parameters(F)

    verify_system = System(
        [F(x, p + λ * v); (sum(a .* x) - 1) * λ + t],
        variables = [x; λ],
        parameters = [t; p; v; a],
    )
    # perform monodromy computation to compute the additional witnesses
    if show_progress
        @info "Compute additional witnesses for completeness check..."
    end
    # use the verify_system but enforce t = 0 and start with λ ≠ 0
    # this way we stay on a different irreducible component

    # Let's compute some start solutions for some parameters
    Y, base_params = let
        # We sample some random parameters to set v = qq - q
        qq = randn(ComplexF64, length(q))
        # Its good to have more than 1 start solution, we can construct n
        # by performing a parameter homotopy to qq
        # and then  computing an `a` such that n solutions are on the linear space  a⋅x-1
        qq_res = solve(
            F,
            sols[1:min(n, length(sols))];
            start_parameters = q,
            target_parameters = qq,
            show_progress = show_progress,
            parameter_homotopy_options...,
        )
        a0 = reduce(vcat, transpose.(solutions(qq_res))) \ ones(nsolutions(qq_res))
        map(s -> [s; 1], solutions(qq_res)), [q; qq - q; a0]
    end

    # now use monodromy to find more solutions
    additional_mres = monodromy_solve(
        verify_system,
        Y,
        [0.0; base_params];
        parameter_sampler = p -> [0; randn(ComplexF64, length(p) - 1)],
        show_progress = show_progress,
        monodromy_options...,
    )
    additional_sols = solutions(additional_mres)
    if show_progress
        @info additional_mres
        @info "Computed $(length(additional_sols)) additional witnesses"
        @info "Compute trace using two parameter homotopies..."
    end

    # Perform parameter homotopy for the trace
    S = [map(s -> [s; 0], sols); additional_sols]
    γ = randn(ComplexF64)
    res1 = solve(
        verify_system,
        S;
        start_parameters = [0.0; base_params],
        target_parameters = [0.5 * γ; base_params],
        show_progress = show_progress,
        parameter_homotopy_options...,
    )
    S1 = solutions(res1)
    if length(S1) ≠ length(S)
        if show_progress
            @warn "Lost solution during parameter homotopy. Abort."
        end
        return nothing
    end

    res2 = solve(
        verify_system,
        solutions(res1);
        start_parameters = [0.5 * γ; base_params],
        target_parameters = [1.0 * γ; base_params],
        show_progress = show_progress,
        parameter_homotopy_options...,
    )
    S2 = solutions(res2)
    if length(S2) ≠ length(S)
        if show_progress
            @warn "Lost solution during parameter homotopy. Abort."
        end
        return nothing
    end

    # Compute trace now
    T = sum(S)
    T1 = sum(S1)
    T2 = sum(S2)

    M = [T T1 T2; 1 1 1]
    singvals = LA.svdvals(M)
    trace_norm = singvals[3] / singvals[1]

    if show_progress
        @info "Norm of trace: $trace_norm"
    end

    trace_norm < trace_tol
end


text: using SpectralCRsUtility, Test


text:
@testset "Test Fractional Differences Equations PECE method" begin
    differencefun(x) = 0.5*x+1
    α=0.5;x0=1;
    T=1; h=0.1
    prob = FractionalDifferenceProblem(differencefun, α, x0)
    sol = solve(prob, T, h, PECEDifference())
    @test isapprox(sol.u, [1.0, 2.5, 4.25, 5.125, 5.5625, 5.78125, 5.890625, 5.9453125, 5.97265625, 5.986328125, 5.9931640625]; atol=1e-3)
end

@testset "Test Grunwald Letnikov method for fractional difference equatoins" begin
    function sys!(du, u, p, t)
        du[1] = -0.05*u[2] - 0.05*u[3] + 0.01*tanh(u[2])
        du[2] = 0.05*u[1] + 0.02*u[2] + 0.01*tanh(u[1])
        du[3] = 0.1 - 0.2*u[3] + 0.05*u[1]*u[3] + 0.01*tanh(u[3])
    end

    prob = FractionalDifferenceSystem(sys!, 0.98, [1, -1, 0])
    result = solve(prob, 7, GL())

    @test isapprox(result, [1.0  0.0423841  0.0350316  0.0241793  0.0103323  -0.00594294  -0.0241551
    -1.0  0.0376159  0.0401587  0.0426291  0.0445987   0.0458335    0.0461619
     0.0  0.1        0.179209   0.24285    0.294188    0.335624     0.369018]; atol=1e-5)
end


text:
function ueberpruefekeineLoesung(A1,b1,A3,b3,blick1,blick2)

if size(A3,1) == 3 && size(A3,2) == 2 && size(b3,1) == 3 && size(b3,2) == 1
    y1 = inv(transpose(A1)*A1)*transpose(A1)*b1;
    if norm(transpose(A3[3,:])*y1-b3[3]) < 0.0001
        println("Die Größe eures Gleichungssystems stimmt zwar, allerdings lässt es sich immernoch eindeutig lösen.")
    else
        println("Ihr habt den Vektor korrekt verändert. Das führt dazu, dass es keine Lösung mehr hat.");
        println("Das Problem ist, dass egal, welchen Vektor man als x bei A*x von euch einsetzt, bekommt man nicht euren b-Vektor heraus. Veranschaulicht könnt ihr das in der folgenden Graphik sehen. Alle möglichen b-Vektoren für unterschiedliche x liegen in der roten Ebene. Euer Punkt jedoch liegt außerhalb.");
        println("Ändert solange die Werte für blick1 und blick2, so dass man erkennt, dass der Punkt neben der Ebene liegt.")
        U2,S2,V2 = svd(A3)
        u = U2[:,1];
        v = U2[:,2];
        s = -20*u-20*v;
        a = b3[1];
        b = b3[2];
        c = b3[3];
        plt = surface([(s[1],s[2],s[3]),(30*u[1],30*u[2],30*u[3]),(30*v[1],30*v[2],30*v[3])],alpha=.1, fillalpha=.1)
        plot3d!([a],[b],[c],marker=2,xlim=(-14,14),ylim=(-14,14),zlim=(-14,14), camera = (blick1,blick2),cbar= :none, leg = false,xlabel = "x_1",ylabel = "x_2",zlabel = "x_3")

        plt
    end
else
    println("Ihr habt noch einen Fehler bei der Anzahl der Einträge in der Matrix oder dem Vektor. Stellt noch einmal sicher, dass ihr euer lineares Gleichungssystem aus a) um eine Zeile erweitert habt. Das heißt, dass sowohl der Vektor, als auch die Matrix, um eins in die Länge wachsen muss.")
end
display(plt);
return u,v;
end


text:
export UrbanDriver

mutable struct UrbanDriver <: DriverModel{LatLonAccelDirection}
    rec::SceneRecord
    mlon::LaneFollowingDriver
    mlat::LateralDriverModel
    mlane::LaneChangeModel
    longAcc::Float64
    latAcc::Float64
    direction::Int
    a_cp_max::Float64
    color::Colorant

    function UrbanDriver(
        timestep::Float64;
        rec::SceneRecord = SceneRecord(1, timestep),
        mlon::LaneFollowingDriver=IDMDriver(),
        mlat::LateralDriverModel=ProportionalLaneTracker(),
        mlane::LaneChangeModel=MOBIL(timestep,mlon=IDMDriver()),
        longAcc::Float64 = 0.0,
        latAcc::Float64 = 0.0,
        direction::Int = 1,
        a_cp_max::Float64 = (rand()*0.4+0.5)*9.8,
        color::Colorant=RGB(0, 1, 0)
        )

        retval = new()
        retval.rec = rec
        retval.mlon = mlon
        retval.mlat = mlat
        retval.mlane = mlane
        retval.longAcc = longAcc
        retval.latAcc = latAcc
        retval.direction = direction
        retval.a_cp_max = a_cp_max
        retval.color = color
        retval
    end
end
AutomotiveDrivingModels.get_name(::UrbanDriver) = "UrbanDriver"

function set_desired_speed!(model::UrbanDriver, v_des::Float64)
    AutomotiveDrivingModels.set_desired_speed!(model.mlon, v_des)
    AutomotiveDrivingModels.set_desired_speed!(model.mlane, v_des)
    model
end

function track_longitudinal!(driver::LaneFollowingDriver, scene::Union{Scene,Frame{Entity{VehicleState, BicycleModel, Int}}}, roadway::Roadway, vehicle_index::Int, fore::NeighborLongitudinalResult)
    v_ego = scene[vehicle_index].state.v
    if fore.ind != 0
        headway, v_oth = fore.Δs, scene[fore.ind].state.v
    else
        headway, v_oth = NaN, NaN
    end
    return track_longitudinal!(driver, v_ego, v_oth, headway)
end

function AutomotiveDrivingModels.observe!(driver::UrbanDriver, scene2::Union{Scene,Frame{Entity{VehicleState, BicycleModel, Int}}}, roadway::Roadway, egoid::Int)
    scene = Scene()
    for veh in scene2
        push!(scene,Vehicle(veh.state,veh.def.def,veh.id))
    end
    AutomotiveDrivingModels.update!(driver.rec, scene)
    observe!(driver.mlane, scene, roadway, egoid)
    lane_change_action = rand(driver.mlane)

    vehicle_index = findfirst(scene, egoid)
    #lane_change_action = LaneChangeChoice(DIR_MIDDLE)
    #currentLane = scene[vehicle_index].state.posF.roadind.tag.lane

    laneoffset = get_lane_offset(lane_change_action, driver.rec, roadway, vehicle_index)
    lateral_speed = convert(Float64, get(VELFT, driver.rec, roadway, vehicle_index))

    if lane_change_action.dir == DIR_MIDDLE
        fore = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    elseif lane_change_action.dir == DIR_LEFT
        fore = get_neighbor_fore_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    else
        @assert(lane_change_action.dir == DIR_RIGHT)
        fore = get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    end

    AutomotiveDrivingModels.track_lateral!(driver.mlat, laneoffset, lateral_speed)
    #AutomotiveDrivingModels.track_longitudinal!(driver.mlon, scene, roadway, vehicle_index, fore.ind)
    roadind = scene[vehicle_index].state.posF.roadind
    max_k,distance = get_max_curvature(roadind, roadway, 25.0,direction = driver.direction)
    #println(max_k,distance)
    v_max = sqrt(driver.a_cp_max/abs(max_k))
    v_des_old = driver.mlon.v_des
    set_desired_speed!(driver.mlon, min(v_des_old,v_max))
    track_longitudinal!(driver.mlon, scene, roadway, vehicle_index, fore)
    set_desired_speed!(driver.mlon, v_des_old)

    driver.latAcc = rand(driver.mlat)
    driver.longAcc = rand(driver.mlon).a
    driver
end

function Base.rand(driver::UrbanDriver)
    LatLonAccelDirection(driver.latAcc, driver.longAcc, driver.direction)
end

Distributions.pdf(driver::UrbanDriver, a::LatLonAccelDirection) = pdf(driver.mlat, a.a_lat) * pdf(driver.mlon, a.a_lon)
Distributions.logpdf(driver::UrbanDriver, a::LatLonAccelDirection) = logpdf(driver.mlat, a.a_lat) * logpdf(driver.mlon, a.a_lon)


text:
# Author: Sean Wu
# Last Updated: December 11, 2020

function insensor(xy::Array{Float64,1})::Bool
    x, y = xy
    withinx = abs(x) <= SENSOR_HALF_WIDTH
    withiny = y <= 0
    return withinx && withiny
end

function inwalls(x::Float64, y::Float64)::Bool
    withinx = SENSOR_HALF_WIDTH < abs(x) <= ESCAPE_X
    withiny = y <= 0
    return withinx && withiny
end

function inenz(xy::Array{Float64,1})::Bool
    withinx = ENZYME_LEFT_X <= xy[1] <= ENZYME_RIGHT_X
    withiny = WALL_Y < xy[2] <= ENZYME_MAX_Y
    return withinx && withiny
end

function inwater(xy::Array{Float64,1})::Bool
    x, y = xy
    aboveEnzX = abs(x) <= ESCAPE_X
    aboveEnzY = ENZYME_MAX_Y < y <= ESCAPE_Y
    besideEnzX = ENZYME_RIGHT_X < abs(x) <= ESCAPE_X
    besideEnzY = WALL_Y < y <= ENZYME_MAX_Y
    return (aboveEnzX && aboveEnzY) || (besideEnzX && besideEnzY)
end

function inppd(xy::Array{Float64,1})::Bool
    return false
end

function sensorcases(x::Float64)::String
    return "sensor"
end

function wallcases(initxy::Array{Float64,1},
                   dx::Float64,
                   dy::Float64,
                   ending_step_size::Float64)
    return approach_wall!(initxy, dx, dy, ending_step_size), "no collision"
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
function parse_markdown(document_body; is_pandoc = false)
    header_text, document_body, offset = separate_header_text(document_body)
    header = parse_header(header_text)
    code_start, code_end = if is_pandoc
        r"^<<(?<options>.*?)>>=\s*$",
        r"^@\s*$"
    else
        r"^[`~]{3}(\{?)julia\s*([;,]?)\s*(?<options>.*?)(\}|\s*)$",
        r"^[`~]{3}\s*$"
    end
    return header, parse_markdown_body(document_body, code_start, code_end, offset)
end

# headers
# -------

const HEADER_REGEX = r"^---$(?<header>((?!---).)+)^---$"ms

# TODO: non-Weave headers should keep live in a doc
# separates header section from `text`
function separate_header_text(text)
    m = match(HEADER_REGEX, text)
    isnothing(m) && return "", text, 0
    header_text = m[:header]
    offset = @static if VERSION ≥ v"1.4"
        count("\n", header_text)
    else
        count(c->c==='\n', header_text)
    end
    return header_text, replace(text, HEADER_REGEX => ""; count = 1), offset
end

# HACK:
# YAML.jl can't parse text including ``` characters, so first replace all the inline code
# with these temporary code start/end string
const HEADER_INLINE_START = "<weave_header_inline_start>"
const   HEADER_INLINE_END = "<weave_header_inline_end>"

function parse_header(header_text)
    isempty(header_text) && return Dict()
    pat = INLINE_REGEX => SubstitutionString("$(HEADER_INLINE_START)\\1$(HEADER_INLINE_END)")
    header_text = replace(header_text, pat)
    return YAML.load(header_text)
end

# body
# ----

function parse_markdown_body(document_body, code_start, code_end, offset)
    lines = split(document_body, '\n')

    state = :doc
    doc_no = 0
    code_no = 0
    content = ""
    start_line = offset

    options = OptionDict()
    option_string = ""

    chunks = WeaveChunk[]
    for (line_no, line) in enumerate(lines)
        m = match(code_start, line)
        if !isnothing(m) && state === :doc
            state = :code

            option_string = isnothing(m[:options]) ? "" : strip(m[:options])
            options = parse_options(option_string)
            haskey(options, :label) && (options[:name] = options[:label])
            haskey(options, :name) || (options[:name] = nothing)

            isempty(strip(content)) || push!(chunks, DocChunk(content, doc_no += 1, start_line))

            start_line = line_no + offset
            content = ""
            continue
        end

        if occursin(code_end, line) && state === :code
            push!(chunks, CodeChunk(content, code_no += 1, start_line, option_string, options))

            start_line = line_no + offset
            content = ""
            state = :doc
            continue
        end

        content *= isone(line_no) ? line : string('\n', line)
    end

    # Remember the last chunk
    isempty(strip(content)) || push!(chunks, DocChunk(content, doc_no += 1, start_line))

    return chunks
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

@testset "functions" begin
    shape = NamedTupleShape(
        a = ArrayShape{Real}(3,2),
        b = ArrayShape{Real}(2)
    )

    x_flat = rand(@inferred totalndof(shape))
    x = @inferred stripscalar(shape(x_flat))

    f(x) = sum(x.a * x.b)

    @test @inferred(shape >> f) isa ValueShapes.FuncWithVarShape
    @test @inferred(varshape((shape >> f))) == shape
    @test @inferred(vardof((shape >> f))) == totalndof(shape)

    fws = shape >> f
    @test f(x) == @inferred(fws(x))
    @test f(x) == @inferred(fws(x_flat))
end


text:
function clipdb(s,cutoff)
  # [clipped] = clipdb(s,cutoff)
  # Clip magnitude of s at its maximum + cutoff in dB.
  # Example: clip(s,-100) makes sure the minimum magnitude
  # of s is not more than 100dB below its maximum magnitude.
  # If s is zero, nothing is done.

  as = abs(s);
  mas = maximum(as);
  if mas==0; return; end
  if cutoff >= 0; return; end

  thresh = mas*10^(cutoff/20); # db to linear
  clipped = float(s);
  clipped[ as .< thresh ] = thresh;
  return clipped;
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo

# ## Initial conditions
q1 = nominal_configuration(halfcheetah)
v1 = zeros(halfcheetah.nq)

# ## Time
h = 0.01
T = 100

# ## Simulator
s = Simulator(halfcheetah, T, h=h)

# ## Simulate
simulate!(s, q1, v1)

# ## Visualizer
vis = Visualizer()
render(vis)

# ## Visualize
visualize!(vis, s)