text:
"""
    DataSetExact(x::AbstractArray, y::AbstractArray, Σ_y::AbstractArray)
    DataSetExact(x::AbstractArray, Σ_x::AbstractArray, y::AbstractArray, Σ_y::AbstractArray)
    DataSetExact(xd::Distribution, yd::Distribution, dims::Tuple{Int,Int,Int}=(length(xd),1,1))
A data container which allows for uncertainties in the independent variables, i.e. ``x``-variables.
Moreover, the observed data is stored in terms of two probability distributions over the spaces ``\\mathcal{X}^N`` and ``\\mathcal{Y}^N`` respectively, which also allows for uncertainties in the observations that are non-Gaussian.
For instance, the uncertainties associated with a given observation might follow a Cauchy, t-student, log-normal or some other smooth distribution.

Examples:
```julia
using InformationGeometry, Distributions
X = product_distribution([Normal(0, 1), Cauchy(2, 0.5)])
Y = MvTDist(2, [3, 8.], [1 0.5; 0.5 3])
DataSetExact(X, Y, (2,1,1))
```

!!! note
    Uncertainties in the independent ``x``-variables are optional for `DataSetExact`, and can be set to zero by wrapping the `x`-data in a `InformationGeometry.Dirac` "distribution".
    The following illustrates numerically equivalent ways of encoding a dataset whose uncertainties in the ``x``-variables is zero:
    ```julia
    using InformationGeometry, Distributions, LinearAlgebra
    DS1 = DataSetExact(InformationGeometry.Dirac([1,2]), MvNormal([5,6], Diagonal([0.1, 0.2].^2)))
    DS2 = DataSetExact([1,2], [5,6], [0.1, 0.2])
    DS3 = DataSet([1,2], [5,6], [0.1, 0.2])
    ```
    where `DS1 == DS2 == DS3` will evaluate to `true`.
"""
struct DataSetExact <: AbstractDataSet
    xdist::Distribution
    ydist::Distribution
    dims::Tuple{Int,Int,Int}
    InvCov::AbstractMatrix{<:Number}
    WoundX::Union{AbstractVector,Nothing}
    xnames::AbstractVector{String}
    ynames::AbstractVector{String}
    DataSetExact(DM::AbstractDataModel) = DataSetExact(Data(DM))
    DataSetExact(DS::DataSet) = InformNames(DataSetExact(xDataDist(DS), yDataDist(DS), dims(DS)), xnames(DS), ynames(DS))
    DataSetExact(x::AbstractArray, y::AbstractArray, allsigmas::Real=1.0) = DataSetExact(x, y, allsigmas*ones(length(y)*length(y[1])))
    DataSetExact(x::AbstractArray, allxsigmas::Real=1.0, args...) = DataSetExact(x, allxsigmas*ones(length(x)*length(x[1])), args...)
    DataSetExact(x::AbstractArray, y::AbstractArray, yerr::AbstractArray) = DataSetExact(x, zeros(size(x,1)*length(x[1])), y, yerr)
    DataSetExact(x::AbstractVector{<:Number},y::AbstractVector{<:Measurement}) = DataSetExact(x,[y[i].val for i in 1:length(y)],[y[i].err for i in 1:length(y)])
    function DataSetExact(x::AbstractVector{<:Measurement}, y::AbstractVector{<:Measurement})
        DataSetExact([x[i].val for i in 1:length(x)],[x[i].err for i in 1:length(x)],[y[i].val for i in 1:length(y)],[y[i].err for i in 1:length(y)])
    end
    ###### No Unwinding above here.
    # Offload most of the checking to DataSet
    function DataSetExact(x::AbstractArray, Σ_x::AbstractArray, y::AbstractArray, Σ_y::AbstractArray, Dims::Union{Tuple{Int,Int,Int},Nothing}=nothing)
        DS = Dims isa Nothing ? DataSet(x, y, Σ_y) : DataSet(x, y, Σ_y, Dims)
        Dims = Dims isa Nothing ? dims(DS) : Dims
        dims(DS) != Dims && throw("DataSetExact: Given dims Tuple inconsistent: $Dims.")
        DataSetExact(DS, Σ_x)
    end
    function DataSetExact(DS::DataSet, Σ_x::AbstractArray)
        Σ_x = size(Σ_x,1) != size(Σ_x,2) ? Unwind(Σ_x) : Σ_x
        if (Σ_x == zeros(size(Σ_x,1))) || (Σ_x == Diagonal(zeros(size(Σ_x, 1))))
            return InformNames(DataSetExact(InformationGeometry.Dirac(xdata(DS)), yDataDist(DS), dims(DS)), xnames(DS), ynames(DS))
        else
            return InformNames(DataSetExact(DataDist(xdata(DS),HealthyCovariance(Σ_x)), yDataDist(DS), dims(DS)), xnames(DS), ynames(DS))
        end
    end
    function DataSetExact(xd::Distribution, yd::Distribution)
        println("No information about dimensionality of x-values or y-values given. Assuming that each x and y value has a single component from here on out.")
        DataSetExact(xd, yd, (length(xd),1,1))
    end
    function DataSetExact(xd::Distribution, yd::Distribution, dims::Tuple{Int,Int,Int})
        !(Int(length(xd)/xdim(dims)) == Int(length(yd)/ydim(dims)) == Npoints(dims)) && throw("Dimensions of given distributions are inconsistent with dimensions $dims.")
        Σinv = HealthyCovariance(InvCov(yd))
        if xdim(dims) == 1
            return DataSetExact(xd, yd, dims, Σinv, nothing)
        else
            # return new(xd,yd,dims,InvCov(yd),collect(Iterators.partition(GetMean(xd),xdim(dims))))
            return DataSetExact(xd, yd, dims, Σinv, [SVector{xdim(dims)}(Z) for Z in Windup(GetMean(xd),xdim(dims))])
        end
    end
    function DataSetExact(xd::Distribution, yd::Distribution, dims::Tuple{Int,Int,Int}, InvCov::AbstractMatrix{<:Number}, WoundX::Union{AbstractVector,Nothing},
                    xnames::AbstractVector{String}=CreateSymbolNames(xdim(dims),"x"), ynames::AbstractVector{String}=CreateSymbolNames(xdim(dims),"y"))
        new(xd, yd, dims, InvCov, WoundX, xnames, ynames)
    end
end

# For SciMLBase.remake
DataSetExact(;
xdist::Distribution=Normal(0,1),
ydist::Distribution=Normal(0,1),
dims::Tuple{Int,Int,Int}=(1,1,1),
InvCov::AbstractMatrix{<:Number}=Diagonal([1.]),
WoundX::Union{AbstractVector,Nothing}=nothing,
xnames::AbstractVector{String}=["x"],
ynames::AbstractVector{String}=["y"]) = DataSetExact(xdist, ydist, dims, InvCov, WoundX, xnames, ynames)

dims(DSE::DataSetExact) = DSE.dims
yInvCov(DSE::DataSetExact) = DSE.InvCov
xInvCov(DSE::DataSetExact) = InvCov(xdist(DSE))

WoundX(DS::DataSetExact) = _WoundX(DS, DS.WoundX)

xdist(DSE::DataSetExact) = DSE.xdist
ydist(DSE::DataSetExact) = DSE.ydist

GetMean(P::Product) = [location(P.v[i]) for i in 1:length(P)]
GetMean(P::Distribution) = mean(P)
# GetMean(P::Distribution) = P.μ

xdata(DSE::DataSetExact) = GetMean(xdist(DSE))
ydata(DSE::DataSetExact) = GetMean(ydist(DSE))

Sigma(P::Product) = [P.v[i].σ^2 for i in 1:length(P)] |> Diagonal
# I thought this was faster but it occasionally causes type problems with PDDiagMat etc.
# Sigma(P::Distribution) = P.Σ
Sigma(P::Distribution) = cov(P)
# Sigma(P::Distribution) = try P.Σ catch; cov(P) end
xsigma(DSE::DataSetExact) = Sigma(xdist(DSE))
ysigma(DSE::DataSetExact) = Sigma(ydist(DSE))

xnames(DSE::DataSetExact) = DSE.xnames
ynames(DSE::DataSetExact) = DSE.ynames

# function InformNames(DS::DataSetExact, xnames::AbstractVector{String}, ynames::AbstractVector{String})
#     @assert length(xnames) == xdim(DS) && length(ynames) == ydim(DS)
#     DataSetExact(xdist(DS), ydist(DS), (Npoints(DS),xdim(DS),ydim(DS)), yInvCov(DS), WoundX(DS), xnames, ynames)
# end

InvCov(P::Product) = [P.v[i].σ^(-2) for i in 1:length(P)] |> Diagonal
function InvCov(P::Distributions.GenericMvTDist)
    if P.df < 3
        return inv(P.Σ).mat
    else
        return invcov(P)
    end
end
InvCov(P::Distribution) = invcov(P)

DataMetric(P::Distribution) = InvCov(P)
function DataMetric(P::Distributions.GenericMvTDist)
    if P.df == 1
        return 0.5 .* InvCov(P)
    else
        println("DataMetric: Don't know what to do for t-distribution with dof=$(P.df), just returning usual inverse covariance matrix.")
        return InvCov(P)
    end
end

# Needs testing!!!!
isCauchy(P::Distribution{Univariate,Continuous}) = false;    isCauchy(P::Cauchy) = true
function DataMetric(P::Product)
    icov = InvCov(P).diag
    [isCauchy(P.v[i]) ? 0.5 : 1. for i in 1:length(P)] .* icov |> Diagonal
end

LogLike(DM::AbstractDataSet, x::AbstractVector{<:Number}, y::AbstractVector{<:Number}) = LogLike(Data(DM), x, y)
LogLike(DSE::DataSetExact, x::AbstractVector{<:Number}, y::AbstractVector{<:Number}) = logpdf(xdist(DSE),x) + logpdf(ydist(DSE),y)

import Distributions: loglikelihood
_loglikelihood(DSE::DataSetExact, model::ModelOrFunction, θ::AbstractVector{<:Number}; kwargs...) = LogLike(DSE, xdata(DSE), EmbeddingMap(DSE,model,θ; kwargs...))

function _Score(DSE::DataSetExact, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, Auto::Val{false}; kwargs...)
    transpose(EmbeddingMatrix(DSE,dmodel,θ; kwargs...)) * gradlogpdf(ydist(DSE), EmbeddingMap(DSE,model,θ; kwargs...))
end

# _FisherMetric(DS::DataSetExact, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}; kwargs...) = Pullback(DS,dmodel,DataMetric(DS),θ; kwargs...)


text:
module DBDSQLite
    using DBI, DataArrays, DataFrames

    export SQLite3

    include("types.jl")
    include("consts.jl")

    include(joinpath("api", "connect.jl"))
    include(joinpath("api", "disconnect.jl"))
    include(joinpath("api", "errors.jl"))
    include(joinpath("api", "execute.jl"))
    include(joinpath("api", "fetch.jl"))
    include(joinpath("api", "finish.jl"))
    include(joinpath("api", "metadata.jl"))
    include(joinpath("api", "misc.jl"))
    include(joinpath("api", "prepare.jl"))

    include("dbi.jl")
end

ext: jl


text:
#
# LazyMaps.jl -
#
# Implementation of lazy maps for Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LazyMaps.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017, Éric Thiébaut.
#

isdefined(Base, :__precompile__) && __precompile__(true)

module LazyMaps

export
    LazyMap,
    lazymap

# Singleton type to indicate that the element type is unknown.
struct Unknown; end
const unknown = Unknown()

struct LazyMapArray{F,A<:AbstractArray,T,N} <: AbstractArray{T,N}
    f::F # function or callable
    a::A # input array/collection
    v::T # first value or `unknown`
end

struct LazyMapOther{F,A,T,N}
    f::F # function or callable
    a::A # input array/collection
end

const LazyMap{F,A,T,N} = Union{LazyMapArray{F,A,T,N},
                               LazyMapOther{F,A,T,N}}

"""
    lazymap(f, a [, guess=true])

yields a read-only collection whose items are given by the function `f`
applied *on-the-fly* to the items of `a`.  This is useful to provide an
elementwise filtered array or collection to some method without creating an
intermediate array or collection.

If `guess` is `true`, the method attempts to evaluate `f` for the first
item of `a` to determine the type of the elements.  Set `guess` to `false`
to avoid this at the cost of making some optimizations not possible.  Note
that if the element type can be determined, it is assumed (and asserted)
that all elements have the same type, *i.e.* the elements of result are
**type stable**.

See also [`map`](@ref).

"""
lazymap(f, a::AbstractArray) = _lazymap(f, a, f(first(a)))
lazymap(f, a::AbstractArray, guess::Bool) =
    (guess ? lazymap(f, a) : _lazymap(f, a, unknown))
lazymap(f::F, a::A) where {F,A} = LazyMapOther{F,A,Unknown,0}(f, a)
lazymap(f::F, a::A, ::Bool) where {F,A} = LazyMapOther{F,A,Unknown,0}(f, a)

_lazymap(f::F, a::A, v::T) where {F,A<:AbstractArray{E,N},T} where {E,N} =
    LazyMapArray{F,A,T,N}(f, a, v)

Base.getindex(m::LazyMapArray{F,A,T,N}, key...) where {F,A,T,N} =
    m.f(m.a[key...]) :: T
Base.getindex(m::LazyMap{F,A,Unknown,N}, key...) where {F,A,N} =
    m.f(m.a[key...])

Base.eltype(m::LazyMap{F,A,T,N}) where {F,A,T,N} = T
Base.eltype(m::LazyMap{F,A,Unknown,N}) where {F,A,N} =
    error("element type is unknown")
Base.length(m::LazyMap) = length(m.a)
Base.ndims(m::LazyMapArray{F,A,T,N}) where {F,A,T,N} = N
Base.size(m::LazyMapArray) = size(m.a)
Base.size(m::LazyMapArray, d) = size(m.a, d)
Base.stride(m::LazyMapArray, k) = stride(m.a, k)
Base.strides(m::LazyMapArray) = strides(m.a)
Base.indices(m::LazyMapArray) = indices(m.a)
Base.indices(m::LazyMapArray, d) = indices(m.a, d)
Base.first(m::LazyMap{F,A,Unknown,N}) where {F,A,N} = m.f(first(m.a))
Base.first(m::LazyMapArray) = m.v
Base.last(m::LazyMap) = m.f(last(m.a))
Base.endof(m::LazyMap) = endof(m.a)
#Base.eachindex(m::LazyMap) = eachindex(m.a)
#Base.eachindex(m::LazyMap, args...) = eachindex(m.a, args...)

# Make an instance of a LazyMap an iterable.
Base.start(m::LazyMap) = start(m.a)
Base.done(m::LazyMap, state) = done(m.a, state)
function Base.next(m::LazyMap, state)
    item, state = next(m.a, state)
    return m.f(item), state
end

# Provide traits for LazyMap instances.
Base.IndexStyle(::Type{LazyMapArray{F,A,T,N}}) where {F,A,T,N} =
    Base.IndexStyle(A)
Base.iteratorsize(::Type{LazyMap{F,A,T,N}}) where {F,A,T,N} =
    Base.iteratorsize(A)
Base.iteratoreltype(::Type{LazyMap{F,A,Unknown,N}}) where {F,A,N} =
    Base.EltypeUnknown()
Base.iteratoreltype(::Type{LazyMap{F,A,T,N}}) where {F,A,T,N} =
    Base.HasEltype()

end


text:
export monodromy_solve,
    MonodromyResult,
    real_solutions,
    is_success,
    is_heuristic_stop,
    nreal,
    parameters,
    verify_solution_completeness,
    solution_completeness_witnesses

#####################
# Monodromy Options #
#####################
const monodromy_options_supported_keywords = [
    :distance,
    :identical_tol,
    :done_callback,
    :group_action,
    :group_actions,
    :group_action_on_all_nodes,
    :parameter_sampler,
    :equivalence_classes,
    :complex_conjugation,
    :check_startsolutions,
    :target_solutions_count,
    :timeout,
    :min_solutions,
    :max_loops_no_progress,
    :reuse_loops,
]

struct MonodromyOptions{F<:Function,F1<:Function,F2<:Tuple,F3<:Function}
    distance_function::F
    identical_tol::Float64
    done_callback::F1
    group_actions::GroupActions{F2}
    group_action_on_all_nodes::Bool
    parameter_sampler::F3
    equivalence_classes::Bool
    complex_conjugation::Bool
    check_startsolutions::Bool
    # stopping heuristic
    target_solutions_count::Int
    timeout::Float64
    min_solutions::Int
    max_loops_no_progress::Int
    reuse_loops::Symbol
end

function MonodromyOptions(
    is_real_system::Bool,
    accuracy::Float64,
    distance;
    identical_tol::Float64 = sqrt(accuracy),
    done_callback = always_false,
    group_action = nothing,
    group_actions = group_action === nothing ? nothing : GroupActions(group_action),
    group_action_on_all_nodes = false,
    parameter_sampler = independent_normal,
    equivalence_classes = true,
    complex_conjugation = is_real_system,
    check_startsolutions = true,
    # stopping heuristic
    target_solutions_count = nothing,
    timeout = float(typemax(Int)),
    min_solutions::Int = default_min_solutions(target_solutions_count),
    max_loops_no_progress::Int = 10,
    reuse_loops::Symbol = :all,
)

    if group_actions isa GroupActions
        actions = group_actions
    else
        if group_actions === nothing
            equivalence_classes = false
        end
        actions = GroupActions(group_actions)
    end

    MonodromyOptions(
        distance,
        identical_tol,
        done_callback,
        actions,
        group_action_on_all_nodes,
        parameter_sampler,
        equivalence_classes,
        complex_conjugation,
        check_startsolutions,
        target_solutions_count === nothing ? typemax(Int) : target_solutions_count,
        float(timeout),
        min_solutions,
        max_loops_no_progress,
        reuse_loops,
    )
end

default_min_solutions(::Nothing) = 1
function default_min_solutions(target_solutions_count::Int)
    div(target_solutions_count, 2)
end

always_false(x, sols) = false

has_group_actions(options::MonodromyOptions) =
    !(options.group_actions isa GroupActions{Tuple{}})

"""
    independent_normal(p::AbstractVector{T}) where {T}

Sample a vector where each entries is drawn independently from the univariate normal distribution.
"""
independent_normal(p::SVector{N,T}) where {N,T} = @SVector randn(T, N)
independent_normal(p::AbstractVector{T}) where {T} = randn(T, length(p))

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

mutable struct MonodromyStatistics
    ntrackedpaths::Int
    ntrackingfailures::Int
    nreal::Int
    nparametergenerations::Int
    nsolutions_development::Vector{Int}
end

MonodromyStatistics(nsolutions::Int) = MonodromyStatistics(0, 0, 0, 1, [nsolutions])
function MonodromyStatistics(solutions)
    stats = MonodromyStatistics(length(solutions))
    for s in solutions
        if is_real_vector(s)
            stats.nreal += 1
        end
    end
    stats
end

Base.show(io::IO, S::MonodromyStatistics) = print_fieldnames(io, S)
Base.show(io::IO, ::MIME"application/prs.juno.inline", S::MonodromyStatistics) = S

# update routines
function trackedpath!(stats::MonodromyStatistics, retcode)
    if is_success(retcode)
        stats.ntrackedpaths += 1
    else
        stats.ntrackingfailures += 1
    end
end

function generated_parameters!(stats::MonodromyStatistics, nsolutions::Int)
    stats.nparametergenerations += 1
    push!(stats.nsolutions_development, nsolutions)
end

function finished!(stats, nsolutions)
    push!(stats.nsolutions_development, nsolutions)
end

function n_completed_loops_without_change(stats, nsolutions)
    k = 0
    for i = length(stats.nsolutions_development):-1:1
        if stats.nsolutions_development[i] != nsolutions
            return max(k - 1, 0)
        end
        k += 1
    end
    return max(k - 1, 0)
end

function n_solutions_current_loop(statistics, nsolutions)
    nsolutions - statistics.nsolutions_development[end]
end

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

export Triangle, Petal

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

struct LoopEdge
    p₁::Vector{ComplexF64}
    p₀::Vector{ComplexF64}
    weights::Union{Nothing,NTuple{2,ComplexF64}}
end

function LoopEdge(p₁, p₀; weights = false)
    if weights
        γ = (randn(ComplexF64), randn(ComplexF64))
    else
        γ = nothing
    end
    LoopEdge(convert(Vector{ComplexF64}, p₁), convert(Vector{ComplexF64}, p₀), γ)
end

struct MonodromyLoop
    edges::Vector{LoopEdge}
end

function MonodromyLoop(base_p, nnodes::Int, options::MonodromyOptions; weights = true)
    edges = LoopEdge[]
    p₁ = base_p
    for i = 2:nnodes
        p₀ = options.parameter_sampler(p₁)
        push!(edges, LoopEdge(p₁, p₀; weights = weights))
        p₁ = p₀
    end
    push!(edges, LoopEdge(p₁, base_p; weights = weights))
    MonodromyLoop(edges)
end

##############
# Loop styles
##############
"""
    LoopStyle

Abstract type defining a style of a loop.
"""
abstract type LoopStyle end

"""
    Triangle(;useweights=true)

A triangle is a loop consisting of the main node and two addtional nodes.
If `weights` is true the edges are equipped with additional random weights.
Note that this is usually only necessary for real parameters.
"""
struct Triangle <: LoopStyle
    useweights::Bool
end
Triangle(; useweights = true) = Triangle(useweights)

function MonodromyLoop(strategy::Triangle, p, options::MonodromyOptions)
    MonodromyLoop(p, 3, options, weights = strategy.useweights)
end

"""
    Petal()

A petal is a loop consisting of the main node and one other node connected
by two edges with different random weights.
"""
struct Petal <: LoopStyle end
function MonodromyLoop(strategy::Petal, p, options::MonodromyOptions)
    MonodromyLoop(p, 2, options, weights = true)
end

"""
    regenerate!(loop::MonodromyLoop, options::MonodromyOptions, stats::MonodromyStatistics)

Regenerate all random parameters in the loop in order to introduce a new monodromy action.
"""
function regenerate!(
    loop::MonodromyLoop,
    options::MonodromyOptions,
    stats::MonodromyStatistics,
)
    main = mainnode(loop)

    # The first node is the main node and doesn't get touched
    store_points = options.group_action_on_all_nodes && has_group_actions(options)
    for i ∈ 2:length(loop.nodes)
        loop.nodes[i] = Node(
            options.parameter_sampler(main.p),
            loop.nodes[i],
            options;
            store_points = store_points, is_main_node = false,
        )
    end
    loop.edges .= Edge.(loop.edges)
    generated_parameters!(stats, length(main.points)) # bookkeeping
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(
    tracker::CoreTracker,
    x::AbstractVector,
    loop::MonodromyLoop,
    ntrackedpaths::Base.RefValue{Int},
    ntrackingfailures::Base.RefValue{Int},
)
    H = basehomotopy(tracker.homotopy)
    local retcode::CoreTrackerStatus.states
    y = x
    for e in loop.edges
        set_parameters!(H, (e.p₁, e.p₀), e.weights)
        retcode = track!(tracker, y)
        if is_success(retcode)
            ntrackedpaths[] += 1
        else
            ntrackingfailures[] += 1
        end
        is_success(retcode) || break
        y = current_x(tracker)
    end
    retcode
end

#############
## Results ##
#############
"""
    MonodromyResult

The monodromy result contains the result of the `monodromy_solve` computation.
"""
struct MonodromyResult{AV<:AbstractVector,T}
    returncode::Symbol
    solutions::Vector{AV}
    parameters::Vector{T}
    statistics::MonodromyStatistics
    equivalence_classes::Bool
    seed::Int
end

Base.iterate(R::MonodromyResult) = iterate(R.solutions)
Base.iterate(R::MonodromyResult, state) = iterate(R.solutions, state)

Base.show(io::IO, ::MIME"application/prs.juno.inline", x::MonodromyResult) = x
function Base.show(io::IO, result::MonodromyResult{N,T}) where {N,T}
    println(io, "MonodromyResult")
    println(io, "==================================")
    if result.equivalence_classes
        println(
            io,
            "• ",
            nsolutions(result),
            " classes of solutions (modulo group action) (",
            nreal(result),
            ") real)",
        )
    else
        println(io, "• $(nsolutions(result)) solutions ($(nreal(result)) real)")
    end
    println(io, "• return code → $(result.returncode)")
    println(io, "• $(result.statistics.ntrackedpaths) tracked paths")
    println(io, "• seed → $(result.seed)")
end

TreeViews.hastreeview(::MonodromyResult) = true
TreeViews.numberofnodes(::MonodromyResult) = 6
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)) classes of solutions (modulo group action)")
        else
            print(io, "$(nsolutions(x)) solutions")
        end
    elseif i == 2
        if x.equivalence_classes
            print(io, "$(nreal(x)) classes of real solutions")
        else
            print(io, "$(nreal(x)) real solutions")
        end
    elseif i == 3
        print(io, "Return code")
    elseif i == 4
        print(io, "Statistics")
    elseif i == 5
        print(io, "Parameters")
    elseif i == 6
        print(io, "Seed")
    end
end
function TreeViews.treenode(r::MonodromyResult, i::Integer)
    if i == 1
        return r.solutions
    elseif i == 2
        return real_solutions(r)
    elseif i == 3
        return r.returncode
    elseif i == 4
        return r.statistics
    elseif i == 5
        return r.parameters
    elseif i == 6
        return r.seed
    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

"""
    mapresults(f, result::MonodromyResult; only_real=false, real_tol=1e-6)

Apply the function `f` to all entries of `MonodromyResult` for which the given conditions apply.

## Example
```julia
# This gives us all solutions considered real (but still as a complex vector).
real_solutions = mapresults(solution, R, only_real=true)
```
"""
function mapresults(f, R::MonodromyResult; only_real = false, real_tol = 1e-6)
    [f(r) for r in R.solutions if (!only_real || is_real_vector(r, real_tol))]
end

"""
    solutions(result::MonodromyResult; only_real=false, real_tol=1e-6)

Return all solutions (as `SVector`s) for which the given conditions apply.

## Example
```julia
real_solutions = solutions(R, only_real=true)
```
"""
function solutions(R::MonodromyResult; kwargs...)
    mapresults(identity, R; kwargs...)
end

"""
    nsolutions(result::MonodromyResult)

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

"""
    real_solutions(res::MonodromyResult; tol=1e-6)

Returns the solutions of `res` whose imaginary part has norm less than 1e-6.
"""
function real_solutions(res::MonodromyResult; tol = 1e-6)
    map(r -> real_vector(r), filter(r -> is_real_vector(r, tol), res.solutions))
end

"""
    nreal(res::MonodromyResult; tol=1e-6)

Counts how many solutions of `res` have imaginary part norm less than 1e-6.
"""
function nreal(res::MonodromyResult; tol = 1e-6)
    count(r -> is_real_vector(r, tol), res.solutions)
end

"""
    parameters(r::MonodromyResult)

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

#####################
## monodromy solve ##
#####################
struct MonodromySolver{
    T<:Number,
    UP<:UniquePoints,
    MO<:MonodromyOptions,
    Tracker<:CoreTracker,
    Prob<:Problem,
}
    parameters::Vector{T}
    solutions::UP
    solutions_lock::ReentrantLock
    loops::Vector{MonodromyLoop}
    options::MO
    statistics::MonodromyStatistics
    trackers::Vector{Tracker}
    problem::Prob
end

function MonodromySolver(
    F,
    solution::AbstractVector{<:Number},
    p₀::AbstractVector{TP};
    kwargs...,
) where {TP}
    MonodromySolver(F, [solution], p₀; kwargs...)
end
function MonodromySolver(
    F,
    startsolutions::Vector{<:AbstractVector{<:Number}},
    p::AbstractVector{TP};
    distance = default_distance(eltype(startsolutions)),
    parameters = nothing,
    strategy = nothing,
    show_progress = true,
    showprogress = nothing, #deprecated
    accuracy = 1e-6, # set a higher refinement_accuracy than the default
    precision = PRECISION_ADAPTIVE,
    kwargs...,
) where {TP,NVars}

    @deprecatekwarg showprogress show_progress

    if parameters !== nothing && length(p) ≠ length(parameters)
        throw(ArgumentError("Number of provided parameters doesn't match the length of initially provided parameter `p₀`."))
    end

    x₀ = complex.(first(startsolutions))
    p₀ = Vector{promote_type(Float64, TP)}(p)

    optionskwargs, restkwargs = splitkwargs(kwargs, monodromy_options_supported_keywords)

    supported, rest = splitkwargs(restkwargs, problem_startsolutions_supported_keywords)
    prob, _ = problem_startsolutions(
        F,
        [x₀];
        parameters = parameters, generic_parameters = p₀, γ₁ = 1, γ₀ = 1, supported...,
    )
    tracker =
        CoreTracker(prob, x₀; accuracy = accuracy, precision = PRECISION_ADAPTIVE, rest...)
    # Check whether homotopy is real
    is_real_system = numerically_check_real(tracker.homotopy, x₀)
    options = MonodromyOptions(is_real_system, accuracy, distance; optionskwargs...)
    # construct UniquePoints

    uniquepoints = UniquePoints(
        typeof(x₀),
        distance;
        group_actions = options.equivalence_classes ? options.group_actions : nothing,
        check_real = true,
    )

    # add only unique points that are true solutions
    for s in startsolutions
        if !options.check_startsolutions || is_valid_start_value(tracker, s, 0.0)
            add!(uniquepoints, s; tol = options.identical_tol)
        end
    end

    statistics = MonodromyStatistics(uniquepoints)

    if strategy === nothing
        strategy = default_strategy(F, parameters, p; is_real_system = is_real_system)
    end

    # construct Loop
    loops = [MonodromyLoop(strategy, p₀, options)]

    MonodromySolver(
        p₀,
        uniquepoints,
        ReentrantLock(),
        loops,
        options,
        statistics,
        [tracker],
        prob,
    )
end

"""
    solutions(MS::MonodromySolver)

Get the solutions of the loop.
"""
solutions(MS::MonodromySolver) = MS.solutions.points

"""
    nsolutions(loop::MonodromyLoop)

Get the number solutions of the loop.
"""
nsolutions(MS::MonodromySolver) = length(solutions(MS))

"""
    monodromy_solve(F, [sols, p]; parameters=..., options..., pathtrackerkwargs...)

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

## Options

* `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.
* `distance_function=euclidean_distance`: The distance function used for
  [`UniquePoints`](@ref).
* `done_callback=always_false`: A callback to end the computation early. This function takes
  2 arguments. The first one is the new solution `x` and the second one are all current
  solutions (including `x`). Return `true` if the compuation is done.
* `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.
* `identical_tol::Float64=1e-6`: The tolerance with which it is decided whether two
  solutions are identical.
* `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::Int=10`: 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 this is half of `target_solutions_count` if applicable otherwise 1.
* `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 the univariate normal distribution.
* `resuse_loops::Symbol=: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.
* `show_progress=true`: Enable a progress meter.
* `strategy`: The strategy used to create loops. If `F` only depends linearly on `p` this
  will be [`Petal`](@ref). Otherwise this will be [`Triangle`](@ref) with weights if `F`
  is a real system.
* `target_solutions_count=nothing`: The computations are stopped if this number of
  solutions is reached.
* `threading = true`: Enable multithreading of the path tracking.
* `timeout=float(typemax(Int))`: The maximal number of *seconds* the computation is allowed
  to run.
"""
function monodromy_solve(
    F::Union{MPPolys,System};
    parameters = F isa MPPolys ?
                 throw(ArgumentError("Necessary to set `parameters = ...`.")) : nothing,
    variable_ordering = nothing,
    seed = randseed(),
    show_progress = true,
    threading::Bool = true,
    kwargs...,
)
    Random.seed!(seed)
    x₀ = p₀ = nothing
    try
        x₀, p₀ = find_start_pair(F, parameters; variable_ordering = variable_ordering)
    catch
        throw(ArgumentError("Cannot compute a start pair (x, p) automatically."))
    end
    monodromy_solve!(
        MonodromySolver(
            F,
            [x₀],
            p₀;
            parameters = parameters,
            variable_ordering = variable_ordering,
            kwargs...,
        ),
        seed;
        show_progress = show_progress,
        threading = threading,
    )
end

function monodromy_solve(
    F::Vector{<:ModelKit.Expression};
    parameters = throw(ArgumentError("Necessary to set `parameters = ...`.")),
    variable_ordering = throw(ArgumentError("Necessary to set `variable_ordering = ...`.")),
    kwargs...,
)
    monodromy_solve(System(F, variable_ordering, parameters))
end
function monodromy_solve(
    F::Vector{<:ModelKit.Expression},
    S,
    p;
    parameters = throw(ArgumentError("Necessary to set `parameters = ...`.")),
    variable_ordering = throw(ArgumentError("Necessary to set `variable_ordering = ...`.")),
    kwargs...,
)
    monodromy_solve(System(F, variable_ordering, parameters), S, p)
end

function find_start_pair(F::System, ::Nothing; variable_ordering = nothing)
    vars = polyvar.(F.variables)
    params = polyvar.(F.parameters)
    P = evaluate(F.expressions, F.variables => vars, F.parameters => params)
    find_start_pair(P, params; variable_ordering = vars)
end

function find_start_pair(F::MPPolys, parameters; variable_ordering = nothing)
    if variable_ordering !== nothing
        vars = variable_ordering
    else
        vars = variables(F, parameters)
    end
    x₀ = randn(ComplexF64, length(vars))
    if F isa Composition
        F = expand(F)
    end
    F_x₀ = [subs(f, vars => x₀) for f in F]
    # let's use the helper function above to get the matrix form
    A, b = linear_system(F_x₀)
    if iszero(b)
        p_end = randn(ComplexF64)
        p₀ = [LA.qr(A[:, 1:end-1], Val(true)) \ ((-p_end) .* A[:, end]); p_end]
    else
        p₀ = LA.qr(A, Val(true)) \ b
    end

    all(
        f -> abs(f(parameters => p₀)) < 1e-8,
        F_x₀,
    ) || throw(ArgumentError("Cannot compute start pair."))

    (x₀, p₀)
end

function monodromy_solve(
    args...;
    seed = randseed(),
    show_progress = true,
    threading::Bool = true,
    kwargs...,
)
    Random.seed!(seed)
    monodromy_solve!(
        MonodromySolver(args...; kwargs...),
        seed;
        show_progress = show_progress,
        threading = threading,
    )
end

function default_strategy(
    F::Union{System,MPPolyInputs},
    parameters,
    p₀::AbstractVector{TP};
    is_real_system = false,
) where {TC,TP}
    # If F depends only linearly on the parameters a petal is sufficient
    vars = variables(F, parameters)
    if all(d -> d ≤ 1, maxdegrees(F; parameters = vars))
        Petal()
    # For a real system we should introduce some weights to avoid the discriminant
    elseif is_real_system
        Triangle(useweights = true)
    else
        Triangle(useweights = false)
    end
end
function default_strategy(F::AbstractSystem, parameters, p₀; is_real_system = false)
    # For a real system we should introduce some weights to avoid the discriminant
    if is_real_system
        Triangle(useweights = true)
    else
        Triangle(useweights = false)
    end
end

function numerically_check_real(H::AbstractHomotopy, x)
    y = copy(x)
    Random.randn!(y)
    for i in eachindex(y)
        y[i] = real(y[i]) + 0.0im
    end
    t = rand()
    r = evaluate(H, y, t)
    isapprox(LinearAlgebra.norm(imag.(r), Inf), 0.0, atol = 1e-14)
end

#################
## Actual work ##
#################

"""
    MonodromyJob

A `MonodromyJob` is consisting of a solution id and a loop id.
"""
struct MonodromyJob
    id::Int
    loop_id::Int
end

function monodromy_solve!(
    MS::MonodromySolver,
    seed;
    show_progress::Bool = true,
    threading::Bool = true,
)
    if nsolutions(MS) == 0
        @warn "None of the provided solutions is a valid start solution."
        return MonodromyResult(
            :invalid_startvalue,
            similar(MS.solutions.points, 0),
            MS.parameters,
            MS.statistics,
            MS.options.equivalence_classes,
            seed,
        )
    end
    # solve
    retcode = :not_assigned
    if show_progress
        if !MS.options.equivalence_classes
            desc = "Solutions found:"
        else
            desc = "Classes of solutions (modulo group action) found:"
        end
        progress =
            ProgressMeter.ProgressUnknown(desc; delay = 0.5, clear_output_ijulia = true)
    else
        progress = nothing
    end

    n_blas_threads = single_thread_blas()

    retcode = monodromy_solve!(MS, seed, progress; threading = threading)

    n_blas_threads > 1 && set_num_BLAS_threads(n_blas_threads)
    finished!(MS.statistics, nsolutions(MS))
    MonodromyResult(
        retcode,
        solutions(MS),
        MS.parameters,
        MS.statistics,
        MS.options.equivalence_classes,
        seed,
    )
end

function monodromy_solve!(
    MS::MonodromySolver,
    seed,
    progress::Union{Nothing,ProgressMeter.ProgressUnknown};
    threading::Bool = true,
)
    Threads.resize_nthreads!(MS.trackers)
    t₀ = time_ns()
    iterations_without_progress = 0 # stopping heuristic
    # intialize job queue
    queue = MonodromyJob.(1:nsolutions(MS), 1)
    thread_queues = [MonodromyJob[] for _ = 1:Threads.nthreads()]

    n = nsolutions(MS)
    retcode = :none
    while n < MS.options.target_solutions_count
        retcode =
            empty_queue!(queue, thread_queues, MS, t₀, progress; threading = threading)

        if retcode == :done
            retcode = :success
            break
        elseif retcode == :timeout ||
               retcode == :invalid_startvalue || retcode == :interrupt
            break
        end

        for q in thread_queues
            append!(queue, q)
            empty!(q)
        end
        isempty(queue) || continue

        # Iterations heuristic
        n_new = nsolutions(MS)
        if n == n_new
            iterations_without_progress += 1
        else
            iterations_without_progress = 0
            n = n_new
        end
        if iterations_without_progress == MS.options.max_loops_no_progress &&
           n_new ≥ MS.options.min_solutions
            retcode = :heuristic_stop
            break
        end

        regenerate_loop_and_schedule_jobs!(queue, MS)
    end
    update_progress!(progress, nsolutions(MS), MS.statistics; finish = true)

    retcode
end

function empty_queue!(
    queue,
    thread_queues,
    MS::MonodromySolver,
    t₀::UInt,
    progress;
    threading::Bool = true,
)
    if !threading
        ntracked = Ref(0)
        nfailures = Ref(0)
        n_since_update = 0
        retcode = :incomplete
        try
            while !isempty(queue)
                job = pop!(queue)
                status = process!(queue, job, MS, progress, ntracked, nfailures, 1)
                if status == :done || status == :invalid_startvalue || status == :interrupt
                    retcode = status
                    break
                end
                n_since_update += 1

                MS.statistics.ntrackedpaths += ntracked[]
                MS.statistics.ntrackingfailures += nfailures[]
                ntracked[] = nfailures[] = 0
                if n_since_update ≥ 4
                    n_since_update = 0
                    update_progress!(progress, nsolutions(MS), MS.statistics)
                end
                # check timeout
                if (time_ns() - t₀) > MS.options.timeout * 1e9 # convert s to ns
                    retcode = :timeout
                    break
                end
            end
        catch e
            if isa(e, InterruptException)
                return :interrupt
            else
                rethrow()
            end
        end
        n_since_update > 0 && update_progress!(progress, nsolutions(MS), MS.statistics)
        retcode
    else
        stop = Ref(false)
        return_status = Ref(:incomplete)
        start_ntrackedpaths = MS.statistics.ntrackedpaths
        start_ntrackingfailures = MS.statistics.ntrackingfailures
        ntrackedpaths = [Ref(0) for _ = 1:Threads.nthreads()]
        ntrackingfailures = [Ref(0) for _ = 1:Threads.nthreads()]
        Threads.@threads for job in queue
            try
                !stop[] || continue
                tid = Threads.threadid()
                status = process!(
                    thread_queues[tid],
                    job,
                    MS,
                    progress,
                    ntrackedpaths[tid],
                    ntrackingfailures[tid],
                    tid,
                )
                if status == :done
                    return_status[] = :done
                    stop[] = true
                elseif status == :invalid_startvalue
                    return_status[] = :invalid_startvalue
                    stop[] = true
                end
                if tid == 1
                    n_tracked = start_ntrackedpaths + sum(getindex, ntrackedpaths)
                    n_failures = start_ntrackingfailures + sum(getindex, ntrackingfailures)
                    MS.statistics.ntrackedpaths = n_tracked
                    MS.statistics.ntrackingfailures = n_failures
                    update_progress!(progress, nsolutions(MS), MS.statistics)
                end
            # check timeout
                if (time_ns() - t₀) > MS.options.timeout * 1e9 # convert s to ns
                    return_status[] = :invalid_startvalue
                    stop[] = true
                end
                if stop[]
                    break
                end
            catch e
                if isa(e, InterruptException)
                    return_status[] = :interrupt
                    stop[] = true
                    break
                else
                    rethrow()
                end
            end
        end
        empty!(queue)
        n_tracked = start_ntrackedpaths + sum(getindex, ntrackedpaths)
        n_failures = start_ntrackingfailures + sum(getindex, ntrackingfailures)
        MS.statistics.ntrackedpaths = n_tracked
        MS.statistics.ntrackingfailures = n_failures
        update_progress!(progress, nsolutions(MS), MS.statistics)

        return_status[]
    end
end

function process!(
    queue,
    job::MonodromyJob,
    MS::MonodromySolver,
    progress,
    ntrackedpaths::Base.RefValue{Int},
    ntrackingfailures::Base.RefValue{Int},
    tid::Integer,
)
    x = solutions(MS)[job.id]
    loop = MS.loops[job.loop_id]
    retcode = track(MS.trackers[tid], x, loop, ntrackedpaths, ntrackingfailures)
    is_success(retcode) || return :incomplete

    y = pull_back(MS.problem, current_x(MS.trackers[tid]))
    add_and_schedule!(MS, queue, y, job) && return :done

    if MS.options.complex_conjugation
        add_and_schedule!(MS, queue, conj.(y), job) && return :done
    end

    if !MS.options.equivalence_classes
        apply_actions(MS.options.group_actions, y) do yᵢ
            add_and_schedule!(MS, queue, yᵢ, job)
        end
    end
    return :incomplete
end

"""
    add_and_schedule!(MS, queue, y, job)

Add `y` to the current `node` (if it not already exists) and schedule a new job to the
`queue`. Returns `true` if we are done. Otherwise `false`.
"""
function add_and_schedule!(MS::MonodromySolver, queue, y, job::MonodromyJob) where {N,T}
    lock(MS.solutions_lock)
    k = add!(MS.solutions, y, Val(true); tol = MS.options.identical_tol)
    unlock(MS.solutions_lock)
    if k == NOT_FOUND || k == NOT_FOUND_AND_REAL
        # Check if we are done
        isdone(MS.solutions, y, MS.options) && return true
        push!(queue, MonodromyJob(nsolutions(MS), job.loop_id))
        # Schedule also on other loops
        if MS.options.reuse_loops == :random && length(MS.loops) > 1
            r_loop_id = rand(2:length(MS.loops))
            if r_loop_id == job.loop_id
                r_loop_id -= 1
            end
            push!(queue, MonodromyJob(nsolutions(MS), r_loop_id))
        elseif MS.options.reuse_loops == :all
            for r_loop_id = 1:length(MS.loops)
                r_loop_id != job.loop_id || continue
                push!(queue, MonodromyJob(nsolutions(MS), r_loop_id))
            end
        end
    end
    MS.statistics.nreal += (k == NOT_FOUND_AND_REAL)
    false
end

function update_progress!(
    ::Nothing,
    nsolutions,
    statistics::MonodromyStatistics;
    finish = false,
)
    nothing
end
function update_progress!(
    progress,
    nsolutions,
    statistics::MonodromyStatistics;
    finish = false,
)
    ProgressMeter.update!(
        progress,
        nsolutions,
        showvalues = (
            ("# paths tracked", statistics.ntrackedpaths),
            ("# loops generated", statistics.nparametergenerations),
            (
                "# completed loops without change",
                n_completed_loops_without_change(statistics, nsolutions),
            ),
            (
                "# solutions in current loop",
                n_solutions_current_loop(statistics, nsolutions),
            ),
            ("# real solutions", statistics.nreal),
        ),
    )
    if finish
        ProgressMeter.finish!(progress)
    end
    nothing
end

function isdone(solutions::UniquePoints, x, options::MonodromyOptions)
    options.done_callback(
        x,
        solutions.points,
    ) || length(solutions) ≥ options.target_solutions_count
end

function regenerate_loop_and_schedule_jobs!(queue, MS::MonodromySolver)
    es = MS.loops[1].edges
    loop = MonodromyLoop(
        MS.parameters,
        length(es),
        MS.options,
        weights = !isnothing(es[1].weights),
    )
    push!(MS.loops, loop)
    for id = 1:nsolutions(MS)
        push!(queue, MonodromyJob(id, length(MS.loops)))
    end
    generated_parameters!(MS.statistics, nsolutions(MS))
    nothing
end

##################
## VERIFICATION ##
##################
"""
    verify_solution_completeness(F, res::MonodromyResult;
                                 parameters=..., trace_tol=1e-6, options...)

Verify that the monodromy computation found all solutions by [`monodromy_solve`](@ref).
This uses a trace test as described in [^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. The `options...` arguments can be everything which is accepted by `solve` and
`monodromy_solve`.

### Example

```
julia> @polyvar x y a b c;

julia> f = x^2+y^2-1;

julia> l = a*x+b*y+c;

julia> res = monodromy_solve([f,l], [-0.6-0.8im, -1.2+0.4im], [1,2,3]; parameters=[a,b,c])
MonodromyResult
==================================
• 2 solutions (0 real)
• return code → heuristic_stop
• 44 tracked paths
• seed → 367230

julia> verify_solution_completeness([f,l], res; parameters=[a,b,c], trace_tol = 1e-8)
[ Info: Compute additional witnesses for completeness...
[ Info: Found 2 additional witnesses
[ Info: Compute trace...
[ Info: Norm of trace: 1.035918995391323e-15
true
```

    verify_solution_completeness(F, S, p; parameters=..., kwargs...)

Verify the solution completeness using the computed solutions `S` to the parameter `p`.

    verify_solution_completeness(TTS, S, W₁₀, p₀::Vector{<:Number}, l₀)

Use the already computeded additional witnesses `W₁₀`. You want to obtain
`TTS`, `W₁₀` and `l₀` as the output from [`solution_completeness_witnesses`](@ref).

[^LRS18]:
    Leykin, Anton, Jose Israel Rodriguez, and Frank Sottile. "Trace test."
    Arnold Mathematical Journal 4.1 (2018): 113-125.
"""
function verify_solution_completeness(F, R::MonodromyResult; kwargs...)
    verify_solution_completeness(F, solutions(R), R.parameters; kwargs...)
end

function verify_solution_completeness(
    F,
    W₀₁::AbstractVector{<:AbstractVector},
    p₀::Vector{<:Number};
    show_progress = true,
    parameters = nothing,
    trace_tol = 1e-6,
    kwargs...,
)
    W₁₀, TTS, l₀ = solution_completeness_witnesses(
        F,
        W₀₁,
        p₀;
        show_progress = show_progress, parameters = parameters, kwargs...,
    )
    verify_solution_completeness(TTS, W₀₁, W₁₀, p₀, l₀; show_progress = show_progress)
end

function verify_solution_completeness(
    TTS,
    W₀₁::AbstractVector{<:AbstractVector},
    W₁₀::AbstractVector{<:AbstractVector},
    p₀::Vector{<:Number},
    l₀;
    show_progress = true,
    trace_tol = 1e-6,
    kwargs...,
)
    # Combine W₀₁ and W₁₀.
    S = append!([[x; 0.0] for x in W₀₁], W₁₀)
    # To verify that we found all solutions we need move in the pencil
    if show_progress
        @info("Compute trace...")
    end

    trace = track_and_compute_trace(TTS, S, l₀; kwargs...)
    if isnothing(trace)
        return nothing
    else
        show_progress && @info("Norm of trace: $(LinearAlgebra.norm(trace))")
        LinearAlgebra.norm(trace) < trace_tol
    end
end

"""
    solution_completeness_witnesses(F, S, p; parameters=..., kwargs...)

Compute the additional necessary witnesses. Returns a triple `(W₁₀, TTS, l)`
containing the additional witnesses `W₁₀`, a trace test system `TTS` and
the parameters `l` for `TTS`.
"""
function solution_completeness_witnesses(
    F,
    W₀₁,
    p₀::Vector{<:Number};
    parameters = nothing,
    show_progress = true,
    kwargs...,
)
    # generate another start pair
    q₀ = randn(ComplexF64, length(p₀))
    # Construct the trace test system
    TTS = TraceTestSystem(SPSystem(F; parameters = parameters), p₀, q₀ - p₀)

    y₁ = solution(solve(F, W₀₁[1]; p₁ = p₀, p₀ = q₀, parameters = parameters, kwargs...)[1])
    # construct an affine hyperplane l(x) going through y₀
    l₁ = cis.(2π .* rand(length(y₁)))
    push!(l₁, -sum(l₁ .* y₁))
    # This is numerically sometimes not so nice. Let's move to a truly generic one.
    l₀ = randn(ComplexF64, length(l₁))
    y₀ = solution(solve(TTS, [y₁; 1]; p₁ = l₁, p₀ = l₀)[1])

    if show_progress
        @info("Compute additional witnesses for completeness...")
    end

    R₁₀ = monodromy_solve(TTS, y₀, l₀; max_loops_no_progress = 5, kwargs...)
    best_result = R₁₀
    best_params = l₀
    result_agreed = 0
    for i = 1:10
        k₀ = randn(ComplexF64, length(l₀))
        S_k₀ = solutions(solve(
            TTS,
            solutions(R₁₀);
            start_parameters = l₀, target_parameters = k₀,
        ))
        new_result = monodromy_solve(TTS, S_k₀, k₀; max_loops_no_progress = 5)
        if nsolutions(new_result) == nsolutions(best_result)
            result_agreed += 1
        elseif nsolutions(new_result) > nsolutions(best_result)
            best_result = new_result
            best_params = k₀
        end
        if result_agreed > 2
            break
        end
    end

    W₁₀ = solutions(best_result)

    if show_progress
        @info("Found $(length(W₁₀)) additional witnesses")
    end
    W₁₀, TTS, best_params
end

function track_and_compute_trace(TTS::TraceTestSystem, S, l₀; kwargs...)
    for i = 1:3
        TTP = TraceTestPencil(TTS, l₀)
        R₁ = solve(TTP, S, start_parameters = [0.0], target_parameters = [0.1], kwargs...)
        R₂ = solve(TTP, S, start_parameters = [0.0], target_parameters = [-.1], kwargs...)
        if nsolutions(R₁) ≠ length(S) || nsolutions(R₂) ≠ length(S)
            if i == 3
                printstyled("Lost solutions $i times. Abort.\n", color = :red, bold = true)
                return nothing
            end
            printstyled(
                "Lost solutions, need to recompute trace...\n",
                color = :yellow,
                bold = true,
            )
        end
        s₁ = sum(solutions(R₁))
        s = sum(S)
        s₂ = sum(solutions(R₂))
        # Due to floating point errors, we compute the RELATIVE trace
        trace = ((s₁ .- s) .- (s .- s₂)) ./ max.(abs.(s₁ .- s), 1.0)
        return trace
    end
    nothing
end


text:
"Types for arbitrary samplers, and a simple Metropolis sampling algorithm."

"""
    AbstractSampler

All samplers should subtype this, and implement the following methods
    `set_temp!(sampler::MySampler, kT::Float64)`
    `get_temp(sampler::MySampler) :: Float64`
    `get_system(sampler::MySampler) :: SpinSystem`
    `sample!(sampler::MySampler)`

    Optionally, can override behavior of these, which by default do full-system
     recalculations on each call.
    `running_energy(sampler::MySampler)`
    `running_mag(sampler::MySampler)`
"""
abstract type AbstractSampler end

running_energy(sampler::S) where {S <: AbstractSampler} = energy(get_system(sampler))
running_mag(sampler::S) where {S <: AbstractSampler} = sum(get_system(sampler))
reset_running_energy!(sampler::S) where {S <: AbstractSampler} = nothing
reset_running_mag!(sampler::S) where {S <: AbstractSampler} = nothing

"""
    thermalize!(sampler, num_samples)

`sample!` a sampler a given number of times.
"""
@inline function thermalize!(sampler::S, num_samples::Int) where {S <: AbstractSampler}
    for _ in 1:num_samples
        sample!(sampler)
    end
end

"""
    anneal!(sampler, temp_schedule, step_schedule)

`sample!` a sampler at a series of temperatures, staying at each temperature
  for the number of steps in `step_schedule`.
"""
@inline function anneal!(sampler::S,
                         temp_schedule,
                         step_schedule) where {S <: AbstractSampler}
    for (temp, num_steps) in zip(temp_schedule, step_schedule)
        set_temp!(sampler, temp)
        thermalize!(sampler, num_steps)
    end
end

"""
    anneal!(sampler, temp_function, num_samples)

`sample!` a sampler `num_samples` times, with the sample at timestep `n`
 drawn at a temperature `temp_function(n)`.
"""
@inline function anneal!(sampler::S,
                         temp_function::Function,
                         num_samples::Int) where {S <: AbstractSampler}
    for t in 1:num_steps
        set_temp!(sampler, temp_function(t))
        sample!(sampler)
    end
end

"""
    MetropolisSampler(sys::SpinSystem, kT::Float64, nsweeps::Int)

A sampler which performs the standard Metropolis Monte Carlo algorithm to
 sample a `SpinSystem` at the requested temperature.

Each single-spin update attempts to move the spin to a random position on
 the unit sphere. One call to `sample!` will attempt to flip each spin
 `nsweeps` times.
"""
mutable struct MetropolisSampler <: AbstractSampler
    system     :: SpinSystem
    β          :: Float64
    nsweeps    :: Int
    E          :: Float64
    M          :: Vec3
    function MetropolisSampler(sys::SpinSystem, kT::Float64, nsweeps::Int)
        @assert kT != 0. "Temperature must be nonzero!"
        new(sys, 1.0 / kT, nsweeps, energy(sys), sum(sys))
    end
end

"""
    IsingSampler(sys::SpinSystem, kT::Float64, nsweeps::Int)

A sampler which performs the standard Metropolis Monte Carlo algorithm to
 sample a `SpinSystem` at the requested temperature.

This version differs from `MetropolisSampler` in that each single-spin update
only attempts to completely flip the spin. One call to `sample!` will attempt
to flip each spin `nsweeps` times.

Before construting, be sure that your `SpinSystem` is initialized so that each
spin points along its "Ising-like" axis.
"""
mutable struct IsingSampler <: AbstractSampler
    system     :: SpinSystem
    β          :: Float64
    nsweeps    :: Int
    E          :: Float64
    M          :: Vec3
    function IsingSampler(sys::SpinSystem, kT::Float64, nsweeps::Int)
        @assert kT != 0. "Temperature must be nonzero!"
        new(sys, 1.0 / kT, nsweeps, energy(sys), sum(sys))
    end
end

"""
    set_temp!(sampler, kT)

Changes the temperature of the sampler to `kT`.
"""
function set_temp!(sampler::MetropolisSampler, kT)
    sampler.β = 1 / kT
end
function set_temp!(sampler::IsingSampler, kT)
    sampler.β = 1 / kT
end

"""
    get_temp(sampler) :: Float64

Returns the temperature of the sampler, as `kT`.
"""
get_temp(sampler::MetropolisSampler) = 1 / sampler.β
get_temp(sampler::IsingSampler) = 1 / sampler.β

"""
    get_system(sampler)

Returns the `SpinSystem` being updated by the `sampler`.
"""
get_system(sampler::MetropolisSampler) = sampler.system
get_system(sampler::IsingSampler) = sampler.system

@inline function _random_spin() :: Vec3
    n = randn(Vec3)
    return n / norm(n)
end

"""
    sample!(sampler)

Samples `sampler.system` to a new state, under the Boltzmann distribution
 as defined by `sampler.system.hamiltonian`.
"""
function sample!(sampler::MetropolisSampler)
    for _ in 1:sampler.nsweeps
        for idx in CartesianIndices(sampler.system)
            # Try to rotate this spin to somewhere randomly on the unit sphere
            newspin = _random_spin()
            ΔE = local_energy_change(sampler.system, idx, newspin)

            if rand() < exp(-sampler.β * ΔE)
                sampler.E += ΔE
                sampler.M += (newspin - sampler.system[idx])
                sampler.system[idx] = newspin
            end
        end
    end
end

function sample!(sampler::IsingSampler)
    for _ in 1:sampler.nsweeps
        for idx in CartesianIndices(sampler.system)
            # Try to completely flip this spin
            newspin = -sampler.system[idx]
            ΔE = local_energy_change(sampler.system, idx, newspin)

            if rand() < exp(-sampler.β * ΔE)
                sampler.E += ΔE
                sampler.M += 2 * newspin
                sampler.system[idx] = newspin
            end
        end
    end
end

@inline running_energy(sampler::MetropolisSampler) = sampler.E
@inline running_mag(sampler::MetropolisSampler) = sampler.M
@inline reset_running_energy!(sampler::MetropolisSampler) = (sampler.E = energy(sampler.system); nothing)
@inline reset_running_mag!(sampler::MetropolisSampler) = (sampler.M = sum(sampler.system); nothing)
@inline running_energy(sampler::IsingSampler) = sampler.E
@inline running_mag(sampler::IsingSampler) = sampler.M
@inline reset_running_energy!(sampler::IsingSampler) = (sampler.E = energy(sampler.system); nothing)
@inline reset_running_mag!(sampler::IsingSampler) = (sampler.M = sum(sampler.system); nothing)

"""
    local_energy_change(sys, idx, newspin)

Computes the change in energy if we replace the spin at `sys[idx]`
  with `newspin`.
"""
function local_energy_change(sys::SpinSystem, idx, newspin::Vec3)
    ℋ = sys.hamiltonian
    ΔE = 0.0
    oldspin = sys[idx]
    spindiff = newspin - oldspin
    (i, cell) = splitidx(idx)

    if !isnothing(ℋ.ext_field)
        ΔE -= ℋ.ext_field.effBs[i] ⋅ spindiff
    end
    for heisen in ℋ.heisenbergs
        J = first(heisen.bondtable.data)
        for (bond, _) in sublat_bonds(heisen.bondtable, i)
            if bond.i == bond.j && iszero(bond.n)
                ΔE += J * (newspin⋅newspin - oldspin⋅oldspin)
            else
                Sⱼ = sys[bond.j, offset(cell, bond.n, sys.lattice.size)]
                ΔE += J * (spindiff ⋅ Sⱼ)
            end
        end
    end
    for diag_coup in ℋ.diag_coups
        for (bond, J) in sublat_bonds(diag_coup.bondtable, i)
            if bond.i == bond.j && iszero(bond.n)
                ΔE += newspin⋅(J.*newspin) - oldspin⋅(J.*oldspin)
            else
                Sⱼ = sys[bond.j, offset(cell, bond.n, sys.lattice.size)]
                ΔE += (J .* spindiff) ⋅ Sⱼ
            end
        end
    end
    for gen_coup in ℋ.gen_coups
        for (bond, J) in sublat_bonds(gen_coup.bondtable, i)
            if bond.i == bond.j && iszero(bond.n)
                ΔE += dot(newspin, J, newspin) - dot(oldspin, J, oldspin)
            else
                Sⱼ = sys[bond.j, offset(cell, bond.n, sys.lattice.size)]
                ΔE += dot(spindiff, J, Sⱼ)
            end
        end
    end
    if !isnothing(ℋ.dipole_int)
        throw("Local energy changes not implemented yet for dipole interactions")
    end
    return ΔE
end


text:
## general data and market data #################################

α = 0.995  ## confidence level

lob_names = [:motor_liab, :motor_other, :mat,
             :fire, :liability, :credit,
             :legal, :assist, :misc,
             :re_np_cas, :re_np_mat, :re_np_prop]
lob_ids = 1:12

## risk free interest curve:
rf = [0.03, 0.031, 0.0315, 0.032, 0.032, 0.032]

## data for premium risk and reserve risk
df_lobs_prem_risk = DataFrame(
  id =                lob_ids,
  abbrev =            lob_names,
  prem_gross_vc =     [0.10,         0.08,           0.15,
                       0.08,         0.14,           0.12,
                       0.07,         0.09,           0.13,
                       0.17,         0.17,           0.17],
  re_np_prem =        [0.8,          1.0,            1.0,
                       0.8,          0.8,            1.0,
                       1.0,          1.0,            1.0,
                       1.0,          1.0,            1.0],
  res_net_vc =        [0.09,         0.8,            0.11,
                       0.10,         0.11,           0.19,
                       0.12,         0.20,           0.20,
                       0.20,         0.20,           0.20]
  )

## correlation matrix for premium risk and reserve risk
corr_lob =
  [1.00 0.50 0.50 0.25 0.50 0.25 0.50 0.25 0.50 0.25 0.25 0.25;
   0.50 1.00 0.25 0.25 0.25 0.25 0.50 0.50 0.50 0.25 0.25 0.25;
   0.50 0.25 1.00 0.25 0.25 0.25 0.25 0.50 0.50 0.25 0.50 0.25;
   0.25 0.25 0.25 1.00 0.25 0.25 0.25 0.50 0.50 0.25 0.50 0.50;
   0.50 0.25 0.25 0.25 1.00 0.50 0.50 0.25 0.50 0.50 0.25 0.25;
   0.25 0.25 0.25 0.25 0.50 1.00 0.50 0.25 0.50 0.50 0.25 0.25;
   0.50 0.50 0.25 0.25 0.50 0.50 1.00 0.25 0.50 0.50 0.25 0.25;
   0.25 0.50 0.50 0.50 0.25 0.25 0.25 1.00 0.50 0.25 0.25 0.50;
   0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 1.00 0.25 0.50 0.25;
   0.25 0.25 0.25 0.25 0.50 0.50 0.50 0.25 0.25 1.00 0.25 0.25;
   0.25 0.25 0.50 0.50 0.25 0.25 0.25 0.25 0.50 0.25 1.00 0.25;
   0.25 0.25 0.25 0.50 0.25 0.25 0.25 0.50 0.25 0.25 0.25 1.00]

## correlation between premium and reserve risk
corr_prem_res = [1.0 0.5; 0.5 1.0]

## data for catastrophe liability risk
df_cat_liability =
  DataFrame(index = collect(1:6),
            description =
              [:l_malpractice,
               :l_employer,
               :l_director_officer,
               :l_personal,
               :l_other,
               :l_nprob_re],
            f_liab = [1.0, 1.6, 1.6, 0.0, 1.0, 2.1])
## correlations for types of catastrophe liability risk
corr_liab =
  [1.00 0.00 0.50 0.00 0.25 0.50;
   0.00 1.00 0.00 0.00 0.25 0.50;
   0.50 0.00 1.00 0.00 0.25 0.50;
   0.00 0.00 0.00 1.00 0.00 0.00;
   0.25 0.25 0.25 0.00 1.00 0.50;
   0.50 0.50 0.50 0.00 0.50 1.00]

## correlations for premium-reserve-risk, lapse-risk, cat-risk
corr_scr = [1.00 0.00 0.25;
            0.00 1.00 0.00;
            0.25 0.00 1.00]

## company data #################################################

## theft -> misc (we classify theft insurance as miscellaneous)
select_lob_names = [:fire, :liability, :misc]

df_lobs = DataFrame(
  name = select_lob_names,              # lobs of X-AG
  index = findin(lob_names, select_lob_names), # ids of lobs
  prem_gross_w_py = [500., 250.,  50.], # gross written prem. PY
  prem_gross_w_cy = [600., 300., 100.], # gross written prem. CY
  upr_py = [50., 20., 5.],              # unearned prem. res. PY
  upr_cy = [70., 40., 12.],             # unearned prem. res. CY
  upr_ny = [ 75., 45., 15.],            # unearned prem. res. NY
  re_prop_q_py = [0.25, 0.20, 0.20],    # prop. reins. ceded PY
  re_prop_q_cy = [0.25, 0.20, 0.20],    # prop. reins. ceded CY
  re_prop_q_ny = [0.25, 0.20, 0.20],    # prop. reins. ceded NY
  re_np = [true, false, false],         # has non-prop. reins.
  res_undisc_py = [200., 150., 20.],    # undisc. res. PY (EoY)
  )

β_raw = DataFrame(
  fire =      [1.0, 0.6, 0.1, 0.0, 0.0, 0.0],
  liability = [1.0, 0.9, 0.6, 0.2, 0.1, 0.0],
  misc =      [1.0, 0.4, 0.0, 0.0, 0.0, 0.0]
  )

## for fire catastrophe risk: total insured sum within the
## the circle of rations 200m that has the largest total
## insured sum.
cat_fire_is = [20, 30, 10]

## 50% employer's liability, 10% D&O, 40% personal liability
cat_liability_mix = [0.0, 0.5, 0.1, 0.4, 0.0, 0.0]


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
import Pkg.Types: VersionSpec

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, we can increment the minor version
# e.g. go from 200.600.300 to 200.601.300.
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 2.7.0 to 200.690.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.

name = "Calcium"
upstream_version = v"0.4.1"
build_for_julia16_or_newer = true
version_offset = build_for_julia16_or_newer ? v"0.0.1" : v"0.0.0"
version = VersionNumber(upstream_version.major * 100 + version_offset.major,
                        upstream_version.minor * 100 + version_offset.minor,
                        upstream_version.patch * 100 + version_offset.patch)

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/fredrik-johansson/calcium/archive/refs/tags/$(upstream_version).tar.gz",
                  "5fbc997e8c9e76c88cd85c12a86f0f14c4ebe602e9f7f11e11f0ca1f89c5d81c"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd calcium*/

if [[ ${target} == *musl* ]]; then
   export CFLAGS=-D_GNU_SOURCE=1
elif [[ ${target} == *mingw* ]]; then
   # /lib is hardcoded in many places
   sed -i -e "s#/lib\>#/$(basename ${libdir})#g" configure
   # MSYS_NT-6.3 is not detected as MINGW
   extraflags=--build=MINGW${nbits}
fi

./configure --prefix=$prefix --disable-static --enable-shared --with-gmp=$prefix --with-mpfr=$prefix --with-flint=$prefix --with-arb=$prefix --with-antic=$prefix ${extraflags}
make -j${nproc}
make install LIBDIR=$(basename ${libdir})
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=build_for_julia16_or_newer)

# The products that we will ensure are always built
products = [
    LibraryProduct("libcalcium", :libcalcium)
]

# Dependencies that must be installed before this package can be built
f(v::VersionNumber) = build_for_julia16_or_newer ?  VersionNumber(v.major, v.minor, v.patch + 1) : v
dependencies = [
    Dependency("FLINT_jll"; compat = "~$(f(v"200.800.100"))"),
    Dependency("Arb_jll", compat = "~$(f(v"200.2000"))"),
    Dependency("Antic_jll", compat = "~$(f(v"0.200.500"))"),
    Dependency("GMP_jll", build_for_julia16_or_newer ? v"6.2.0" : v"6.1.2"),
    Dependency("MPFR_jll", build_for_julia16_or_newer ? v"4.1.1" : v"4.0.2"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat = build_for_julia16_or_newer ? "1.6" : "1.0-1.5")

max_stars_repo_name: cometscome/Yggdrasil

max_issues_repo_name: cometscome/Yggdrasil

max_forks_repo_name: cometscome/Yggdrasil


text:
xprecision = 3
yprecision = 5
x = round.([1, 2, 3, 1e11],xprecision)
y = round.([1, 1.4142135623730951, 1.7320508075688772, 316227.76601683791],yprecision)
writedlm("filename", [x y], '\t')

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test
using AlgebraicNumbers
import AlgebraicNumbers.inv_totient

function test1(n)
	coeff = rand(1:10,n+1)
	a = AlgebraicNumber(coeff, Complex{BigFloat}(0.0), BigFloat(0.0))
	a.apprx = prec_roots(a.coeff)[rand(1:n)]
	calc_precision!(a)
	a = simplify(a)
	@show a.coeff, convert(Complex{Float64},a.apprx), a.prec
	b = root(a,2)
	#b = a*a
	@show b.coeff, convert(Complex{Float64},b.apprx), b.prec
	#c = root(b,2)
	c = b*b
	#c = pow2(b)
	@show c.coeff, convert(Complex{Float64},c.apprx), c.prec
end

function test2(n)
	coeff = rand(1:10,n+1)
	a = AlgebraicNumber(coeff, BigFloat(0.0), BigFloat(0.0))
	a.apprx = roots(a.coeff)[rand(1:n)]
	calc_precision!(a)

	coeff = rand(1:10,n+1)
	b = AlgebraicNumber(coeff, BigFloat(0.0), BigFloat(0.0))
	b.apprx = roots(b.coeff)[rand(1:n)]
	calc_precision!(b)

	c = a*b
	abs(roots(c.coeff) .- c.apprx)
end

function test3()
	a = sqrt(AlgebraicNumber(2))
	b = sqrt(AlgebraicNumber(3))
	axb = a*b
	# @show axb, axb.coeff
	@test abs(AlgebraicNumbers.confirm_algnumber(axb)) < 1e-10
	apb = a+b
	# @show apb, apb.coeff
	@test abs(AlgebraicNumbers.confirm_algnumber(apb)) < 1e-10
end

function test4()
	# simple test
	sqrt2 = sqrt(AlgebraicNumber(2))
	@test sqrt2^2 == 2

	# Golden ratio
	ϕ = 1//2 + sqrt(AlgebraicNumber(5)/4)

	# As we all know, this has the property that:
	@test 1+1/ϕ == ϕ
end

function plastic_constant_test()
	# see http://mathworld.wolfram.com/PlasticConstant.html
	a = sqrt(AlgebraicNumber(69))
	n = cbrt(9-a) + cbrt(9+a)
	p = n*inv(cbrt(AlgebraicNumber(18)))

	@test p-1==1/(p^4)
	@test p+1==p^3
end

function test5()
	# just an answer to a stackoverflow Q.
	# http://math.stackexchange.com/questions/422233/how-to-find-a-minimal-polynomial-field-theory
	n = sqrt(AlgebraicNumber(9*5))-sqrt(AlgebraicNumber(4*7))+sqrt(AlgebraicNumber(35))
	d = 1-sqrt(AlgebraicNumber(5))+sqrt(AlgebraicNumber(7))
	α=n/d
	@test α.coeff == BigInt[3596, 2312, -280, -156, 19]
end

function test_abs()
	ii = sqrt(AlgebraicNumber(-1))
	@test conj(ii) == -ii
	@test abs(ii) == 1
	@test abs(AlgebraicNumber(-7//8))==7//8
end

function test_real_imag()
	a = root(AlgebraicNumber(-1),5)
	alg_im = sqrt(AlgebraicNumber(-1))
	@test real(a) + alg_im*imag(a) == a

	x = rand(0:20)//10
	@test cos_alg(x - 1//2) == sin_alg(x)
	x = 2//1
	@test cos_alg(x - 1//2) == sin_alg(x)
end

# TODO: add some more tests
function test_pow2()
	a = AlgebraicNumber(3//2)
	@test pow2(a) == AlgebraicNumber(9//4)
end

function test_show()
	a = IOBuffer()
	show(a, sqrt(AlgebraicNumber(-1))+1)
	@test String(take!(a)) == "1.0 + 1.0im"

	a = IOBuffer()
	show(a, sqrt(AlgebraicNumber(2)))
	@test String(take!(a))[1] == '≈'
end

function test_log_alg()
	@test log_alg(exp_alg(1//9)) == 1//9

	@test log_alg(AlgebraicNumber(2)) == Nothing
end

function test_trig_alg()
	@test acos_alg(cos_alg(3//7)) == 3//7
	@test asin_alg(sin_alg(3//7)) == 3//7

	@test acos_alg(AlgebraicNumber(1)) == 0//1
	@test asin_alg(AlgebraicNumber(1)) == 1//2

	@test asin_alg(AlgebraicNumber(3//2)) == Nothing
	@test acos_alg(AlgebraicNumber(3//2)) == Nothing
end

function totient(x::T) where T <: Integer
	prod([(fac.first^(fac.second-1))*(fac.first-1) for fac in Nemo.factor(x)])
end

# Test the correctness of inv_totient for all totients up to m
function check_inv_totient(m::T) where T <: Integer
    # Lower bound on phi(n)==m,
    # And thus worst-case maximum range we need to consider.
    n = 2*m^2

    tots = [totient(i) for i=1:n]

    for i = 1:m
        _gold = findall(==(i), tots)
        _test = sort(collect(inv_totient(i)))

        length(_gold) == length(_test) || return false

        all(_gold .== _test) || return false
    end

	return true
end

function test_inv_totient(m::T) where T <: Integer
	@test check_inv_totient(m)
end

test3()
test4()
test5()
plastic_constant_test()
test_abs()
test_real_imag()
test_pow2()
test_show()
test_log_alg()
test_trig_alg()

# testcase of issue #5
@test AlgebraicNumber(1)+sqrt(AlgebraicNumber(-1)) != AlgebraicNumber(2)

# sqrt2 = root(AlgebraicNumber(2),2)
# an.p = (x^2-2)*(x^2-3)
# calc_precision!(an)
# an = simplify(an)
# @show an.p  (should be x^2-2)

# test multiplication of square roots
sqrt2 = root(AlgebraicNumber(2),2)
sqrt3 = root(AlgebraicNumber(3),2)
sqrt6=sqrt2*sqrt3
sqrt6_ = root(AlgebraicNumber(6),2)
@test sqrt6 == sqrt6_

#
an=root(root(AlgebraicNumber(3),2) + AlgebraicNumber(-1),2)
b = an*an
@test abs(AlgebraicNumbers.confirm_algnumber(b)) < 1e-10


text:
export
    InstantAlgorithm,
    infer

"""
    InstantAlgorithm

    Algorithm that runs once on an `InstantNetwork`.
"""
abstract type InstantAlgorithm <: Algorithm end

# FIXME:
# placeholder_beliefs is intended to provide information about the distribution of placeholders to an algorithm.
# The current interface allows marginal distributions over placeholders to be passed, which is fine for an
# algorithm like BP, but if an algorithm needs joint information it has to get it in a different way.
# In particular, during particle filtering, the importance sampling algorithm needs joint samples of the previous state.
# Fortunately, it can get it out of the :particles state in the runtime, but this is clunky and not in line with the
# intentions of the general infer method.
# We should introduce sfuncs that define distributions over dictionaries, and provide operators to produce joint samples
# or compute marginals over individual variables, as they are able to support them. placeholder_beliefs should be
# such an sfunc.
"""
    infer(algorithm::InstantAlgorithm, runtime::InstantRuntime,
        evidence::Dict{Symbol, Score},
        interventions::Dict{Symbol, Dist},
        placeholder_beliefs::Dict{Symbol, Dist})

Run the inference algorithm.

Stores the results in `runtime`. The format of these results is up to
`algorithm`, but they should be usable by queries with this algorithm.

# Arguments
- `algorithm`: The instant algorithm to run.
- `runtime`: The runtime in which to run the algorithm.
- `evidence`: The supplied evidence, which defaults to `Dict()`.
- `interventions`: The supplied interventions, which defaults to `Dict()`.
- `placeholder_beliefs`: The beliefs associated with the placeholders in the
network, which default to `Dict()`. Instant algorithms might require that a belief be
supplied for every placeholder in `network`.
"""


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AirSpaceImpl

export
    initialize,
    update,

    updateObjectState,
    getAllObjectStates,
    getAll,
    updateAll,

    AirSpace,
    ASWMState

using AbstractWorldModelImpl
using AbstractWorldModelInterfaces
using CommonInterfaces

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractWorldModelInterfaces.updateObjectState
import AbstractWorldModelInterfaces.getAllObjectStates
import AbstractWorldModelInterfaces.getAll
import AbstractWorldModelInterfaces.updateAll

type ASWMState

    t::Float64
    x::Float64
    y::Float64
    h::Float64
    vx::Float64
    vy::Float64
    vh::Float64
end

type AirSpace <: AbstractWorldModel

    number_of_aircraft::Int

    states::Vector{ASWMState}
    states_update::Vector{ASWMState}

    function AirSpace(number_of_aircraft::Int)

        obj = new()

        obj.number_of_aircraft = number_of_aircraft

        obj.states = Array(ASWMState, number_of_aircraft)
        obj.states_update = Array(ASWMState, number_of_aircraft)

        return obj
    end
end

function initialize(as::AirSpace, aircraft_number::Int, state)
  initialize(as,aircraft_number,convert(ASWMState, state))
end

function initialize(as::AirSpace, aircraft_number::Int, state::ASWMState)

    as.states[aircraft_number] = state
end

function updateObjectState(as::AirSpace, aircraft_number::Int, state::ASWMState)

    as.states_update[aircraft_number] = state
end

update(as::AirSpace, aircraft_number::Int, state) = update(as, aircraft_number, convert(ASWMState,state))

update(as::AirSpace, aircraft_number::Int, state::ASWMState) = updateObjectState(as, aircraft_number, state)

function getAllObjectStates(as::AirSpace)

    return as.states
end

getAll(as::AirSpace) = getAllObjectStates(as)

function updateAll(as::AirSpace)

    as.states = copy(as.states_update)
end

end


text:
using APIKeyHandlers
using Test

const servicename = "APIAPI"
const serviceurl = "https://api.api.com"
const servicekey = "mykeyiscool123"

g = GenericKeyHandler(servicename,serviceurl,servicekey)

@testset "Getting" begin
    @test get_name(g) == servicename
    @test get_api_url(g) == serviceurl
    @test get_key(g) == servicekey
end

@testset "Setting" begin
    set_name!(g,"API1")
    set_api_url!(g,"https://api2.api.com")
    set_key!(g,"coolerkey2")
    @test get_name(g) == "API1"
    @test get_api_url(g) == "https://api2.api.com"
    @test get_key(g) == "coolerkey2"
end

@testset "IO" begin
    g = GenericKeyHandler(servicename,serviceurl,servicekey)
    save!(g,abspath(joinpath(@__DIR__, "..", "test", "config.ini")))

    g2 = GenericKeyHandler("name","url","key")

    load!(g2,abspath(joinpath(@__DIR__, "..", "test", "config.ini")))
    @test get_name(g2) == servicename
    @test get_api_url(g2) == serviceurl
    @test get_key(g2) == servicekey

    rm(abspath(joinpath(@__DIR__, "..", "test", "config.ini")))
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libsvm"
version = v"3.24.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/cjlin1/libsvm/archive/v324.tar.gz", "3ba1ac74ee08c4dd57d3a9e4a861ffb57dab88c6a33fd53eac472fc84fbb2a8f"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libsvm-324/
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
mkdir -p ${bindir}
mkdir -p ${libdir}
if [[ "${target}" == *-freebsd* ]] || [[ "${target}" == *-apple-* ]]; then
    CC=gcc
    CXX=g++
fi
make
make lib
cp libsvm.${dlext} ${libdir}
cp svm-train${exeext} ${bindir}
cp svm-predict${exeext} ${bindir}
cp svm-scale${exeext} ${bindir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libsvm", :libsvm),
    ExecutableProduct("svm-scale", :svm_scale),
    ExecutableProduct("svm-train", :svm_train),
    ExecutableProduct("svm-predict", :svm_predict)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Flutes
using Test

@test round(Flutes.soundspeed(0.0); sigdigits=6) == 331.316
@test round(Flutes.soundspeed(25.0); sigdigits=6) == 346.146
@test Flutes.tone("A4") == 440.0
@test Flutes.tone("A5"; A=442.0) == 884.0
@test Flutes.tone("C♯0") == 17.323914
@test Flutes.tone("C♭8") == 3951.06641

max_stars_repo_name: isaacwhanson/scad-flute

max_issues_repo_name: isaacwhanson/scad-flute

max_forks_repo_name: isaacwhanson/scad-flute


text:
function set_matrixterms()

end


text:
"""
    nll(L::LPDO{MPS},data::Matrix{Pair{String,Pair{String, Int}}})

Compute the negative log-likelihood of process data using a MPO ansatz
for the vectorized unitary operator.
"""
function nll(L::LPDO{MPS}, data::Matrix{Pair{String,Pair{String,Int}}})
  data_in = first.(data)
  data_out = convertdatapoints(last.(data))

  ψ = L.X
  N = length(ψ)
  loss = 0.0
  s_in = [firstind(ψ[j]; tags="Input") for j in 1:length(ψ)]
  s_out = [firstind(ψ[j]; tags="Output") for j in 1:length(ψ)]

  for n in 1:size(data_in)[1]
    x_in = data_in[n, :]
    x_out = data_out[n, :]

    ψx = dag(ψ[1]) * dag(state(x_in[1], s_in[1]))
    ψx = ψx * state(x_out[1], s_out[1])
    for j in 2:N
      ψ_r = dag(ψ[j]) * dag(state(x_in[j], s_in[j]))
      ψ_r = ψ_r * state(x_out[j], s_out[j])
      ψx = ψx * ψ_r
    end
    prob = abs2(ψx[])
    loss -= log(prob) / size(data_in)[1]
  end
  return loss
end

nll(ψ::MPS, data::Matrix{Pair{String,Pair{String,Int}}}) = nll(LPDO(ψ), data)

"""
    nll(L::LPDO{MPO},data::Matrix{Pair{String,Pair{String, Int}}})

Compute the negative log-likelihood of process data using a LPDO ansatz
for the Choi matrix.
"""
function nll(L::LPDO{MPO}, data::Matrix{Pair{String,Pair{String,Int}}})
  # if the MPO in LPDO{MPO} is a unitary MPO (instead of a MPO with Kraus index)
  # then transform the MPO to MPS and run the nll on that
  #if !hastags(L.X[1], default_purifier_tags)
  #  return nll(_UnitaryMPOtoMPS(copy(L.X)), data)
  #end

  data_in = first.(data)
  data_out = convertdatapoints(last.(data))

  ρ = L.X
  N = length(ρ)
  loss = 0.0
  s_in = [firstind(ρ[j]; tags="Input") for j in 1:N]
  s_out = [firstind(ρ[j]; tags="Output") for j in 1:N]

  for n in 1:size(data_in)[1]
    x_in = data_in[n, :]
    x_out = data_out[n, :]
    ρdag = dag(copy(ρ))
    for j in 1:N
      ρdag[j] = ρdag[j] * dag(state(x_in[j], s_in[j]))
      ρdag[j] = ρdag[j] * state(x_out[j], s_out[j])
    end
    prob = inner(ρdag, ρdag)
    loss -= log(real(prob)) / size(data_in)[1]
  end
  return loss
end

"""
    TP(L::LPDO)

Γ = 1/√D * √(Tr[Φ²] - 2*Tr[Φ] + D)
"""
function TP(L::LPDO)
  Λ = copy(L)
  normalize!(Λ; localnorm=2)
  Φ = tr(Λ, "Output")

  D = 2^length(Φ)
  @assert D ≈ tr(Φ)
  Γ = (1 / sqrt(D)) * sqrt(inner(Φ, Φ) - D)
  return real(Γ)
end

"""
    gradnll(L::LPDO{MPS},
            data::Matrix{Pair{String,Pair{String, Int}}};
            sqrt_localnorms = nothing)

Compute the gradients of the negative log-likelihood for process data
using a MPS representation of the vectorized unitary operation.
"""
function gradnll(
  L::LPDO{MPS}, data::Matrix{Pair{String,Pair{String,Int}}}; sqrt_localnorms=nothing
)
  data_in = first.(data)
  data_out = convertdatapoints(last.(data))

  ψ = L.X
  N = length(ψ)

  s_in = [firstind(ψ[j]; tags="Input") for j in 1:length(ψ)]
  s_out = [firstind(ψ[j]; tags="Output") for j in 1:length(ψ)]

  links = [linkind(ψ, n) for n in 1:(N - 1)]

  ElT = eltype(ψ[1])

  nthreads = Threads.nthreads()

  L = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Lψ = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  R = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Rψ = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  P = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  for nthread in 1:nthreads
    for n in 1:(N - 1)
      L[nthread][n] = ITensor(ElT, undef, links[n])
      Lψ[nthread][n] = ITensor(ElT, undef, s_in[n], s_out[n], links[n])
    end
    Lψ[nthread][N] = ITensor(ElT, undef, s_in[N], s_out[N])

    for n in N:-1:2
      R[nthread][n] = ITensor(ElT, undef, links[n - 1])
      Rψ[nthread][n] = ITensor(ElT, undef, links[n - 1], s_in[n], s_out[n])
    end
    Rψ[nthread][1] = ITensor(ElT, undef, s_in[1], s_out[1])

    for n in 1:N
      P[nthread][n] = ITensor(ElT, undef, s_in[n], s_out[n])
    end
  end

  if isnothing(sqrt_localnorms)
    sqrt_localnorms = ones(N)
  end

  ψdag = dag(ψ)

  gradients = [[ITensor(ElT, inds(ψ[j])) for j in 1:N] for _ in 1:nthreads]

  grads = [[ITensor(ElT, undef, inds(ψ[j])) for j in 1:N] for _ in 1:nthreads]

  loss = zeros(nthreads)

  Threads.@threads for n in 1:size(data_in)[1]
    nthread = Threads.threadid()

    x_in = data_in[n, :]
    x_out = data_out[n, :]

    """ LEFT ENVIRONMENTS """
    P[nthread][1] = dag(state(x_in[1], s_in[1])) * state(x_out[1], s_out[1])
    L[nthread][1] .= ψdag[1] .* P[nthread][1]
    for j in 2:(N - 1)
      P[nthread][j] = dag(state(x_in[j], s_in[j])) * state(x_out[j], s_out[j])
      Lψ[nthread][j] .= L[nthread][j - 1] .* ψdag[j]
      L[nthread][j] .= Lψ[nthread][j] .* P[nthread][j]
    end
    P[nthread][N] = dag(state(x_in[N], s_in[N])) * state(x_out[N], s_out[N])
    Lψ[nthread][N] .= L[nthread][N - 1] .* ψdag[N]
    ψx = (Lψ[nthread][N] * P[nthread][N])[]
    prob = abs2(ψx)
    loss[nthread] -= log(prob) / size(data_in)[1]

    #""" RIGHT ENVIRONMENTS """
    R[nthread][N] .= ψdag[N] .* P[nthread][N]
    for j in reverse(2:(N - 1))
      Rψ[nthread][j] .= ψdag[j] .* R[nthread][j + 1]
      R[nthread][j] .= Rψ[nthread][j] .* P[nthread][j]
    end

    """ GRADIENTS """
    # TODO: fuse into one call to mul!
    grads[nthread][1] .= P[nthread][1] .* R[nthread][2]
    gradients[nthread][1] .+= (1 / (sqrt_localnorms[1] * ψx)) .* grads[nthread][1]
    for j in 2:(N - 1)
      Rψ[nthread][j] .= L[nthread][j - 1] .* P[nthread][j]
      # TODO: fuse into one call to mul!
      grads[nthread][j] .= Rψ[nthread][j] .* R[nthread][j + 1]
      gradients[nthread][j] .+= (1 / (sqrt_localnorms[j] * ψx)) .* grads[nthread][j]
    end
    grads[nthread][N] .= L[nthread][N - 1] .* P[nthread][N]
    gradients[nthread][N] .+= (1 / (sqrt_localnorms[N] * ψx)) .* grads[nthread][N]
  end

  for nthread in 1:nthreads
    for g in gradients[nthread]
      g .= (-2 / size(data_in)[1]) .* g
    end
  end

  gradients_tot = [ITensor(ElT, inds(ψ[j])) for j in 1:N]
  loss_tot = 0.0
  for nthread in 1:nthreads
    gradients_tot .+= gradients[nthread]
    loss_tot += loss[nthread]
  end

  return gradients_tot, loss_tot
end

"""
    gradnll(L::LPDO{MPO},
            data::Matrix{Pair{String,Pair{String, Int}}};
            sqrt_localnorms = nothing)

Compute the negative log-likelihood of process data using a LPDO ansatz
for the Choi matrix.
"""
function gradnll(
  L::LPDO{MPO}, data::Matrix{Pair{String,Pair{String,Int}}}; sqrt_localnorms=nothing
)
  data_in = first.(data)
  data_out = convertdatapoints(last.(data))

  ρ = L.X
  N = length(ρ)

  s_in = [firstind(ρ[j]; tags="Input") for j in 1:length(ρ)]
  s_out = [firstind(ρ[j]; tags="Output") for j in 1:length(ρ)]

  links = [linkind(ρ, n) for n in 1:(N - 1)]

  ElT = eltype(ρ[1])

  kraus = Index[]
  for j in 1:N
    push!(kraus, firstind(ρ[j], default_purifier_tags))
  end

  nthreads = Threads.nthreads()

  L = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Lρ = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Lgrad = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  R = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Rρ = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  Agrad = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  T = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Tp = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  grads = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  gradients = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  P = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  for nthread in 1:nthreads
    for n in 1:(N - 1)
      L[nthread][n] = ITensor(ElT, undef, links[n]', links[n])
    end
    for n in 2:(N - 1)
      Lρ[nthread][n] = ITensor(ElT, undef, kraus[n], links[n]', links[n - 1])
    end
    for n in 1:(N - 2)
      Lgrad[nthread][n] = ITensor(ElT, undef, links[n], kraus[n + 1], links[n + 1]')
    end
    Lgrad[nthread][N - 1] = ITensor(ElT, undef, links[N - 1], kraus[N])

    for n in N:-1:2
      R[nthread][n] = ITensor(ElT, undef, links[n - 1]', links[n - 1])
    end
    for n in (N - 1):-1:2
      Rρ[nthread][n] = ITensor(ElT, undef, links[n - 1]', kraus[n], links[n])
    end

    Agrad[nthread][1] = ITensor(ElT, undef, kraus[1], links[1]', s_in[1], s_out[1])
    for n in 2:(N - 1)
      Agrad[nthread][n] = ITensor(
        ElT, undef, links[n - 1], kraus[n], links[n]', s_in[n], s_out[n]
      )
    end

    T[nthread][1] = ITensor(ElT, undef, kraus[1], links[1])
    Tp[nthread][1] = prime(T[nthread][1], "Link")
    for n in 2:(N - 1)
      T[nthread][n] = ITensor(ElT, undef, kraus[n], links[n], links[n - 1])
      Tp[nthread][n] = prime(T[nthread][n], "Link")
    end
    T[nthread][N] = ITensor(ElT, undef, kraus[N], links[N - 1])
    Tp[nthread][N] = prime(T[nthread][N], "Link")

    grads[nthread][1] = ITensor(ElT, undef, links[1], kraus[1], s_in[1], s_out[1])
    gradients[nthread][1] = ITensor(ElT, links[1], kraus[1], s_in[1], s_out[1])
    for n in 2:(N - 1)
      grads[nthread][n] = ITensor(
        ElT, undef, links[n], links[n - 1], kraus[n], s_in[n], s_out[n]
      )
      gradients[nthread][n] = ITensor(
        ElT, links[n], links[n - 1], kraus[n], s_in[n], s_out[n]
      )
    end
    grads[nthread][N] = ITensor(ElT, undef, links[N - 1], kraus[N], s_in[N], s_out[N])
    gradients[nthread][N] = ITensor(ElT, links[N - 1], kraus[N], s_in[N], s_out[N])

    for n in 1:N
      P[nthread][n] = ITensor(ElT, undef, s_in[n], s_out[n])
    end
  end

  if isnothing(sqrt_localnorms)
    sqrt_localnorms = ones(N)
  end

  loss = zeros(nthreads)

  Threads.@threads for n in 1:size(data_in)[1]
    nthread = Threads.threadid()

    x_in = data_in[n, :]
    x_out = data_out[n, :]

    """ LEFT ENVIRONMENTS """
    P[nthread][1] = dag(state(x_in[1], s_in[1])) * state(x_out[1], s_out[1])
    P[nthread][1] = dag(P[nthread][1])
    T[nthread][1] .= ρ[1] .* P[nthread][1]
    L[nthread][1] .= prime(T[nthread][1], "Link") .* dag(T[nthread][1])
    for j in 2:(N - 1)
      P[nthread][j] = dag(state(x_in[j], s_in[j])) * state(x_out[j], s_out[j])
      P[nthread][j] = dag(P[nthread][j])
      T[nthread][j] .= ρ[j] .* P[nthread][j]
      Lρ[nthread][j] .= prime(T[nthread][j], "Link") .* L[nthread][j - 1]
      L[nthread][j] .= Lρ[nthread][j] .* dag(T[nthread][j])
    end
    P[nthread][N] = dag(state(x_in[N], s_in[N])) * state(x_out[N], s_out[N])
    P[nthread][N] = dag(P[nthread][N])
    T[nthread][N] .= ρ[N] .* P[nthread][N]
    prob = L[nthread][N - 1] * prime(T[nthread][N], "Link")
    prob = prob * dag(T[nthread][N])
    prob = real(prob[])
    loss[nthread] -= log(prob) / size(data_in)[1]

    """ RIGHT ENVIRONMENTS """
    R[nthread][N] .= prime(T[nthread][N], "Link") .* dag(T[nthread][N])
    for j in reverse(2:(N - 1))
      Rρ[nthread][j] .= prime(T[nthread][j], "Link") .* R[nthread][j + 1]
      R[nthread][j] .= Rρ[nthread][j] .* dag(T[nthread][j])
    end

    """ GRADIENTS """

    Tp[nthread][1] .= prime(ρ[1], "Link") .* P[nthread][1]
    Agrad[nthread][1] .= Tp[nthread][1] .* dag(P[nthread][1])
    grads[nthread][1] .= R[nthread][2] .* Agrad[nthread][1]
    gradients[nthread][1] .+= (1 / (sqrt_localnorms[1] * prob)) .* grads[nthread][1]
    for j in 2:(N - 1)
      Tp[nthread][j] .= prime(ρ[j], "Link") .* P[nthread][j]
      Lgrad[nthread][j - 1] .= L[nthread][j - 1] .* Tp[nthread][j]
      Agrad[nthread][j] .= Lgrad[nthread][j - 1] .* dag(P[nthread][j])
      grads[nthread][j] .= R[nthread][j + 1] .* Agrad[nthread][j]
      gradients[nthread][j] .+= (1 / (sqrt_localnorms[j] * prob)) .* grads[nthread][j]
    end
    Tp[nthread][N] .= prime(ρ[N], "Link") .* P[nthread][N]
    Lgrad[nthread][N - 1] .= L[nthread][N - 1] .* Tp[nthread][N]
    grads[nthread][N] .= Lgrad[nthread][N - 1] .* dag(P[nthread][N])
    gradients[nthread][N] .+= (1 / (sqrt_localnorms[N] * prob)) .* grads[nthread][N]
  end

  for nthread in 1:nthreads
    for g in gradients[nthread]
      g .= (-2 / size(data_in)[1]) .* g
    end
  end

  gradients_tot = Vector{ITensor}(undef, N)
  gradients_tot[1] = ITensor(ElT, links[1], kraus[1], s_in[1], s_out[1])
  for n in 2:(N - 1)
    gradients_tot[n] = ITensor(ElT, links[n], links[n - 1], kraus[n], s_in[n], s_out[n])
  end
  gradients_tot[N] = ITensor(ElT, links[N - 1], kraus[N], s_in[N], s_out[N])

  loss_tot = 0.0
  for nthread in 1:nthreads
    gradients_tot .+= gradients[nthread]
    loss_tot += loss[nthread]
  end

  return gradients_tot, loss_tot
end

"""
    gradTP(L::LPDO, gradlogZ::Vector{<:ITensor},
           logZ::Float64; sqrt_localnorms = nothing)

Compute the gradients of the trace-preserving regularization.
"""
function gradTP(
  L::LPDO, gradlogZ::Vector{<:ITensor}, logZ::Float64; sqrt_localnorms=nothing
)
  N = length(L)
  D = 2^N

  gradients_TrΦ², trΦ² = grad_TrΦ²(L; sqrt_localnorms=sqrt_localnorms)

  trΦ = exp(logZ)
  @assert D ≈ trΦ
  Γ = (1 / sqrt(D)) * sqrt(trΦ² - D)

  gradients = Vector{ITensor}(undef, N)

  for j in 1:N
    grad = gradients_TrΦ²[j] - 2 * trΦ² * gradlogZ[j]
    gradients[j] = (1 / D) * grad / (2.0 * Γ)
  end
  return gradients, Γ
end

function grad_TrΦ²(L::LPDO{MPS}; sqrt_localnorms=nothing)
  N = length(L)
  Ψ = copy(L.X)
  Ψdag = dag(Ψ)

  if isnothing(sqrt_localnorms)
    sqrt_localnorms = ones(N)
  end

  L = Vector{ITensor}(undef, N - 1)
  R = Vector{ITensor}(undef, N)

  L[1] = Ψdag[1] * prime(prime(Ψ[1], "Link"), "Input")
  L[1] = L[1] * prime(prime(Ψdag[1]), "Link")
  L[1] = L[1] * prime(prime(Ψ[1], "Output"), 3, "Link")

  for j in 2:(N - 1)
    L[j] = L[j - 1] * Ψdag[j]
    L[j] = L[j] * prime(prime(Ψ[j], "Link"), "Input")
    L[j] = L[j] * prime(prime(Ψdag[j]), "Link")
    L[j] = L[j] * prime(prime(Ψ[j], "Output"), 3, "Link")
  end
  trΦ² = L[N - 1] * Ψdag[N]
  trΦ² = trΦ² * prime(prime(Ψ[N], "Link"), "Input")
  trΦ² = trΦ² * prime(prime(Ψdag[N]), "Link")
  trΦ² = trΦ² * prime(prime(Ψ[N], "Output"), 3, "Link")
  trΦ² = real(trΦ²[])

  R[N] = Ψdag[N] * prime(prime(Ψ[N], "Link"), "Input")
  R[N] = R[N] * prime(prime(Ψdag[N]), "Link")
  R[N] = R[N] * prime(prime(Ψ[N], "Output"), 3, "Link")

  for j in reverse(2:(N - 1))
    R[j] = R[j + 1] * Ψdag[j]
    R[j] = R[j] * prime(prime(Ψ[j], "Link"), "Input")
    R[j] = R[j] * prime(prime(Ψdag[j]), "Link")
    R[j] = R[j] * prime(prime(Ψ[j], "Output"), 3, "Link")
  end

  gradients = Vector{ITensor}(undef, N)
  tmp = prime(Ψ[1], 3, "Link") * R[2]
  tmp = tmp * prime(prime(Ψdag[1], 2, "Link"), "Input")
  gradients[1] = (prime(prime(Ψ[1], "Link"), "Input") * tmp) / (sqrt_localnorms[1])

  for j in 2:(N - 1)
    tmp = prime(Ψ[j], 3, "Link") * L[j - 1]
    tmp = tmp * prime(prime(Ψdag[j], 2, "Link"), "Input")
    tmp = prime(prime(Ψ[j], "Link"), "Input") * tmp
    gradients[j] = (tmp * R[j + 1]) / (sqrt_localnorms[j])
  end
  tmp = prime(Ψ[N], 3, "Link") * L[N - 1]
  tmp = prime(prime(Ψdag[N], 2, "Link"), "Input") * tmp
  gradients[N] = (prime(prime(Ψ[N], "Link"), "Input") * tmp) / (sqrt_localnorms[N])

  return 4 * gradients, trΦ²
end

function grad_TrΦ²(Λ::LPDO{MPO}; sqrt_localnorms=nothing)
  N = length(Λ)

  if isnothing(sqrt_localnorms)
    sqrt_localnorms = ones(N)
  end

  L = Vector{ITensor}(undef, N - 1)
  R = Vector{ITensor}(undef, N)
  L[1] = bra(Λ, 1) * noprime(ket(Λ, 1); tags="Output")
  L[1] = L[1] * prime(bra(Λ, 1)', "Link")
  L[1] = L[1] * prime(prime(noprime(ket(Λ, 1), "Input"), 2, "Link"), default_purifier_tags)
  for j in 2:(N - 1)
    L[j] = L[j - 1] * bra(Λ, j)
    L[j] = L[j] * noprime(ket(Λ, j); tags="Output")
    L[j] = L[j] * prime(bra(Λ, j)', "Link")
    L[j] =
      L[j] * prime(prime(noprime(ket(Λ, j), "Input"), 2, "Link"), default_purifier_tags)
  end
  trΦ² = L[N - 1] * bra(Λ, N)
  trΦ² = trΦ² * noprime(ket(Λ, N); tags="Output")
  trΦ² = trΦ² * prime(bra(Λ, N)', "Link")
  trΦ² = trΦ² * prime(prime(noprime(ket(Λ, N), "Input"), 2, "Link"), default_purifier_tags)
  trΦ² = real(trΦ²[])

  R[N] = bra(Λ, N) * noprime(ket(Λ, N); tags="Output")
  R[N] = R[N] * prime(bra(Λ, N)', "Link")
  R[N] = R[N] * prime(prime(noprime(ket(Λ, N), "Input"), 2, "Link"), default_purifier_tags)

  for j in reverse(2:(N - 1))
    R[j] = R[j + 1] * bra(Λ, j)
    R[j] = R[j] * noprime(ket(Λ, j); tags="Output")
    R[j] = R[j] * prime(bra(Λ, j)', "Link")
    R[j] =
      R[j] * prime(prime(noprime(ket(Λ, j), "Input"), 2, "Link"), default_purifier_tags)
  end

  gradients = Vector{ITensor}(undef, N)
  tmp = prime(noprime(ket(Λ, 1)), 3, "Link") * R[2]
  tmp = tmp * prime(prime(bra(Λ, 1), "Input"), 2, "Link")
  gradients[1] = (noprime(ket(Λ, 1), "Output") * tmp) / (sqrt_localnorms[1])

  for j in 2:(N - 1)
    tmp = prime(noprime(ket(Λ, j)), 3, "Link") * L[j - 1]
    tmp = tmp * prime(prime(bra(Λ, j), "Input"), 2, "Link")
    tmp = noprime(ket(Λ, j), "Output") * tmp
    gradients[j] = (tmp * R[j + 1]) / (sqrt_localnorms[j])
  end
  tmp = prime(noprime(ket(Λ, N)), 3, "Link") * L[N - 1]
  tmp = prime(prime(bra(Λ, N), "Input"), 2, "Link") * tmp
  gradients[N] = (noprime(ket(Λ, N), "Output") * tmp) / (sqrt_localnorms[N])

  return 4 * gradients, trΦ²
end

"""
    PastaQ.gradients(L::LPDO, data::Array; sqrt_localnorms = nothing)
    PastaQ.gradients(ψ::MPS, data::Array; localnorms = nothing)

Compute the gradients of the cost function:
`C = log(Z) - ⟨log P(σ)⟩_data + TP`
"""
function gradients(
  L::LPDO,
  data::Matrix{Pair{String,Pair{String,Int}}};
  sqrt_localnorms=nothing,
  trace_preserving_regularizer=nothing,
  kwargs...
)
  g_logZ, logZ = gradlogZ(L; sqrt_localnorms=sqrt_localnorms)
  g_nll, NLL = gradnll(L, data; sqrt_localnorms=sqrt_localnorms)
  g_TP, TP_distance = gradTP(L, g_logZ, logZ; sqrt_localnorms=sqrt_localnorms)

  grads = g_logZ + g_nll
  loss = logZ + NLL

  # trace-preserving regularization
  if !isnothing(trace_preserving_regularizer)
    grads += trace_preserving_regularizer * g_TP
  end

  # TODO: check if this slows down the training
  # permute the gradients
  for j in 1:length(L)
    grads[j] = permute(grads[j], inds(L.X[j])...)
  end

  return grads, loss
end


text:
using GLMakie, DifferentialEquations, ParameterizedFunctions
GLMakie.activate!()
let
  g = @ode_def begin
    dx = σ*(y - x)
    dy = x*(ρ - z) - y
    dz = x*y - β*z
  end σ ρ β
  u0 = [1.0; 0.0; 0.0]
  tspan = (0.0, 80.0)
  p = [10.0,28.0,8/3]
  prob = ODEProblem(g, u0, tspan, p)
  sol = solve(prob, Tsit5(), saveat = 0.005)
  x, y, z = sol[1,:], sol[2,:], sol[3,:]
  tempo = sol.t
  fig = Figure(resolution=(1000,600), fontsize = 20)
  ax = Axis3(fig, aspect = (1,1,0.5), azimuth = -0.3π, elevation = π/9)
  pltobj = lines!(ax, x,y,z, color = tempo, overdraw = false)
  cbar = Colorbar(fig, pltobj, label = "time", width = 15, ticksize=15,
                  tickalign = 1, height = Relative(0.5))
  fig[1,1] = ax
  fig[1,2] = cbar
  fig
  save(joinpath(@__DIR__, "output", "lorenzAttractor.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie", "DifferentialEquations", "ParameterizedFunctions"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
function parse_spec()
    paths = Dict{String, Any}()
    if !isdir(TMP)
        return paths
    end
    for filename in readdir(TMP)
        ext = split(filename, ".")[end]
        if (ext == "yml")
            for (path, val) in YAML.load_file(joinpath(TMP, filename), dicttype=Dict{String,Any})
                if !haskey(paths, path)
                    paths[path] = Dict{String, Any}()
                end
                for (method, v) in val
                    paths[path][method] = v
                end
            end
        end
    end
    # clean up
    rm(TMP, force=true, recursive=true)
    return paths
end