text:
export monodromy_solve, realsolutions, nreal, parameters

#####################
# 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,
    :minimal_number_of_solutions, :maximal_number_of_iterations_without_progress]

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
    minimal_number_of_solutions::Int
    maximal_number_of_iterations_without_progress::Int
end

function MonodromyOptions(isrealsystem;
    distance=euclidean_distance,
    identical_tol::Float64=1e-6,
    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=isrealsystem,
    check_startsolutions=true,
    # stopping heuristic
    target_solutions_count=nothing,
    timeout=float(typemax(Int)),
    minimal_number_of_solutions::Int=default_minimal_number_of_solutions(target_solutions_count),
    maximal_number_of_iterations_without_progress::Int=10)

    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),
        minimal_number_of_solutions,
        maximal_number_of_iterations_without_progress)
end

default_minimal_number_of_solutions(::Nothing) = 1
function default_minimal_number_of_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 isrealvector(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 retcode == CoreTrackerStatus.success
        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 in length(stats.nsolutions_development):-1:1
        if stats.nsolutions_development[i] != nsolutions
            return k
        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

#####################
# Nodes of the Loop
#####################

"""
    Node(p, x::UniquePoints; store_points=true, is_main_node=false)

Create a node with a parameter from the same type as `p` and
points from `x`. If `stores_points` is `true` a `UniquePoints`
data structure is allocated to keep track of all known solutions of this node.
"""
struct Node{T, UP<:UniquePoints}
    p::Vector{T}
    points::Union{Nothing, UP}
    # Metadata / configuration
    main_node::Bool
end

function Node(p::AbstractVector{T}, x::UP; store_points=true, is_main_node=false) where {T, UP<:UniquePoints}
    if store_points == false
        Node{T, typeof(x)}(Vector(p), nothing, is_main_node)
    else
        Node(Vector(p), x, is_main_node)
    end
end
function Node(p::AbstractVector{T}, node::Node{T,UP}, options::MonodromyOptions; store_points=true, is_main_node=false) where {T, UP}
    if store_points == false
        Node{T, UP}(Vector(p), nothing, is_main_node)
    else
        Node{T, UP}(Vector(p), similar(UP), is_main_node)
    end
end

"""
    add!(node::Node, x; kwargs...)

Calls [`add!`](@ref) on the points of the Node with option `Val(true)`.
"""
function add!(node::Node, x; kwargs...)
    if node.points === nothing
        nothing
    else
        add!(node.points, x, Val(true); kwargs...)
    end
end

#####################
# Edges of the Loop
#####################
"""
    Edge(start::Int, target::Int; usegamma=false)

Create an edge between two nodes references by the index. `usegamma` refers
to an optional weight of two random complex numbers γ=(γ₁, γ₀). These are the same
gammas as in [`ParameterHomotopy`](@ref).

    Edge(edge::Edge)

Construct an indentical `Edge` as `edge`, but with newly samples γ (if applicable).
"""
struct Edge
    start::Int
    target::Int
    γ::Union{Nothing, NTuple{2, ComplexF64}}
end
function Edge(start::Int, target::Int; usegamma=false)
    if usegamma
        γ = (randn(ComplexF64), randn(ComplexF64))
    else
        γ = nothing
    end
    Edge(start, target, γ)
end
function Edge(edge::Edge)
    if edge.γ === nothing
        edge
    else
        γ = (randn(ComplexF64), randn(ComplexF64))
        Edge(edge.start, edge.target, γ)
    end
end

#######################
# Loop data structure
#######################
"""
    Loop(p, x::UniquePoints, nnodes::Int, options::MonodromyOptions; usegamma=true)

Construct a loop using `nnodes` nodes with parameters of the type of `p` and solutions `x`. `usegamma` refers to the use weights on the edges. See also
[`Edge`](@ref).

    Loop(style::LoopStyle, p, x::UniquePoints, options::MonodromyOptions)

Construct a loop using the defined style `style` with parameters of the type of `p` and solutions `x`.
"""
struct Loop{N<:Node}
    nodes::Vector{N}
    edges::Vector{Edge}
end

function Loop(p₁, x₁::UP, nnodes::Int, options::MonodromyOptions; usegamma=true) where {UP<:UniquePoints}
    n₁ = Node(p₁, x₁, is_main_node = true)
    nodes = [n₁]
    store_points = options.group_action_on_all_nodes && has_group_actions(options)
    for i = 2:nnodes
        p = options.parameter_sampler(p₁)
        push!(nodes, Node(p, n₁, options; store_points=store_points, is_main_node=false))
    end

    loop = map(i -> Edge(i - 1, i; usegamma=usegamma), 2:nnodes)
    push!(loop, Edge(nnodes, 1; usegamma=usegamma))

    Loop(nodes, loop)
end

"""
    solutions(loop::Loop)

Get the solutions of the loop.
"""
solutions(loop::Loop) = loop.nodes[1].points

"""
    nsolutions(loop::Loop)

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

"""
    mainnode(loop::Loop)

Get the main [`Node`](@ref) of the loop, i.e., the one with the original provided parameter.
"""
mainnode(loop::Loop) = loop.nodes[1]

nextedge(loop::Loop, edge::Edge) = loop.edges[edge.target]

##############
# 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 Loop(strategy::Triangle, p, x::UP, options::MonodromyOptions) where {UP<:UniquePoints}
    Loop(p, x, 3, options, usegamma=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 Loop(strategy::Petal, p, x::UP, options::MonodromyOptions) where {UP<:UniquePoints}
    Loop(p, x, 2, options, usegamma=true)
end

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

Regenerate all random parameters in the loop in order to introduce a new monodromy action.
"""
function regenerate!(loop::Loop, 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::AbstractVector, edge::Edge, loop::Loop, stats::MonodromyStatistics)

Track `x` along the edge `edge` in the loop `loop` using `tracker`. Record statistics
in `stats`.
"""
function track(tracker::CoreTracker, x::AbstractVector, edge::Edge, loop::Loop, stats::MonodromyStatistics)
    set_parameters!(tracker, edge, loop)
    track(tracker, x, stats)
end
function track(tracker::CoreTracker, x::AbstractVector, stats::MonodromyStatistics)
    retcode = track!(tracker, x, 1.0, 0.0)
    trackedpath!(stats, retcode)
    retcode
end

"""x
    set_parameters!(tracker::CoreTracker, e::Edge, loop::Loop)

Setup the parameters in the ParameterHomotopy in `tracker` to fit the edge `e`.
"""
function set_parameters!(tracker::CoreTracker, e::Edge, loop::Loop)
    H = basehomotopy(tracker.homotopy)
    if !(H isa ParameterHomotopy)
        error("Base homotopy is not a ParameterHomotopy")
    end
    p = (loop.nodes[e.start].p, loop.nodes[e.target].p)
    set_parameters!(H, p, e.γ)
end

#############
## Results ##
#############
struct MonodromyResult{N, T1, T2}
    returncode::Symbol
    solutions::Vector{SVector{N, T1}}
    parameters::Vector{T2}
    statistics::MonodromyStatistics
    equivalence_classes::Bool
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")
end

TreeViews.hastreeview(::MonodromyResult) = true
TreeViews.numberofnodes(::MonodromyResult) = 5
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, "Tracked paths")
    elseif i == 5
        print(io, "Parameters")
    end
end
function TreeViews.treenode(r::MonodromyResult, i::Integer)
    if i == 1
        return r.solutions
    elseif i == 2
        return realsolutions(r)
    elseif i == 3
        return r.returncode
    elseif i == 4
        return r.statistics.ntrackedpaths
    elseif i == 5
        return r.parameters
    end
    missing
end

"""
    mapresults(f, result::MonodromyResult; onlyreal=false, realtol=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).
realsolutions = mapresults(solution, R, onlyreal=true)
```
"""
function mapresults(f, R::MonodromyResult;
    onlyreal=false, realtol=1e-6)
    [f(r) for r in R.solutions if
        (!onlyreal || isrealvector(r, realtol))]
end

"""
    solutions(result::MonodromyResult; onlyreal=false, realtol=1e-6)

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

## Example
```julia
realsolutions = solutions(R, onlyreal=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)

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

Returns the solutions of `res` whose imaginary part has norm less than 1e-6.
"""
function realsolutions(res::MonodromyResult; tol=1e-6)
    map(r -> real.(r), filter(r -> LinearAlgebra.norm(imag.(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 -> LinearAlgebra.norm(imag.(r)) < tol, res.solutions)
end

"""
    parameters(r::MonodromyResult)

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

#####################
## monodromy solve ##
#####################
"""
MonodromyCache{FT<:FixedHomotopy, Tracker<:CoreTracker, NC<:AbstractNewtonCache}

Cache for monodromy loops.
"""
struct MonodromyCache{FT<:FixedHomotopy, Tracker<:CoreTracker, NC<:AbstractNewtonCache, AV<:AbstractVector}
    F::FT
    tracker::Tracker
    newton_cache::NC
    out::AV
end

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

## Options
* `target_solutions_count=nothing`: The computations are stopped if this number of solutions is reached.
* `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.
* `maximal_number_of_iterations_without_progress::Int=10`: The maximal number of iterations (i.e. loops generated) without any progress.
* `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.
* `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.
* `show_progress=true`: Enable a progress meter.
* `distance_function=euclidean_distance`: The distance function used for [`UniquePoints`](@ref).
* `identical_tol::Float64=1e-6`: The tolerance with which it is decided whether two solutions are identical.
* `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.
* `group_action_on_all_nodes=false`: By default the group_action(s) are only applied on the solutions with the main parameter `p`. If this is enabled then it is applied for every parameter `q`.
* `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.
* `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.
* `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.
* `timeout=float(typemax(Int))`: The maximal number of *seconds* the computation is allowed to run.
* `minimal_number_of_solutions`: The minimal number of solutions before a stopping heuristic is applied. By default this is half of `target_solutions_count` if applicable otherwise 2.
"""
function monodromy_solve(F::Inputs, solution::Vector{<:Number}, p₀::AbstractVector{TP}; kwargs...) where {TP}
    monodromy_solve(F, [solution], p₀; kwargs...)
end
function monodromy_solve(F::Inputs, solutions::Vector{<:AbstractVector{<:Number}}, p₀::AbstractVector{TP}; kwargs...) where {TP}
    monodromy_solve(F, static_solutions(solutions), p₀; kwargs...)
end
function monodromy_solve(F::Inputs,
        startsolutions::Vector{<:SVector{NVars, <:Complex}},
        p::AbstractVector{TP};
        parameters=nothing,
        strategy=nothing,
        show_progress=true,
        showprogress=nothing, #deprecated
        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

    p₀ = Vector{promote_type(Float64, TP)}(p)

    optionskwargs, restkwargs = splitkwargs(kwargs, monodromy_options_supported_keywords)

    # construct tracker
    tracker = coretracker(F, startsolutions; parameters=parameters, generic_parameters=p₀, restkwargs...)
    if affine_tracking(tracker)
        HC = HomotopyWithCache(tracker.homotopy, tracker.state.x, 1.0)
        F₀ = FixedHomotopy(HC, 0.0)
    else
        # Force affine newton method
        patch_state = state(EmbeddingPatch(), tracker.state.x)
        HC = HomotopyWithCache(PatchedHomotopy(tracker.homotopy, patch_state), tracker.state.x, 1.0)
        F₀ = FixedHomotopy(HC, 0.0)
    end

    # Check whether homotopy is real
    isrealsystem = numerically_check_real(tracker.homotopy, tracker.state.x)

    options = MonodromyOptions(isrealsystem; optionskwargs...)
    # construct cache
    newt_cache = newton_cache(F₀, tracker.state.x)
    C =  MonodromyCache(F₀, tracker, newt_cache, copy(tracker.state.x))
    # construct UniquePoints
    if options.equivalence_classes
        uniquepoints = UniquePoints(eltype(startsolutions), options.distance_function;
                                    group_actions = options.group_actions,
                                    check_real = true)
    else
        uniquepoints = UniquePoints(eltype(startsolutions), options.distance_function; check_real = true)
    end
    # add only unique points that are true solutions
    if options.check_startsolutions
        for s in startsolutions
            y = verified_affine_vector(C, s, s, options)
            if y !== nothing
                add!(uniquepoints, y; tol=options.identical_tol)
            end
        end
    else
        add!(uniquepoints, startsolutions; tol=options.identical_tol)
    end
    statistics = MonodromyStatistics(uniquepoints)
    if  length(points(uniquepoints)) == 0
        @warn "None of the provided solutions is a valid start solution."
        return MonodromyResult(:invalid_startvalue,
                Vector{SVector{length(startsolutions[1]),ComplexF64}}(), p₀, statistics,
                options.equivalence_classes)
    end

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

    # construct Loop
    loop = Loop(strategy, p₀, uniquepoints, options)

    # solve
    retcode = :not_assigned
    if show_progress
        if !options.equivalence_classes
            progress = make_monodromy_progress("Solutions found:")
        else
            progress = make_monodromy_progress("Classes of solutions (modulo group action) found:")
        end
    else
        progress = nothing
    end
    try
        retcode = monodromy_solve!(loop, C, options, statistics, progress)
    catch e
        if (e isa InterruptException)
            retcode = :interrupt
        else
            rethrow(e)
        end
    end
    finished!(statistics, nsolutions(loop))
    MonodromyResult(retcode, points(solutions(loop)), p₀, statistics, options.equivalence_classes)
end

function make_monodromy_progress(desc::String; dt::Real=0.1, output::IO=stderr)
    tfirst = time()
    tlast = tfirst + 0.5 # half a second delay before something is displayed
    printed = false
    ProgressMeter.ProgressUnknown(false, dt, 0, false, tfirst, tlast, printed, desc, :green, output, 0)
end

function default_strategy(F::MPPolyInputs, parameters, p₀::AbstractVector{TP}; isrealsystem=false) where {TC,TP}
    # If F depends only linearly on the parameters a petal is sufficient
    vars = variables(F; parameters=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 isrealsystem
        Triangle(useweights=true)
    else
        Triangle(useweights=false)
    end
end
function default_strategy(F::AbstractSystem, parameters, p₀; isrealsystem=false)
    # For a real system we should introduce some weights to avoid the discriminant
    if isrealsystem
        Triangle(useweights=true)
    else
        Triangle(useweights=false)
    end
end

# convert vector of vectors to vector of svectors
static_solutions(V::Vector) = static_solutions(V, Val(length(V[1])))
function static_solutions(V::Vector, ::Val{N}) where {N}
    map(v -> complex.(float.(SVector{N}(v))), V)
end
function static_solutions(V::Vector{<:AbstractVector{<:Complex{<:AbstractFloat}}}, ::Val{N}) where {N}
    SVector{N}.(V)
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 ##
#################

"""
    Job{N, T}

A `Job` is consisting of an `Edge` and a solution to the start node of this edge.
"""
struct Job{N, T}
    x::SVector{N, T}
    edge::Edge
end

function monodromy_solve!(loop::Loop, C::MonodromyCache, options::MonodromyOptions,
    stats::MonodromyStatistics, progress)

    t₀ = time_ns()
    iterations_without_progress = 0 # stopping heuristic
    # intialize job queue
    queue = map(x -> Job(x, loop.edges[1]), solutions(loop))

    n = nsolutions(loop)
    while n < options.target_solutions_count
        retcode = empty_queue!(queue, loop, C, options, t₀, stats, progress)

        if retcode == :done
            update_progress!(progress, loop, stats; finish=true)
            break
        elseif retcode == :timeout
            return :timeout
        elseif retcode == :invalid_startvalue
            return :invalid_startvalue
        end

        # Iterations heuristic
        n_new = nsolutions(loop)
        if n == n_new
            iterations_without_progress += 1
        else
            iterations_without_progress = 0
            n = n_new
        end
        if iterations_without_progress == options.maximal_number_of_iterations_without_progress &&
            n_new ≥ options.minimal_number_of_solutions
            return :heuristic_stop
        end

        regenerate_loop_and_schedule_jobs!(queue, loop, options, stats)
    end

    :success
end

function empty_queue!(queue, loop::Loop, C::MonodromyCache, options::MonodromyOptions,
        t₀::UInt, stats::MonodromyStatistics, progress)
    while !isempty(queue)
        job = pop!(queue)
        status = process!(queue, job, C, loop, options, stats, progress)
        if status == :done
            return :done
        elseif status == :invalid_startvalue
            return :invalid_startvalue
        end
        update_progress!(progress, loop, stats)
        # check timeout
        if (time_ns() - t₀) > options.timeout * 1e9 # convert s to ns
            return :timeout
        end
    end
    :incomplete
end

function verified_affine_vector(C::MonodromyCache, ŷ, x, options)
    # We distinguish solutions which have a distance larger than identical_tol
    # But due to the numerical error in the evaluation of the distance, we need to be a little bit
    # careful. Therefore, we require that the solutions should be one magnitude closer to
    # the true solutions as necessary
    tol = 0.1 * options.identical_tol
    result = newton!(C.out, C.F, ŷ, euclidean_norm, C.newton_cache,
                tol=tol, miniters=1, maxiters=3, simplified_last_step=false)

    if isconverged(result)
        return affine_chart(x, C.out)
    else
        return nothing
    end
end

affine_chart(x::SVector, y::PVector) = ProjectiveVectors.affine_chart!(x, y)
affine_chart(x::SVector{N, T}, y::AbstractVector) where {N, T} = SVector{N,T}(y)

function process!(queue::Vector{<:Job}, job::Job, C::MonodromyCache, loop::Loop, options::MonodromyOptions, stats::MonodromyStatistics, progress)
    retcode = track(C.tracker, job.x, job.edge, loop, stats)
    if retcode ≠ CoreTrackerStatus.success
        if retcode == CoreTrackerStatus.terminated_invalid_startvalue && stats.ntrackedpaths == 0
            return :invalid_startvalue
        end
        return :incomplete
    end

    node = loop.nodes[job.edge.target]

    if node.main_node
        y = verified_affine_vector(C, currx(C.tracker), job.x, options)
        # is the solution at infinity?
        if y === nothing
            return :incomplete
        end
    else
        y = affine_chart(job.x, currx(C.tracker))
    end
    next_edge = nextedge(loop, job.edge)
    add_and_schedule!(node, queue, y, options, stats, next_edge) && return :done

    if options.complex_conjugation
        add_and_schedule!(node, queue, conj.(y), options, stats, next_edge) && return :done
    end

    if !options.equivalence_classes && node.points !== nothing
        apply_actions(options.group_actions, y) do yᵢ
            add_and_schedule!(node, queue, yᵢ, options, stats, next_edge)
        end
    end
    return :incomplete
end

"""
    add_and_schedule!(node, queue::Vector{Job{N,T}}, y, options, stats, nextedge)

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!(node, queue::Vector{Job{N,T}}, y, options, stats, next_edge) where {N,T}
    k = add!(node, y; tol=options.identical_tol)
    if k == NOT_FOUND
        # Check if we are done
        node.main_node && isdone(node, y, options) && return true
        push!(queue, Job(SVector{N,T}(y), next_edge))
    elseif k == NOT_FOUND_AND_REAL
        # If we are on the main node check whether we have a real root.
        if node.main_node
            stats.nreal +=1
        end
        # Check if we are done
        node.main_node && isdone(node, y, options) && return true
        push!(queue, Job(SVector{N,T}(y), next_edge))
    elseif k === nothing
        push!(queue, Job(SVector{N,T}(y), next_edge))
    end
    false
end

function update_progress!(::Nothing, loop::Loop, statistics::MonodromyStatistics; finish=false)
    nothing
end
function update_progress!(progress, loop::Loop, statistics::MonodromyStatistics; finish=false)
    nsolutions = length(solutions(loop))
    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(node::Node, x, options::MonodromyOptions)
    !node.main_node && return false

    options.done_callback(x, node.points) ||
    length(node.points) ≥ options.target_solutions_count
end

function regenerate_loop_and_schedule_jobs!(queue, loop::Loop, options::MonodromyOptions, stats::MonodromyStatistics)
    sols = solutions(loop)
    # create a new loop by regenerating the parameters (but don't touch our
    # main node)
    regenerate!(loop, options, stats)
    for x in sols
        push!(queue, Job(x, loop.edges[1]))
    end
    nothing
end


text:
"""
Test ideas for eXtensions of the GapLas/MGT idea
"""
module TestGapLaX

using Test
import GraphIdx
import GraphIdx: num_edges, IncidenceIndex
import TreeLas

function sort_edges(nidx::IncidenceIndex)
    local m = num_edges(nidx)
    edge_rank = zeros(Int, 2m)
    for (i, (v, ei)) in enumerate(nidx.pi)
        if edge_rank[2ei-1] == 0
            edge_rank[2ei-1] = i
        else
            edge_rank[2ei-0] = i
        end
    end
    return edge_rank
end

function non_tree_idx(nidx::IncidenceIndex, selected)
    m = num_edges(nidx)
    nonidx = deepcopy(nidx)
    edge_rank = sort_edges(nidx)

    for (i, s) in enumerate(selected)
        s <= 0 && continue
        ri, rj = edge_rank[2s-1], edge_rank[2s-0]
        (u, ei) = nonidx.pi[ri]
        (v, ej) = nonidx.pi[rj]
        @assert ei == ej
        nonidx.pi[ri] = (u, 0)
        nonidx.pi[rj] = (v, 0)
    end
    error("not finished")
    return nonidx
end

function non_tree_idx(n::Int, edges::Vector{E}, selected) where {E}
    edges = copy(edges)
    for s in selected
        if s > 0
            edges[s] = (0, 0)
        end
    end
    return IncidenceIndex(n, edges)
end

function swap_to_front(f::Function, selected::Vector{E}, root::Integer) where {E}
    try
        selected[root], selected[1] = selected[1], selected[root]
        f(selected)
    finally
        selected[root], selected[1] = selected[1], selected[root]
    end
end

function non_tree_edges(edges::Vector{E}, selected::Vector{I}, root::Int) where {E,I}
    swap_to_front(selected, root) do selected::Vector{I}
        @assert selected[1] == -1
        nt = copy(edges)
        for s in @view selected[2:end]
            nt[s] = (0, 0)
        end

        local m2::Int = length(edges)
        local i::Int = 1
        while true
            while nt[i][1] > 0 && i < m2
                i += 1
            end
            while nt[m2][1] <= 0 && i <= m2
                m2 -= 1
            end
            if i < m2
                nt[i], nt[m2] = nt[m2], nt[i]
                i += 1
                m2 -= 1
            else
                return resize!(nt, m2)
            end
        end
    end
end

@testset "GapLax: demo 3x7               " begin
    grid = GraphIdx.Grid.GridGraph(3, 7)
    edges, lam = GraphIdx.Grid.collect_edges(grid)
    y = [0.62 0.73 0.71 1.5 1.17 0.43 1.08
         0.62 1.73 0.95 1.46 1.6 1.16 0.38
         0.9 0.32 -0.48 0.95 1.08 0.02 0.4]
    n, m = length(y), length(edges)
    root_node = 1
    α = zeros(Float64, m)
    γ = Vector{Float64}(undef, m)
    TreeLas.Dual.gap_vec!(γ, y, α, grid, -1.0)
    nidx = GraphIdx.IncidenceIndex(n, edges)
    pmem = GraphIdx.PrimMstMem(nidx)
    pi = GraphIdx.prim_mst_edges(γ, root_node, pmem)
    selected = pmem.selected
    cidx = GraphIdx.Tree.ChildrenIndex(pi)

    ntree_edges = non_tree_edges(edges, selected, root_node)
    @test minimum(map(e -> min(e...), ntree_edges)) >= 1
    @test maximum(map(e -> max(e...), ntree_edges)) <= m
    @test length(ntree_edges) == m - n + 1

    lcas = GraphIdx.Tree.lowest_common_ancestors(cidx, pi, ntree_edges)

    @test sortperm(nidx.pi, by=e -> e[2]) == sort_edges(nidx)

    nonidx = non_tree_idx(n, edges, selected)
    @test num_edges(nonidx) == m - n + 1
end

end

ext: jl

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
# collect(::Vector) does a copy.
# lazy_collect avoid this copy in case `v` is a Vector
lazy_collect(v::Vector) = v
lazy_collect(v) = collect(v)

# H-representation

"""
    hrep(hyperplanes::HyperPlaneIt, halfspaces::HalfSpaceIt; d::FullDim)

Creates an H-representation for the polyhedron of full dimension `d` equal to
the intersection of the hyperplanes `hyperplanes` and halfspaces `halfspaces`.

### Examples
For instance, the simplex
```math
\\begin{align*}
  x_1 + x_2 &= 1 \\\\
  x_1 &\\geq 0 \\\\
  x_2 &\\geq 0
\\end{align*}
```
can be created as follows:
```jldoctest
julia> hrep([HyperPlane([1, 1], 1)], [HalfSpace([0, -1], 0), HalfSpace([-1, 0], 0)])
H-representation Polyhedra.Intersection{Int64,Array{Int64,1},Int64}:
1-element iterator of HyperPlane{Int64,Array{Int64,1}}:
 HyperPlane([1, 1], 1),
2-element iterator of HalfSpace{Int64,Array{Int64,1}}:
 HalfSpace([0, -1], 0)
 HalfSpace([-1, 0], 0)
```
"""
function hrep(hyperplanes::HyperPlaneIt, halfspaces::HalfSpaceIt;
              d::FullDim=FullDim_rec(hyperplanes, halfspaces))
    return Intersection(d, hyperplanes, halfspaces)
end

"""
    hrep(halfspaces::HalfSpaceIt; d::FullDim)

Creates an H-representation for the polyhedron of full dimension `d` equal to
the intersection of the halfspaces `halfspaces`.

### Examples
For instance, the polytope
```math
\\begin{align*}
  x_1 + x_2 &\\leq 1 \\\\
  x_1 - x_2 &\\leq 0 \\\\
  x_1 & \\geq 0.
\\end{align*}
```
can be created as follows:
```julia
hrep([HalfSpace([1, 1], 1), HalfSpace([1, -1], 0), HalfSpace([-1, 0], 0)])
```
"""
function hrep(halfspaces::ElemIt{HalfSpace{T, AT}}; kws...) where {T, AT}
    return hrep(HyperPlane{T, AT}[], halfspaces; kws...)
end

mutable struct Intersection{T, AT, D<:FullDim} <: HRepresentation{T}
    hyperplanes::HyperPlanesIntersection{T, AT, D}
    halfspaces::Vector{HalfSpace{T, AT}}
    function Intersection{T, AT, D}(d::FullDim, hyperplanes::HyperPlaneIt{T},
                                    halfspaces::HalfSpaceIt{T}) where {T, AT, D}
        new{T, AT, D}(HyperPlanesIntersection{T, AT, D}(d, hyperplanes),
                      lazy_collect(halfspaces))
    end
end
function Intersection(d::FullDim, hyperplanes::ElemIt{HyperPlane{T, AT}},
                      halfspaces::ElemIt{HalfSpace{T, AT}}) where {T, AT}
    return Intersection{T, AT, typeof(d)}(d, hyperplanes, halfspaces)
end

FullDim(h::Intersection) = FullDim(h.hyperplanes)
hvectortype(::Type{<:Intersection{T, AT}}) where {T, AT} = AT
similar_type(PT::Type{<:Intersection}, d::FullDim, ::Type{T}) where {T} = Intersection{T, similar_type(hvectortype(PT), d, T), typeof(d)}

Intersection(h::HRepresentation{T}) where {T} = Intersection{T}(h)
Intersection{T}(h::HRepresentation) where {T} = convert(Intersection{T, similar_type(hvectortype(typeof(h)), T), typeof(FullDim(h))}, h)

function Base.intersect!(h::Intersection, hp::HyperPlane)
    intersect!(h.hyperplanes, hp)
    return h
end
function Base.intersect!(h::Intersection, hs::HalfSpace)
    push!(h.halfspaces, hs)
    return h
end

@subrepelem Intersection HyperPlane hyperplanes
@vecrepelem Intersection HalfSpace halfspaces

fulltype(::Type{<:Union{Intersection{T, AT, D}, HyperPlanesIntersection{T, AT, D}}}) where {T, AT, D} = Intersection{T, AT, D}

# V-representation

#"""
#    vrep(sympoints::SymPointIt)
#
#Creates a V-representation for the symmetric polytope equal to the convex hull of the symmetric points `sympoints`.
#
#### Examples
#The following creates a square
#```julia
#vrep([SymPoint([1, 1])], [SymPoint([1, -1])])
#```
#"""
#vrep(sympoints::SymPointIt) = SymPointsHull(sympoints)
#
#mutable struct SymPointsHull{T, AT} <: VSymPolytope{T}
#    sympoints::Vector{SymPoint{T, AT}}
#    function SymPointsHull{T, AT}(sympoints::SymPointIt{T}) where {T, AT}
#        new{T, AT}(lazy_collect(sympoints))
#    end
#end
#SymPointsHull(ps::ElemIt{SymPoint{T, AT}}) where {T, AT<:AbstractVector{T}} = SymPointsHull{T, AT}(collect(ps))
#vectortype(::Union{SymPointsHull{T, AT}, Type{SymPointsHull{T, AT}}}) where {T, AT} = AT
#similar_type(PT::Type{<:SymPointsHull}, d::FullDim, ::Type{T}) where {T} = SymPointsHull{T, similar_type(vectortype(PT), d, T)}
#
#SymPointsHull{T, AT}(sympoints::SymPointIt, points::PointIt, lines::LineIt, rays::RayIt) where {T, AT} = Hull{T, AT}(sympoints, points, lines, rays)
#SymPointsHull{T, AT}(sympoints::SymPointIt, lines::LineIt, rays::RayIt) where {T, AT} = Hull{T, AT}(sympoints, AT[], lines, rays)
#
#@vecrepelem SymPointsHull SymPoint sympoints
#
## SymPoint's can be split
#removevredundancy(vrep::SymPointsHull, hrep::HRep; kws...) = removevredundancy(PointsHull(sympoints(vrep), points(vrep)), hrep; kws...)

#"""
#    vrep(sympoints::SymPointIt, points::PointIt)
#
#Creates a V-representation for the polytope equal to the convex hull of the symmetric points `sympoints` and points `points`.
#
#### Examples
#The convex hull of ``(0, -1)``, ``(0, 1)`` and ``(1/2, 1/2)`` can be created as follows:
#```julia
#vrep([SymPoint([0, 1])], [[1/2, 1/2]])
#```
#"""
#vrep(sympoints::SymPointIt, points::PointIt) = PointsHull(sympoints, points)

"""
    vrep(points::PointIt; d::FullDim)

Creates a V-representation for the polytope of full dimension `d` equal to the
convex hull of the points `points`.

### Examples
The convex hull of ``(0, 0)``, ``(0, 1)`` and ``(1/2, 1/2)`` can be created as follows using exact arithmetic
```julia
vrep([[0, 0], [0, 1], [1//2, 1//2]])
```
or as follows using floating point arithmetic
```julia
vrep([[0, 0], [0, 1], [1/2, 1/2]])
```
"""
vrep(points::PointIt; d::FullDim = FullDim_rec(points)) = PointsHull(d, points)

mutable struct PointsHull{T, AT, D<:FullDim} <: VPolytope{T}
    d::D
    points::Vector{AT}
    function PointsHull{T, AT, D}(d::FullDim, points::PointIt) where {T, AT, D}
        new{T, AT, D}(FullDim_convert(D, d), lazy_collect(points))
    end
end
function PointsHull(d::FullDim, points::PointIt)
    return PointsHull{coefficient_type(eltype(points)), eltype(points),
                      typeof(d)}(d, points)
end
FullDim(v::PointsHull) = v.d
vvectortype(::Type{<:PointsHull{T, AT}}) where {T, AT} = AT
similar_type(PT::Type{<:PointsHull}, d::FullDim, ::Type{T}) where {T} = PointsHull{T, similar_type(vvectortype(PT), d, T), typeof(d)}

vreptype(::Type{PointsHull{T, AT, D}}) where {T, AT, D} = Hull{T, AT, D}

@vecrepelem PointsHull Point points

"""
    vrep(lines::LineIt, rays::RayIt; d::FullDim)

Creates a V-representation for the polyhedral cone of full dimension `d` equal
to the conic hull of the lines `lines` and rays `rays`.

### Examples
```julia
vrep([Line([0, 1])], [Ray([1, 0])])
```
creates a V-representation for the halfspace ``x_1 \\ge 0``.
"""
function vrep(lines::LineIt, rays::RayIt; d = FullDim_rec(lines, rays))
    return RaysHull(d, lines, rays)
end

"""
    vrep(rays::RayIt)

Creates a V-representation for the polyhedral cone of full dimension `d` equal
to the conic hull of the rays `rays`.

### Examples
```julia
vrep([Ray([1, 0]), Ray([0, 1])])
```
creates a V-representation for positive orthant.
"""
function vrep(rays::ElemIt{Ray{T, AT}}; kws...) where {T, AT}
    return vrep(Line{T, AT}[], rays; kws...)
end

mutable struct RaysHull{T, AT, D<:FullDim} <: VCone{T}
    lines::LinesHull{T, AT, D}
    rays::Vector{Ray{T, AT}}
    function RaysHull{T, AT, D}(d::FullDim, ls::LineIt{T},
                                rs::RayIt{T}) where {T, AT, D}
        new{T, AT, D}(LinesHull{T, AT, D}(d, ls), lazy_collect(rs))
    end
end
function RaysHull(d::FullDim, ls::ElemIt{Line{T, AT}}, rs::ElemIt{Ray{T, AT}}) where {T, AT}
    RaysHull{T, AT, typeof(d)}(d, ls, rs)
end
FullDim(v::RaysHull) = FullDim(v.lines)
vvectortype(::Type{<:RaysHull{T, AT}}) where {T, AT} = AT
similar_type(PT::Type{<:RaysHull}, d::FullDim, ::Type{T}) where {T} = RaysHull{T, similar_type(vvectortype(PT), d, T), typeof(d)}

@vecrepelem RaysHull Ray rays
@subrepelem RaysHull Line lines

vreptype(::Type{RaysHull{T, AT, D}}) where {T, AT, D} = Hull{T, AT, D}

"""
    vrep(points::PointIt, lines::LineIt, rays::RayIt)

Creates a V-representation for the polyhedron of full dimension `d` equal to the
minkowski sum of the convex hull of `points` with the conic hull of `lines` and
`rays`.
"""
function vrep(points::PointIt, lines::LineIt, rays::RayIt;
              d::FullDim = FullDim_rec(points, lines, rays))
    return Hull(d, points, lines, rays)
end

function vrep(points::ElemIt{AT}, lines::ElemIt{Line{T, AT}};
              kws...) where {T, AT}
    return vrep(points, lines, Ray{T, AT}[]; kws...)
end

mutable struct Hull{T, AT, D<:FullDim} <: VRepresentation{T}
    points::PointsHull{T, AT, D}
    rays::RaysHull{T, AT, D}
    function Hull{T, AT, D}(d::FullDim, vits::VIt{T}...) where {T, AT, D}
        points, lines, rays = fillvits(d, vits...)
        new{T, AT, D}(PointsHull(d, points), RaysHull(d, lines, rays))
    end
end
function Hull(d::FullDim, points::ElemIt{AT}, lines::ElemIt{Line{T, AT}},
              rays::ElemIt{Ray{T, AT}}) where {T, AT}
    Hull{T, AT, typeof(d)}(d, points, lines, rays)
end
FullDim(v::Hull) = FullDim(v.points)
vvectortype(::Type{<:Hull{T, AT}}) where {T, AT} = AT
similar_type(PT::Type{<:Hull}, d::FullDim, ::Type{T}) where {T} = Hull{T, similar_type(vvectortype(PT), d, T), typeof(d)}

Hull(v::VRepresentation{T}) where {T} = Hull{T}(v)
Hull{T}(v::VRepresentation) where {T} = convert(Hull{T, similar_type(vvectortype(typeof(v)), T), typeof(FullDim(v))}, v)

@subrepelem Hull Point points
@subrepelem Hull Line rays
@subrepelem Hull Ray rays

fulltype(::Type{<:Union{Hull{T, AT, D}, PointsHull{T, AT, D}, LinesHull{T, AT, D}, RaysHull{T, AT, D}}}) where {T, AT, D} = Hull{T, AT, D}

dualtype(::Type{<:Intersection{T, BT, D}}, ::Type{AT}) where {T, AT, BT, D} = Hull{T, AT, D}
dualtype(::Type{<:Hull{T, BT, D}}, ::Type{AT}) where {T, AT, BT, D} = Intersection{T, AT, D}
const AnyIntersection{T, AT, D} = Union{Intersection{T, AT, D}, HyperPlanesIntersection{T, AT, D}}
function dualfullspace(h::Union{AnyIntersection, Type{<:AnyIntersection}},
                       d::FullDim, ::Type{T}, ::Type{AT}) where {T, AT}
    Hull{T, AT, typeof(d)}(d, [origin(AT, fulldim(d))],
                           Line{T, AT}.(basis.(AT, Ref(d), 1:fulldim(d))),
                           Ray{T, AT}[])
end


text:
using DReal
using AbstractDomains
using Base.Test
x = Var(Int,"x")
y = Var(Int,"y")
add!((x > 2) & (y < 10) & (x + 2*y == 7))
@test is_satisfiable()

x = Var(Float64,"x",-100.0,100.0)
y = Var(Float64,"y",-100.0,100.0)
add!((x^2 + y^2 > 4) & (x^3 + y < 5))
is_satisfiable()
x,y = model(x,y)
@test (x^2 + y^2 > 4) & (x^3 + y < 5) == t

# Flyspeck
x1 = Var(Float64, "x1", 3.0, 3.14)
x2 = Var(Float64, "x2", -7.0, 5.0)
c = (2π - 2*x1*asin(cos(0.797) * sin(π/x1))) <= (-0.591 - 0.0331 *x2 + 0.506 + 1.0)
add!(c)
@test !is_satisfiable()

p = Bool('p')
q = Bool('q')
r = Bool('r')
add!(p → q & r == !q & (!p | r))
model(p,q,r)


text:
"""
# BusinessModel
## Fields
- `t::Integer`: 1 if model was just initialized, 1 + the number of time periods elapsed since model initialization otherwise. Time is considered in an abstract sense, without units.
- `p::Vector{<:AbstractFloat}`: A vector of individual probabilities (0 <= p <= 1) linked to the access rights, representing the probability that this particular access right be required by a new user.
- `pσ::Number`: The standard deviation of the normal distribution from which random samples will be taken to make `p` evolve from time period to time period. Default: 1%.
- `a::Matrix{<:Bool}`: A matrix of access rights business requirements. Rows = Users. Columns = Access Rights.
"""
mutable struct BusinessModel
    t::Int
    p::Vector{<:Number}
    pσ::Number
    a::BitArray{2}
    function BusinessModel(u:: Int, r:: Int, pσ::Number = 0.01)
        t = 1; # The first instant of the model life is 1.
        p = rand(r); # Initialize right probabilities are randomely select from uniform distribution.
        a = BitArray{2}(undef, 0, r); # Initialize the original array without any users.
        m = new(t, p, pσ, a);
        add_users(m, u);
        return m;
    end;
end;
#function Base.show(io::IO, m::BusinessModel)
    #compact = get(io, :compact, false);
    #if compact
        #print("m(t: ", m.t, ", |u|: ", count_users(m), ", |r|: ", count_rights(m), "pσ: ", m.pσ, ")");
    #else
    #    println("BusinessModel");
    #    println("t: ", m.t);
    #    println("|u|: ", count_users(m));
    #    println("|r|: ", count_rights(m));
    #    println("pσ: ", m.pσ);
    #    println("p: ", m.p);
    #end;
#end;
#function Base.show(io::IO, ::MIME"text/plain", m::BusinessModel)
#    println("MIME: text/plain");
#    show(io, m);
#end;
function count_users(m::BusinessModel)
    return size(m.a)[1];
end;
function count_rights(m::BusinessModel)
    return size(m.a)[2];
end;
function add_users(m::BusinessModel, u::Int = 1)
    # Generate a matrix of random values,
    # that will be used to randomely determine which user is granted what access right.
    up = rand(u, count_rights(m));
    # Reshape the model access right probabilities to match the size of the above random matrix.
    mp = repeat(reshape(m.p, 1, length(m.p)),u);
    # For every matrix value where the user's random value is lower than the probability,
    # consider the right as granted, otherwise not granted.
    ua = broadcast(<, up, mp);
    # Merge the new user's rights into the model.
    m.a = cat(dims=1, m.a, ua);
end;
function sample1()
    m = BusinessModel(3,4,.01);
    count_users(m);
    count_rights(m);
    m;
    m.a;
end;
#sample1();

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
baremodule Atomix

export AtomicRefArray

import Base

macro atomic end
macro atomicreplace end
macro atomicswap end

struct AtomicRef{T,Ptr,Data}
    eltype::Base.Val{T}
    ptr::Ptr
    data::Data
end

struct AtomicRefArray{T,N,Data<:Base.AbstractArray{T,N}} <:
       Base.AbstractArray{AtomicRef{T,Data},N}
    data::Data
end

function asref end

module Internal

import ..Atomix: AtomicRef, AtomicRefArray, @atomic, @atomicswap, @atomicreplace
using ..Atomix: Atomix, asref

using Base.Meta: isexpr
using Base: @propagate_inbounds
using UnsafeAtomics:
    Ordering, UnsafeAtomics, monotonic, acquire, release, acq_rel, seq_cst, right

include("utils.jl")
include("core.jl")
include("sugar.jl")
if isdefined(Base, :replaceproperty!)  # 1.7 or later
    include("properties.jl")
end

function define_docstring()
    path = joinpath(@__DIR__, "..", "README.md")
    include_dependency(path)
    doc = read(path, String)
    doc = replace(doc, r"^```julia"m => "```jldoctest README")
    # Setting `LineNumberNode` to workaround an error from logging(?) `no method
    # matching getindex(::Nothing, ::Int64)`:
    ex = :($Base.@doc $doc Atomix)
    ex.args[2]::LineNumberNode
    ex.args[2] = LineNumberNode(1, Symbol(path))
    Base.eval(Atomix, ex)
end

end  # module Internal

Internal.define_docstring()

end  # baremodule Atomix


text:
using OptiMo, Bazinga
using Printf
using CUTEst

VERBOSE = true
FREQ = 10
TOL = 1e-8
MAXIT = 1000

problems = CUTEst.select( min_var=1, max_var=100, min_con=1, max_con=100 )
problem = problems[6]

@printf "\n========================================\n"
@printf "%s\n" problem
nlp = CUTEstModel( problem )
prob = NLPOptiModel( nlp )
x0 = copy( prob.meta.x0 )
T = eltype( x0 )
nvar = prob.meta.nvar
ncon = prob.meta.ncon

# augmented lagrangian
fx0 = obj( prob, x0 )
cx0 = cons( prob, x0 )
px0 = proj( prob, cx0 )
px0 .= cx0 .- px0

ye = zeros(T,ncon)
mu = zeros(T,ncon)
Bazinga.penalty_estimate!( fx0, px0, mu, 0.1 )

alprob = AugLagOptiModel( prob, mu, ye )

# solver
solver = Bazinga.ZEROFPR( verbose=VERBOSE, freq=FREQ, tol_optim=TOL, max_iter=MAXIT )
@time out = solver( alprob )
print( out )

solver = Bazinga.PANOC( verbose=VERBOSE, freq=FREQ, tol_optim=TOL, max_iter=MAXIT )
@time out = solver( alprob )
print( out )

finalize( nlp )


text:
@doc """
# Description
Runs the CROPS algorithm using a specified cost function for a given minimum and maximum penalty value to find a range of segmentations. Function can also be invoked through @PELT macro.

# Arguments
* `segment_cost::Function`: Calculates cost between two specified indices
* `n::Int`: Length of time series
* `pen::Tuple{Real,Real}`: Tuple containing minimum and maximum penalty values

# Returns
* `out::Dict{ASCIIString,Array{T,N}}`:
  * `out["penalty"]` : vector of penalties between minimum and maximum pen
  * `out["number"]` : vector of number of changepoints for each penalty
  * `out["constrained"]` : vector of optimal segmentation costs for each penalty
  * `out["changepoints"]` : vector of changepoints sets for each penalty

# Example
```
# Sample time series with change points
n = 1000       # Length of time series
λ = 100        # Frequency of change points
μ, σ = Normal(0.0, 10.0), 1.0
sample, cps = @changepoint_sampler n λ Normal(μ, σ)

# Create cost function and run CROPS
seg_cost = NormalMeanSegment(sample)
pen = (4.0,100.0)
crops_output = CROPS(seg_cost, n, pen)
```

# See also
PELT, @PELT

# References
Haynes, K., Eckley. I.A., and Fearnhead, P., (2014) Efficient penalty search for multiple changepoint problems arXiv:1412.3617

""" ->
function CROPS(segment_cost::Function , n::Int64, pen::Tuple{Real,Real} )

    pen_interval = [ minimum(pen) , maximum(pen) ]

    # these are what we output
    out_num_cpts = Array(Int64,0)
    out_max_pen = Array(Float64,0)
    out_constrain = Array(Float64,0)
    out_cpts = Array(Array{Int64},0)

    while true

        stop = 0

        if length(pen_interval) > 0

            # do PELT and record no of chpts and
            # calculates the cost of the constrained problem
            for i in 1:length(pen_interval)

                # do PELT for each pen_interval
                cpts , opt = PELT( segment_cost , n , pen=pen_interval[i] )

                # if there are already same no of chpts in out
                # see if penalty different
                if length(cpts) ∈ out_num_cpts
                    index = findin(out_num_cpts, length(cpts))[1]

                    # check for penalty: is it bigger than max penalty
                    if pen_interval[i] > out_max_pen[index]
                        out_max_pen[index] = pen_interval[i]
                    else
                        # no change
                        stop+=1
                    end

                else
                    # number of chpts not seen before add to output and constrained likelihood

                    push!( out_num_cpts , length(cpts) )
                    push!( out_constrain , opt - (length(cpts)-1)*pen_interval[i] )
                    push!( out_max_pen , pen_interval[i])
                    push!( out_cpts, cpts)

                end

            end
	# end of if length(pen_interval) > 0
        end

        # leave loop and finish subject to output if no changes
        if stop == length(pen_interval)
            break
        end

        # look through out_num_cpts and calc beta_int and put these in pen_interval
        ord_ind = sortperm(out_num_cpts)
        out_num_cpts = out_num_cpts[ord_ind]
        out_constrain = out_constrain[ord_ind]
        out_max_pen = out_max_pen[ord_ind]
        out_cpts = out_cpts[ord_ind]

        # make pen_interval
        pen_interval = Array(Float64,0)

        for i in 1:(length(out_num_cpts)-1)

            # if they differ by 1, then just calc beta_int and extend range of pen for large no of chpts
            if out_num_cpts[i] ==  out_num_cpts[i+1] + 1
                beta_int =  (out_constrain[i+1] - out_constrain[i])
                out_max_pen[i] = beta_int
	    else
                # difference in no of chpts is >= 1 so add to pen_interval to check how many
                beta_int =  ( out_constrain[i+1] - out_constrain[i] )/( out_num_cpts[i] - out_num_cpts[i+1] )
                push!( pen_interval , beta_int )
             end

        end

       if length(pen_interval) > 0
           i = 1
           while i <= length(pen_interval)
               for j in 1:length(out_max_pen)
                   if abs(pen_interval[i] - out_max_pen[j]) < 1e-2
                       splice!(pen_interval,i)
                       i = i - 1
                       break
                   end
               end
               i = i + 1
           end

       end

    # end of while loop
    end

    # Calculate beta intervals
    nb = length(out_max_pen)
    beta_e= Array(Float64,0)
    beta_int = Array(Float64,0)
    sort_out_max_pen = sort(out_max_pen);
    sort_out_constrain = sort(out_constrain);
    sort_out_num_cpts = sort(out_num_cpts, rev = true);

    for k in 1:nb
        if k == 1
            beta_int = push!(beta_int, sort_out_max_pen[1])

        else
            beta_int = push!(beta_int, beta_e[k-1])
        end

        if k == nb
            beta_e = push!(sort_out_max_pen[k])

        else
            beta_e = push!(beta_e,(sort_out_constrain[k] - sort_out_constrain[k+1])/(sort_out_num_cpts[k+1] - sort_out_num_cpts[k]))
        end
    end

    # organise output into a dictionary
    out = Dict{ASCIIString, Array}()
    out["number"] = out_num_cpts
    out["penalty"] = sort(beta_int, rev = true)
    out["constrained"] = out_constrain
    out["changepoints"] = out_cpts
    return out

end


text:
using NumericalMethodsforEngineers
using Test

a = [10. 1. -5;-20. 3. 20.;5. 3. 5.]
b = [1., 2., 6.]

c = a\b
f = lu(a);
(lower, upper) = lufac(a)
subfor!(lower, b)
subbac!(upper, b)

@test round.(b, digits=10) == round.(c, digits=10)


text:
@testset "promote $T" for T in (Double16, Double32, Double64)

f2 = 2.0

@test promote_type(DoubleFloat{Float64}, DoubleFloat{Float32}) == DoubleFloat{Float64}
@test promote_type(DoubleFloat{Float64}, DoubleFloat{Float16}) == DoubleFloat{Float64}
@test promote_type(DoubleFloat{Float32}, DoubleFloat{Float16}) == DoubleFloat{Float32}

@test promote_type(DoubleFloat{Float64}, BigFloat) == BigFloat
@test promote_type(DoubleFloat{Float32}, BigFloat) == BigFloat
@test promote_type(DoubleFloat{Float16}, BigFloat) == BigFloat

@test promote_type(DoubleFloat{Float64}, BigInt) == BigInt
@test promote_type(DoubleFloat{Float32}, BigInt) == BigInt
@test promote_type(DoubleFloat{Float16}, BigInt) == BigInt

@test promote(DoubleFloat{Float64}(f2), DoubleFloat{Float32}(f2))[2] == DoubleFloat{Float64}(f2)
@test promote(DoubleFloat{Float64}(f2), DoubleFloat{Float16}(f2))[2] == DoubleFloat{Float64}(f2)
@test promote(DoubleFloat{Float32}(f2), DoubleFloat{Float16}(f2))[2] == DoubleFloat{Float32}(f2)

@test promote(DoubleFloat{Float64}(f2), BigFloat(f2))[2] == DoubleFloat{Float64}(f2)
@test promote(DoubleFloat{Float32}(f2), BigFloat(f2))[2] == DoubleFloat{Float32}(f2)
@test promote(DoubleFloat{Float16}(f2), BigFloat(f2))[2] == DoubleFloat{Float16}(f2)

@test promote(DoubleFloat{Float64}(f2), BigInt(f2))[2] == DoubleFloat{Float64}(f2)
@test promote(DoubleFloat{Float32}(f2), BigInt(f2))[2] == DoubleFloat{Float32}(f2)
@test promote(DoubleFloat{Float16}(f2), BigInt(f2))[2] == DoubleFloat{Float16}(f2)

end


text:
include("../src/star.jl")
using DistributedArrays

scalar_times_vector(a, x, w) = x[:] = a * w

#compute x=u'vw
utvw_star(u, v, w, x) = Star.scan(scalar_times_vector, sum(Star.reduce(dot, u, v)), x, w)

function utvw_star_test(u, v, w, x)
  #compute x=u'vw
  t = @elapsed utvw_star(u, v, w, x)

  #code to test the results
  a = convert(Array{Float64, 1}, u)
  b = convert(Array{Float64, 1}, v)
  c = convert(Array{Float64, 1}, w)
  d = convert(Array{Float64, 1}, x)

  e = zeros(Float64, length(a))
  t1 = @elapsed begin
    e = dot(a, b) * c
  end

  println("parallel time ", t, " serial time ", t1)
  println("rel norm ", norm(d - e)/norm(e))
end

function utvw_no_star(u, v, w, x)
  remote_ref = [@spawnat p dot(localpart(u), localpart(v)) for p in workers()]
  partial_dot_product = [fetch(r) for r in remote_ref]
  s = sum(partial_dot_product)
  @sync [@spawnat p scalar_times_vector(s, localpart(x), localpart(w)) for p in workers()]
end


text:
#export makeAxisArray

# Support for handling complex datatypes in HDF5 files
function writeComplexArray(file, dataset, A::AbstractArray{Complex{T},D}) where {T,D}
  d_type_compound = HDF5.h5t_create(HDF5.H5T_COMPOUND,2*sizeof(T))
  HDF5.h5t_insert(d_type_compound, "r", 0 , HDF5.hdf5_type_id(T))
  HDF5.h5t_insert(d_type_compound, "i", sizeof(T) , HDF5.hdf5_type_id(T))

  shape = collect(reverse(size(A)))
  space = HDF5.h5s_create_simple(D, shape, shape)

  dset_compound = HDF5.h5d_create(file, dataset, d_type_compound, space,
                                  HDF5.H5P_DEFAULT,HDF5.H5P_DEFAULT,HDF5.H5P_DEFAULT)
  HDF5.h5s_close(space)

  HDF5.h5d_write(dset_compound, d_type_compound, HDF5.H5S_ALL, HDF5.H5S_ALL, HDF5.H5P_DEFAULT, A)

  HDF5.h5d_close(dset_compound)
  HDF5.h5t_close(d_type_compound)
end

function isComplexArray(file, dataset)
  if eltype(file[dataset]) <: Complex
    return true
  end

  # If complex number support in HDF5 is disabled (or version < 0.12.2)
  if eltype(file[dataset]) <: HDF5.HDF5Compound{2}
    if HDF5.h5t_get_member_name(datatype(file[dataset]).id,0) == "r" &&
      HDF5.h5t_get_member_name(datatype(file[dataset]).id,1) == "i"
      return true
    end
  end
  return false
end

function getComplexType(file, dataset)
  T = HDF5.hdf5_to_julia_eltype(
            HDF5Datatype(
              HDF5.h5t_get_member_type( datatype(file[dataset]).id, 0 )
          )
        )
    return Complex{T}
end

function readComplexArray(file::HDF5File, dataset)
  T = getComplexType(file, dataset)
  A = copy(readmmap(file[dataset],Array{getComplexType(file,dataset)}))
  return A
end

function readComplexArray(filename::String, dataset)
  h5open(filename, "r") do file
    return readComplexArray(file, dataset)
  end
end

function makeAxisArray(array::Array{T,5}, pixelspacing, offset, dt) where T
  N = size(array)
  im = AxisArray(array, Axis{:color}(1:N[1]),
		 Axis{:x}(range(offset[1],step=pixelspacing[1],length=N[2])),
		 Axis{:y}(range(offset[2],step=pixelspacing[2],length=N[3])),
		 Axis{:z}(range(offset[3],step=pixelspacing[3],length=N[4])),
		 Axis{:time}(range(0*unit(dt),step=dt,length=N[5])))
  return im
end

# RawFile removed see commit https://github.com/MagneticParticleImaging/MPIFiles.jl/commit/c2b49833cc00127770fb2e1f5342883f213ff002


text:
using DataDeps

for (ver, checksum) in [("1.6", "16814254fe194d55a2fcc24858aa76d71de3c49e495bd98478cc7345e766d8b7"),
            ("1.7", "0495577ac3a87c2a64fe6189798ea046de0f44943dfb7b60fe38cf648d34c421"),
            ("1.7.1", "70b9eb7ca0dc9d67655f9d671d40be10aeff490f0bea4f10cb1946127b74c102"),
            ("2.0", "93fbae725f0125dedb7369403fda1dace85b2dcd8a523ed80af23e863b18ef2c"),
            ("3.0", "a8000014d6fc864f8bd9d83c62be601151cadd617c6554a39a1ad38b4b3f017b")]

    register(DataDep("SemCor $ver",
        """
        Website: http://web.eecs.umich.edu/%7Emihalcea/downloads.html#semcor
        Orignal Author: George A. Miller et al.
        Maintainer: Rada Mihalcea
        For WordNet version $ver

        This is SemCor orginally developed along side WordNet.
        It was automatically updated with the new numbering beyond WordNet 1.6,  by Rada Mihalcea

        Please cite the following publication if you use the corpora:
        George A. Miller, Claudia Leacock, Randee Tengi, and Ross T. Bunker. (1993). "A Semantic Concordance." In: Proceedings of the 3 DARPA Workshop on Human Language Technology.
        """,
        "http://web.eecs.umich.edu/~mihalcea/downloads/semcor/semcor$(ver).tar.gz",
        checksum;
        post_fetch_method = fn -> begin
            unpack(fn)
            innerdir = "semcor$(ver)"
            innerfiles = readdir(innerdir)
            # Move everything to current directory, under same name
            mv.(joinpath.(innerdir, innerfiles), innerfiles)
            rm(innerdir)
        end
    ))
end


text:
function keyword_check{T}(ta::TimeArray{T,1}, prices::Bool, log_transform::Bool)
    if prices == false && log_transform == false # simple returns
        return ta
    elseif prices == false && log_transform == true # log returns
        return basecall(ta,expm1)
    elseif prices == true && log_transform == false # regular prices
        return percentchange(ta)
    elseif prices == true && log_transform == true # log prices
        return percentchange(basecall(ta,exp))
    else
        error("invalid keyword argument")
    end
end


text:
using Libdl
using MPI
using MSTM
using OffsetArrays
using Test

@testset "MSTM" begin
    MPI.Init()

    mstm_lib = ("CI" => "true") in ENV ? "/home/runner/work/MSTM.jl/MSTM.jl/shared/mstm" : "../shared/mstm"
    mstm = Libdl.dlopen(mstm_lib)

    @testset "Constants" begin
        @testset "init!($notd)" for notd in [5, 10, 20, 50, 100]
            @test begin
                ctx = MSTM.Constants.init(notd)
                bcof_julia, _, _, vwh_coef_julia = MSTM.Constants.get_offset_constants(ctx)

                MSTM.Wrapper.init!(mstm, notd)
                bcof_fortran = MSTM.Wrapper.bcof(mstm, notd)
                vwh_coef_fortran = MSTM.Wrapper.vwh_coef(mstm, notd)

                isapprox(bcof_julia, bcof_fortran) && isapprox(vwh_coef_julia, vwh_coef_fortran)
            end
        end
    end

    @testset "SpecialFunctions" begin
        ϵ = 1e-8
        @testset "ricbessel($n, $ds)" for (n, ds) in [(50, 10.0), (100, 100.0), (150, 200.0)]
            @test begin
                Ψ_julia = MSTM.SpecialFunctions.ricbessel(n, ds, ϵ)
                Ψ_fortran = MSTM.Wrapper.ricbessel(mstm, n, ds, ϵ)
                isapprox(Ψ_julia, Ψ_fortran)
            end
        end

        @testset "richankel($n, $ds)" for (n, ds) in [(50, 10.0), (100, 100.0), (150, 200.0)]
            @test begin
                Ξ_julia = MSTM.SpecialFunctions.richankel(n, ds)
                Ξ_fortran = MSTM.Wrapper.richankel(mstm, n, ds)
                isapprox(Ξ_julia, Ξ_fortran)
            end
        end

        @testset "cricebessel($n, $ds)" for (n, ds) in [
            (10, ComplexF64(0, 0)),
            (20, ComplexF64(10, 10)),
            (500, ComplexF64(0.5, 5.1)),
        ]
            @test begin
                Ψ_julia = MSTM.SpecialFunctions.cricbessel(n, ds)
                Ψ_fortran = MSTM.Wrapper.cricbessel(mstm, n, ds)
                isapprox(Ψ_julia, Ψ_fortran)
            end
        end

        @testset "crichankel($n, $ds)" for (n, ds) in [
            (10, ComplexF64(0, 0)),
            (20, ComplexF64(10, 10)),
            (500, ComplexF64(0.5, 5.1)),
        ]
            @test begin
                Ξ_julia = MSTM.SpecialFunctions.crichankel(n, ds)
                Ξ_fortran = MSTM.Wrapper.crichankel(mstm, n, ds)
                isapprox(Ξ_julia, Ξ_fortran)
            end
        end

        @testset "cspherebessel($n, $z)" for (n, z) in [
            (10, ComplexF64(0, 0)),
            (20, ComplexF64(10, 10)),
            (500, ComplexF64(0.5, 5.1)),
        ]
            @test begin
                csj_julia, csy_julia = MSTM.SpecialFunctions.cspherebessel(n, z)
                csj_fortran, csy_fortran = MSTM.Wrapper.cspherebessel(mstm, n, z)
                isapprox(csj_julia, csj_fortran) && isapprox(csy_julia, csy_fortran)
            end
        end

        @testset "vcfunc($m, $n, $k, $l)" for (m, n, k, l) in [(3, 4, 3, 10), (30, 40, 30, 100), (-4, 10, 4, 5)]
            @test begin
                ctx = MSTM.Constants.init()
                vcn_julia = MSTM.SpecialFunctions.vcfunc(ctx, m, n, k, l)
                vcn_fortran = MSTM.Wrapper.vcfunc(mstm, m, n, k, l)
                isapprox(vcn_julia, vcn_fortran)
            end
        end

        @testset "normalizedlegendre($cbe, $mmax, $nmax)" for (cbe, mmax, nmax) in
                                                              [(0.5, 10, 10), (-0.21, 15, 20), (0.98, 30, 34)]
            @test begin
                ctx = MSTM.Constants.init()
                dc_julia = MSTM.SpecialFunctions.normalizedlegendre(ctx, cbe, mmax, nmax)
                dc_fortran = MSTM.Wrapper.normalizedlegendre(mstm, cbe, mmax, nmax)
                isapprox(dc_julia, dc_fortran)
            end
        end

        @testset "rotcoef($cbe, $kmax, $nmax)" for (cbe, kmax, nmax) in [(0.5, 10, 10), (-0.21, 15, 20), (0.98, 30, 34)]
            @test begin
                ctx = MSTM.Constants.init()
                dc_julia = MSTM.SpecialFunctions.rotcoef(ctx, cbe, kmax, nmax)
                dc_fortran = MSTM.Wrapper.rotcoef(mstm, cbe, kmax, nmax)
                isapprox(dc_julia, dc_fortran)
            end
        end

        @testset "taufunc($cb, $nmax)" for (cb, nmax) in [(0.5, 10), (-0.21, 20), (0.98, 34)]
            @test begin
                ctx = MSTM.Constants.init()
                τ_julia = MSTM.SpecialFunctions.taufunc(ctx, cb, nmax)
                τ_fortran = MSTM.Wrapper.taufunc(mstm, cb, nmax)
                isapprox(τ_julia, τ_fortran)
            end
        end

        @testset "pifunc($cb, $E_ϕ, $nmax, $ndim)" for (cb, E_ϕ, nmax, ndim) in [
            (0.5, 1 + 1.0im, 10, 10),
            (-0.21, 0.2 - 0.3im, 14, 20),
            (0.98, -1.2 + 2.5im, 15, 15),
        ]
            @test begin
                ctx = MSTM.Constants.init()
                π_julia = MSTM.SpecialFunctions.pifunc(ctx, cb, E_ϕ, nmax, ndim)
                π_fortran = MSTM.Wrapper.pifunc(mstm, cb, E_ϕ, nmax, ndim)
                isapprox(π_julia, π_fortran)
            end
        end

        @testset "planewavecoef($α, $β, $nodr)" for (α, β, nodr) in
                                                    [(0.5, 0.6, 10), (-0.21, -0.45, 14), (0.98, -0.45, 20)]
            @test begin
                ctx = MSTM.Constants.init()
                pmnp0_julia = MSTM.SpecialFunctions.planewavecoef(ctx, α, β, nodr)
                pmnp0_fortran = MSTM.Wrapper.planewavecoef(mstm, α, β, nodr)
                isapprox(pmnp0_julia, pmnp0_fortran)
            end
        end

        @testset "gaussianbeamcoef($α, $β, $cbeam, $nodr)" for (α, β, cbeam, nodr) in [
            (0.5, 0.6, 15.0, 10),
            (-0.21, -0.45, 0.5, 14),
            (0.98, -0.45, 100.0, 20),
        ]
            @test begin
                ctx = MSTM.Constants.init()
                pmnp0_julia = MSTM.SpecialFunctions.gaussianbeamcoef(ctx, α, β, cbeam, nodr)
                pmnp0_fortran = MSTM.Wrapper.gaussianbeamcoef(mstm, α, β, cbeam, nodr)
                isapprox(pmnp0_julia, pmnp0_fortran)
            end
        end

        @testset "sphereplanewavecoef" begin
            @test begin
                ctx = MSTM.Constants.init()

                nsphere = 3
                nodr = [3, 3, 3]
                α = 0.5
                β = 0.5
                rpos = [[0.0, 0.0, 0.0] [0.0, 0.0, -0.5] [0.0, 0.0, -1.0]]
                hostsphere = [0, 1, 2]
                numberfieldexp = [3, 3, 3]
                rimedium = [1.1 + 1.3im, 1.0 - 0.5im]

                pmnp_julia = MSTM.SpecialFunctions.sphereplanewavecoef(
                    ctx,
                    nodr,
                    α,
                    β,
                    rpos,
                    hostsphere,
                    numberfieldexp,
                    rimedium,
                )

                pmnp_fortran =
                    MSTM.Wrapper.sphereplanewavecoef(mstm, nodr, α, β, rpos, hostsphere, numberfieldexp, rimedium)

                isapprox(pmnp_julia, pmnp_fortran)
            end
        end

        @testset "axialtrancoefrecurrence($itype, $r, $nmax, $lmax)" for (itype, r, nmax, lmax) in [
            (1, 0.0, 10, 10),
            (1, 1.3, 15, 20),
            (3, 0.0, 4, 5),
            (3, 4.3, 20, 15),
        ]
            @test begin
                ri = [1.0 + 0.5im, 1.0 - 0.5im]
                ctx = MSTM.Constants.init()
                ac_julia = MSTM.SpecialFunctions.axialtrancoefrecurrence(ctx, itype, r, ri, nmax, lmax)
                ac_fortran = MSTM.Wrapper.axialtrancoefrecurrence(mstm, itype, r, ri, nmax, lmax)
                isapprox(ac_julia, ac_fortran)
            end
        end

        @testset "tranordertest($r, $ri, $lmax, $ϵ)" for (r, ri, lmax, ϵ) in [
            (0.01, 1.1 + 0.3im, 10, 1e-4),
            (0.1, 1.1 + 0.03im, 20, 1e-6),
            (0.35, 1.1 + 0.003im, 50, 1e-6),
        ]
            @test begin
                ctx = MSTM.Constants.init()
                nmax_julia = MSTM.SpecialFunctions.tranordertest(ctx, r, ri, lmax, ϵ)
                nmax_fortran = MSTM.Wrapper.tranordertest(mstm, r, ri, lmax, ϵ)
                nmax_julia == nmax_fortran
            end
        end

        @testset "gentrancoef($itype, $xptran, $nrow0, $nrow1, $ncol0, $ncol1, $iaddrow0, $iaddcol0)" for (
            itype,
            xptran,
            nrow0,
            nrow1,
            ncol0,
            ncol1,
            iaddrow0,
            iaddcol0,
        ) in [
            (1, [0.0, 0.0, 0.0], 1, 5, 1, 5, 0, 0),
            (3, [0.0, 0.0, 0.0], 1, 5, 1, 5, 0, 0),
            (1, [1.0, -0.05, 0.04], 1, 5, 1, 5, 0, 0),
            (3, [1.0, -0.05, 0.04], 1, 5, 1, 5, 0, 0),
            (1, [1.0, -0.05, 0.04], 1, 5, 1, 5, 5, 5),
            (3, [1.0, -0.05, 0.04], 1, 5, 1, 5, 2, 3),
        ]
            @test begin
                ri = [1.0 + 0.5im, 1.0 - 0.5im]
                ctx = MSTM.Constants.init()
                ac_julia = MSTM.SpecialFunctions.gentrancoef(
                    ctx,
                    itype,
                    xptran,
                    ri,
                    nrow0,
                    nrow1,
                    ncol0,
                    ncol1,
                    iaddrow0,
                    iaddcol0,
                )
                ac_fortran =
                    MSTM.Wrapper.gentrancoef(mstm, itype, xptran, ri, nrow0, nrow1, ncol0, ncol1, iaddrow0, iaddcol0)
                isapprox(ac_julia, ac_fortran)
            end
        end

        @testset "eulerrotation($xp, $eulerangf, $dir)" for (xp, eulerangf, dir) in [
            ([1.02, 2.03, 3.04], [0.95, 2.33, 1.34], 1),
            ([2.43, -23.2, -19.4], [10.2, 12.3, 11.3], -1),
            ([-0.42, -0.2, -0.001], [-23.2, -23.4, 12.2], 1),
        ]
            @test begin
                xprot_julia = MSTM.SpecialFunctions.eulerrotation(xp, eulerangf, dir)
                xprot_fortran = MSTM.Wrapper.eulerrotation(mstm, xp, eulerangf, dir)
                isapprox(xprot_julia, xprot_fortran)
            end
        end

        @testset "planewavetruncationorder($r, $rimedium, $ϵ)" for (r, rimedium, ϵ) in [
            (0.01, [1.1 + 0.3im, 0.8 - 0.2im], 1e-4),
            (0.1, [1.1 + 0.03im, 1.1 - 0.05im], 1e-6),
            (0.35, [1.1 + 0.003im, 1.2 + 1.1im], 1e-8),
        ]
            @test begin
                nodr_julia = MSTM.SpecialFunctions.planewavetruncationorder(r, rimedium, ϵ)
                nodr_fortran = MSTM.Wrapper.planewavetruncationorder(mstm, r, rimedium, ϵ)
                nodr_julia == nodr_fortran
            end
        end

        @testset "vwhcalc($rpos, $nodr, $itype)" for (rpos, nodr, itype) in [
            ([0.0, 0.0, 0.0], 5, 1),
            ([0.0, 0.0, 0.0], 5, 3),
            ([1.0, -0.05, 0.04], 10, 1),
            ([1.0, -0.05, 0.04], 10, 3),
            ([1.0, -0.05, 0.04], 20, 1),
            ([1.0, -0.05, 0.04], 20, 3),
        ]
            @test begin
                ri = [1.0 + 0.5im, 1.0 - 0.5im]
                ctx = MSTM.Constants.init()
                vwh_julia = MSTM.SpecialFunctions.vwhcalc(ctx, rpos, ri, nodr, itype)
                vwh_fortran = MSTM.Wrapper.vwhcalc(mstm, rpos, ri, nodr, itype)
                isapprox(vwh_julia, vwh_fortran)
            end
        end

        @testset "vwhaxialcalc($rpos, $nodr, $itype)" for (rpos, nodr, itype) in [
            ([0.0, 0.0, 0.0], 5, 1),
            ([0.0, 0.0, 0.0], 5, 3),
            ([1.0, -0.05, 0.04], 10, 1),
            ([1.0, -0.05, 0.04], 10, 3),
            ([1.0, -0.05, 0.04], 20, 1),
            ([1.0, -0.05, 0.04], 20, 3),
        ]
            @test begin
                ri = [1.0 + 0.5im, 1.0 - 0.5im]
                ctx = MSTM.Constants.init()
                vwh_julia = MSTM.SpecialFunctions.vwhaxialcalc(ctx, rpos, ri, nodr, itype)
                vwh_fortran = MSTM.Wrapper.vwhaxialcalc(mstm, rpos, ri, nodr, itype)
                isapprox(vwh_julia, vwh_fortran)
            end
        end

        @testset "twobytwoinverse($mat)" for mat in [
            [1.1+1.0im 2.2-1.5im; 3.3-0.2im -2.5+0.4im],
            [1.1+1.0im 1.2+1.5im; 0.3-0.2im -1.5+0.3im],
        ]
            @test begin
                imat_julia = MSTM.SpecialFunctions.twobytwoinverse(mat)
                imat_fortran = MSTM.Wrapper.twobytwoinverse(mstm, mat)
                isapprox(imat_julia, imat_fortran)
            end
        end
    end

    @testset "Data" begin
        @testset "read MSTM config from .toml file" begin
            config = MSTM.Data.MSTMConfig("fixtures/default.toml")
            @test length(config.spheres) == 2
        end
    end

    @testset "Mie" begin
        ri = OffsetArray(
            [[1.1 + 1.3im, 0.1 + 0.0im] [1.2 + 0.2im, 0.1 + 0.0im] [1.5 - 0.4im, 0.1 + 0.0im] [
                1.5 + 1.5im,
                0.1 + 0.0im,
            ]],
            1:2,
            0:3,
        )

        @testset "mieoa($x, $ri, $nodr0, $qeps, $ri_medium)" for (x, ri, nodr0, qeps, ri_medium) in [
            (0.35, [1.1 + 1.0im, 2.2 - 1.5im], 3, 1e-6, nothing),
            (0.15, [1.1 + 1.0im, 2.2 - 1.5im], 3, 1e-6, [1.0 + 0.1im, 1.0 + 0.1im]),
            (0.25, [1.1 + 1.0im, 2.2 - 1.5im], 4, 0.0, [1.0 + 0.1im, 1.0 + 0.1im]),
            (0.45, [1.1 + 1.0im, 2.2 - 1.5im], 4, -6.0, [1.0 + 0.1im, 1.0 + 0.1im]),
        ]
            @test begin
                qext_julia, qsca_julia, qabs_julia, nodr_julia, _, _, _, _, _, _ =
                    MSTM.Mie.mieoa(x, ri, nodr0, qeps, ri_medium = ri_medium)
                qext_fortran, qsca_fortran, qabs_fortran, nodr_fortran =
                    MSTM.Wrapper.mieoa(mstm, x, ri, nodr0, qeps, ri_medium = ri_medium)
                isapprox(qext_julia, qext_fortran) &&
                    isapprox(qsca_julia, qsca_fortran) &&
                    isapprox(qabs_julia, qabs_fortran) &&
                    nodr_julia == nodr_fortran
            end
        end

        @testset "miecoefcalc($xsp, $hostsphere, $numberfieldexp, $qeps) (No T-Matrix files)" for (
            xsp,
            hostsphere,
            numberfieldexp,
            qeps,
        ) in [
            ([1.2, 0.4, 0.1], [0, 1, 2], [2, 2, 2], 1e-4),
            ([1.2, 0.4, 0.1], [0, 1, 2], [3, 3, 3], 1e-6),
            ([1.2, 0.4, 0.1], [0, 1, 2], [3, 3, 3], 1e-8),
            ([1.2, 0.4, 0.1], [0, 1, 2], [4, 4, 4], 1e-6),
        ]
            @test begin
                params = MSTM.Data.MSTMParameters()
                mie_julia = MSTM.Mie.miecoefcalc(params, xsp, ri, hostsphere, numberfieldexp, qeps)
                mie_fortran = MSTM.Wrapper.miecoefcalc(mstm, xsp, ri, hostsphere, numberfieldexp, qeps)
                isapprox(mie_julia.an, mie_fortran.an) &&
                    isapprox(mie_julia.cn, mie_fortran.cn) &&
                    isapprox(mie_julia.qext, mie_fortran.qext) &&
                    isapprox(mie_julia.qabs, mie_fortran.qabs) &&
                    mie_julia.order == mie_fortran.order &&
                    mie_julia.offset == mie_fortran.offset &&
                    mie_julia.block == mie_fortran.block &&
                    mie_julia.block_offset == mie_fortran.block_offset &&
                    mie_julia.number_field_expansions == mie_fortran.number_field_expansions &&
                    mie_julia.is_optically_active == mie_fortran.is_optically_active &&
                    mie_julia.number_eqns == mie_fortran.number_eqns
            end
        end

        @testset "onemiecoeffmult($i, $mie_coefficent)" for (i, mie_coefficent) in
                                                            [(1, 'a'), (2, 'c'), (3, 'd'), (2, 'u'), (3, 'v')]
            @test begin
                xsp = [1.2, 0.4, 0.1]
                hostsphere = [0, 1, 2]
                numberfieldexp = [3, 3, 3]
                qeps = 1e-6
                cb = 0.5
                E_ϕ = 1.0 + 1.0im
                params = MSTM.Data.MSTMParameters()
                mie = MSTM.Mie.miecoefcalc(params, xsp, ri, hostsphere, numberfieldexp, qeps)
                nodr = mie.order[i]
                ctx = MSTM.Constants.init()
                MSTM.Wrapper.miecoefcalc(mstm, xsp, ri, hostsphere, numberfieldexp, qeps)
                cx = MSTM.SpecialFunctions.pifunc(ctx, cb, E_ϕ, nodr, nodr)
                cy_julia = MSTM.Mie.onemiecoeffmult(mie, i, cx, mie_coefficent)
                cy_fortran = MSTM.Wrapper.onemiecoeffmult(mstm, i, nodr, cx, mie_coefficent)
                isapprox(cy_julia, cy_fortran)
            end
        end

        @testset "multmiecoeffmult" begin
            xsp = [1.2, 0.4, 0.1]
            hostsphere = [0, 1, 2]
            numberfieldexp = [3, 2, 1]
            qeps = 1e-6
            cb = 0.5
            E_ϕ = 1.0 + 1.0im
            params = MSTM.Data.MSTMParameters()
            mie = MSTM.Mie.miecoefcalc(params, xsp, ri, hostsphere, numberfieldexp, qeps)
            MSTM.Wrapper.miecoefcalc(mstm, xsp, ri, hostsphere, numberfieldexp, qeps)
            for i in 1:4
                @test begin
                    nrhs = i + 3
                    ain = 1.0im * rand(mie.number_eqns * nrhs) + rand(mie.number_eqns * nrhs)
                    aout_julia = MSTM.Mie.multmiecoeffmult(mie, nrhs, i % 2, ain)
                    aout_fortran = MSTM.Wrapper.multmiecoeffmult(mstm, mie.number_eqns, nrhs, i % 2, ain)
                    isapprox(aout_julia, aout_fortran)
                end
            end
        end
    end

    Libdl.dlclose(mstm)
end


text:
let b1 = (1, 2)
    begin
        function f1(a1 :: Tuple{A, B}) where {A, B}
            b2 = a1[1]
            function f2(a2 :: A)
                A
                b3 = a1[2]
                function f3(a3 :: B)
                   # B # if don't explicitly give a `B` here, err can be raised
                   a2 + a3
                end
                f3(b3)
            end
            f2(b2)
        end
        f1(b1)
    end
end


text:
#
# nonseparable.jl --
#
# Non-separable multidimensional interpolation.
#
#------------------------------------------------------------------------------
#
# This file is part of the LinearInterpolators package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2016-2021, Éric Thiébaut.
#

# FIXME: if axes are aligned, use separable interpolation.

"""
    TwoDimensionalTransformInterpolator(rows, cols, ker1, ker2, R)

yields a linear mapping which interpolates its input of size `cols` to produce
an output of size `rows` by 2-dimensional interpolation with kernels `ker1` and
`ker2` along each dimension and applying the affine coordinate transform
specified by `R`.

As a shortcut:

    TwoDimensionalTransformInterpolator(rows, cols, ker, R)

is equivalent to `TwoDimensionalTransformInterpolator(rows,cols,ker,ker,R)`
that is the same kernel is used along all dimensions.

""" TwoDimensionalTransformInterpolator

struct TwoDimensionalTransformInterpolator{T<:AbstractFloat,
                                           K1<:Kernel{T},
                                           K2<:Kernel{T}} <: LinearMapping
    rows::Dims{2}
    cols::Dims{2}
    ker1::K1
    ker2::K2
    R::AffineTransform2D{T}
end

function TwoDimensionalTransformInterpolator(rows::Dims{2},
                                             cols::Dims{2},
                                             ker::Kernel, R::AffineTransform2D)
    TwoDimensionalTransformInterpolator(rows, cols, ker, ker, R)
end

function TwoDimensionalTransformInterpolator(rows::Dims{2},
                                             cols::Dims{2},
                                             ker1::Kernel{T1},
                                             ker2::Kernel{T2},
                                             R::AffineTransform2D{Tr}) where {T1,T2,Tr}
    T = promote_type(T1, T2, Tr)
    TwoDimensionalTransformInterpolator(rows, cols,
                                        convert(Kernel{T}, ker1),
                                        convert(Kernel{T}, ker2),
                                        convert(AffineTransform2D{T}, R))
end

# TODO: Replace assertions by error messages.
input_size(A::TwoDimensionalTransformInterpolator) = A.cols
output_size(A::TwoDimensionalTransformInterpolator) = A.rows

function vcreate(::Type{Direct}, A::TwoDimensionalTransformInterpolator{T},
                 x::AbstractArray{T,2}, scratch::Bool = false) where {T}
    # Checking x is done by apply!
    Array{T,2}(undef, A.rows)
end

function vcreate(::Type{Adjoint}, A::TwoDimensionalTransformInterpolator{T},
                 x::AbstractArray{T,2}, scratch::Bool = false) where {T}
    # Checking x is done by apply!
    Array{T,2}(undef, A.cols)
end

function apply!(α::Real,
                ::Type{Direct},
                A::TwoDimensionalTransformInterpolator{T},
                x::AbstractArray{T,2},
                scratch::Bool,
                β::Real,
                y::AbstractArray{T,2},) where {T}
    @assert !Base.has_offset_axes(x, y)
    @assert size(x) == A.cols
    @assert size(y) == A.rows
    apply!(α, Direct, A.ker1, A.ker2, A.R, x, β, y)
end

function apply!(α::Real,
                ::Type{Adjoint},
                A::TwoDimensionalTransformInterpolator{T},
                x::AbstractArray{T,2},
                scratch::Bool,
                β::Real,
                y::AbstractArray{T,2},) where {T}
    @assert !Base.has_offset_axes(x, y)
    @assert size(x) == A.rows
    @assert size(y) == A.cols
    apply!(α, Adjoint, A.ker1, A.ker2, A.R, x, β, y)
end

# Provide default Direct operation.
function apply!(dst::AbstractArray{T,2},
                ker::Kernel{T},
                R::AffineTransform2D{T},
                src::AbstractArray{T,2}) where {T}
    apply!(dst, Direct, ker, R, src)
end
function apply!(dst::AbstractArray{T,2},
                ker1::Kernel{T},
                ker2::Kernel{T},
                R::AffineTransform2D{T},
                src::AbstractArray{T,2}) where {T}
    return apply!(dst, Direct, ker1, ker2, R, src)
end

# Provide default α=1 and β=0 factors.
function apply!(dst::AbstractArray{T,2},
                ::Type{P},
                ker::Kernel{T},
                R::AffineTransform2D{T},
                src::AbstractArray{T,2}) where {P<:Operations,T}
    return apply!(1, P, ker, R, src, 0, dst)
end
function apply!(dst::AbstractArray{T,2},
                ::Type{P},
                ker1::Kernel{T},
                ker2::Kernel{T},
                R::AffineTransform2D{T},
                src::AbstractArray{T,2}) where {P<:Operations,T}
    return apply!(1, P, ker1, ker2, R, src, 0, dst)
end

# Provide default pair of kernels (ker1,ker2) = ker.
function apply!(α::Real,
                ::Type{P},
                ker::Kernel{T},
                R::AffineTransform2D{T},
                src::AbstractArray{T,2},
                β::Real,
                dst::AbstractArray{T,2}) where {P<:Operations,T}
    return apply!(α, P, ker, ker, R, src, β, dst)
end

@generated function apply!(α::Real,
                           ::Type{Direct},
                           ker1::Kernel{T,S1},
                           ker2::Kernel{T,S2},
                           R::AffineTransform2D{T},
                           src::AbstractArray{T,2},
                           β::Real,
                           dst::AbstractArray{T,2}) where {T<:AbstractFloat,
                                                           S1,S2}
    # Generate pieces of code.
    J1 = [Symbol(:j1_,s) for s in 1:S1]
    W1 = [Symbol(:w1_,s) for s in 1:S1]
    J2 = [Symbol(:j2_,s) for s in 1:S2]
    W2 = [Symbol(:w2_,s) for s in 1:S2]
    code2 = (:( pos2 = convert(T, i2)  ),
             :( off1 = R.xy*pos2 + R.x ),
             :( off2 = R.yy*pos2 + R.y ))
    code1 = (:( pos1 = convert(T, i1)  ),
             Meta.generate_getcoefs(J1, W1, :ker1, :lim1, :(R.xx*pos1 + off1)),
             Meta.generate_getcoefs(J2, W2, :ker2, :lim2, :(R.yx*pos1 + off2)))
    expr = Meta.generate_interp_expr(:src, J1, W1, J2, W2)

    quote
        if α == 0
            # Just scale destination.
            vscale!(dst, β)
        else
            # Get dimensions and limits.
            n1, n2 = size(dst)
            lim1 = limits(ker1, size(src, 1))
            lim2 = limits(ker2, size(src, 2))

            # Apply the operator considering the specific values of α and β.
            if α == 1 && β == 0
                for i2 in 1:n2
                    $(code2...)
                    @inbounds for i1 in 1:n1
                        $(code1...)
                        dst[i1,i2] = $expr
                    end
                end
            else
                alpha = convert(T, α)
                beta = convert(T, β)
                for i2 in 1:n2
                    $(code2...)
                    @inbounds for i1 in 1:n1
                        $(code1...)
                        dst[i1,i2] = $expr*alpha + dst[i1,i2]*beta
                    end
                end
            end
        end
        return dst
    end
end

@generated function apply!(α::Real,
                           ::Type{Adjoint},
                           ker1::Kernel{T,S1},
                           ker2::Kernel{T,S2},
                           R::AffineTransform2D{T},
                           src::AbstractArray{T,2},
                           β::Real,
                           dst::AbstractArray{T,2}) where{T<:AbstractFloat,
                                                          S1,S2}
    # Generate pieces of code.
    J1 = [Symbol(:j1_,s) for s in 1:S1]
    W1 = [Symbol(:w1_,s) for s in 1:S1]
    J2 = [Symbol(:j2_,s) for s in 1:S2]
    W2 = [Symbol(:w2_,s) for s in 1:S2]
    temp = [Symbol(:tmp_,s) for s in 1:S2]
    code2 = (:( pos2 = convert(T, i2)  ),
             :( off1 = R.xy*pos2 + R.x ),
             :( off2 = R.yy*pos2 + R.y ))
    code1 = [:( pos1 = convert(T, i1)  ),
             Meta.generate_getcoefs(J1, W1, :ker1, :lim1, :(R.xx*pos1 + off1)),
             Meta.generate_getcoefs(J2, W2, :ker2, :lim2, :(R.yx*pos1 + off2))]
    for i2 in 1:S2
        push!(code1, :( $(temp[i2]) = $(W2[i2])*val ))
        for i1 in 1:S1
            push!(code1, :(
                dst[$(J1[i1]),$(J2[i2])] += $(W1[i1])*$(temp[i2])
            ))
        end
    end

    quote
        # Pres-scale or zero destination.
        vscale!(dst, β)

        # Get dimensions and limits.
        n1, n2 = size(src)
        lim1 = limits(ker1, size(dst, 1))
        lim2 = limits(ker2, size(dst, 2))

        # Apply adjoint operator.
        if α == 1
            for i2 in 1:n2
                $(code2...)
                @inbounds for i1 in 1:n1
                    val = src[i1,i2]
                    $(code1...)
                end
            end
        elseif α != 0
            alpha = convert(T, α)
            for i2 in 1:n2
                $(code2...)
                @inbounds for i1 in 1:n1
                    val = alpha*src[i1,i2]
                    $(code1...)
                end
            end
        end
        return dst
    end
end


text:
abstract type TrilinosSolver end

"""
Solver that works directly on Tpetra CrsMatrix and MultiVector, parametrized using a Teuchos ParameterList
"""
struct TpetraSolver{SolverT} <: TrilinosSolver
  solver::SolverT
  parameters::CxxWrap.SmartPointer{Teuchos.ParameterList}
  TpetraSolver{SolverT}(solver::SolverT,params::CxxWrap.SmartPointer{Teuchos.ParameterList}) where SolverT = new(solver,params)
end

function TpetraSolver(A::CxxUnion{Tpetra.CrsMatrix{ST,LT,GT,NT}}, parameterlist = Teuchos.ParameterList()) where {ST,LT,GT,NT}
  params = Teuchos.ParameterListPair(parameterlist, Teuchos.ParameterList())
  linear_solver_type = get(params, "Linear Solver Type", "Belos")
  solver_types = Teuchos.sublist(params, "Linear Solver Types")

  if linear_solver_type == "Belos"
    factory = Belos.SolverFactory{ST, Tpetra.MultiVector{ST,LT,GT,NT}, Tpetra.Operator{ST,LT,GT,NT}}()
    belos_parameters = Teuchos.sublist(solver_types, "Belos")
    solver_name = get(belos_parameters, "Solver Type", "BLOCK GMRES")
    belos_solver_types = Teuchos.sublist(belos_parameters, "Solver Types")
    solver_parameters = Teuchos.ParameterList(Teuchos.sublist(belos_solver_types.input, solver_name))

    solver = Belos.create(factory, solver_name, solver_parameters)
    linprob = Belos.LinearProblem(A)

    prec_type = get(params, "Preconditioner Type", "Ifpack2")
    if prec_type == "Ifpack2"
      prec_list = Teuchos.sublist(params, "Preconditioner Types")
      ifpack_list = Teuchos.sublist(prec_list, "Ifpack2")
      ifpack2_type = get(ifpack_list, "Preconditioner Type", "ILUT")
      ifpack2_params = Teuchos.sublist(ifpack_list, ifpack2_type)
      ifpack2_params_new = Teuchos.ParameterList(Teuchos.sublist(ifpack_list.input, ifpack2_type))
      prec_factory = Ifpack2.Factory()
      M = Ifpack2.create(prec_factory, ifpack2_type, A)
      Ifpack2.setParameters(M, ifpack2_params_new)
      Ifpack2.initialize(M)
      Ifpack2.compute(M)
      Belos.setRightPrec(linprob, M)
      ifpack_list.output[ifpack2_type] = ifpack2_params_new
    elseif prec_type == "MueLu"
      prec_list = Teuchos.sublist(params, "Preconditioner Types")
      muelu_list = Teuchos.sublist(prec_list, "MueLu")
      if isempty(muelu_list.input)
        MueLu.parameters_easy!(muelu_list.output)
      end

      muelu_prec = MueLu.CreateTpetraPreconditioner(A, muelu_list.output)
      Belos.setRightPrec(linprob, muelu_prec)
    elseif prec_type == "None"
      nothing
    else
      error("Unsupported  preconditioner type $prec_type")
    end

    Belos.setProblem(solver,linprob)

    belos_solver_types.output[solver_name] = solver_parameters
    return TpetraSolver{typeof(solver)}(solver, params.output)
  end
  error("Unsupported linear solver type: $linear_solver_type")
end

function default_parameters(solver_type="BLOCK GMRES"; prec_factory = Ifpack2.default_parameters)
  params = Dict("Linear Solver Type" => "Belos", "Linear Solver Types" => Dict("Belos" => Dict("Solver Type" => solver_type, "Solver Types" => Teuchos.ParameterList())))
  params["Preconditioner Type"] = "Ifpack2"

  prec_types = Teuchos.ParameterList("Preconditioner Types")
  prec_params = prec_factory()
  prec_types[Teuchos.name(prec_params)] = prec_params
  params["Preconditioner Types"] = prec_types

  stypes = params["Linear Solver Types"]["Belos"]["Solver Types"]
  sp = Teuchos.ParameterList()
  factory = Belos.SolverFactory{Float64, Tpetra.MultiVector{Float64,Int32,Int64,Kokkos.default_node_type()}, Tpetra.Operator{Float64,Int32,Int64,Kokkos.default_node_type()}}()
  solver = Belos.create(factory, solver_type, sp)
  stypes[solver_type] = sp
  finalize(solver)
  finalize(factory)
  return params
end

import Base: \

function \(A::TrilinosSolver, b::CxxUnion{Tpetra.Vector{ST,LT,GT,NT}}) where {ST,LT,GT,NT}
  prob = Belos.getProblem(A.solver)
  x = Tpetra.Vector(Tpetra.getDomainMap(Belos.getOperator(prob)))
  Belos.setProblem(prob,x,b)
  Belos.solve(A.solver)
  return x
end


text:
# TODO I think here is some potential for faster computation
# However, up to 25 states, the algorithm works fine and fast (main knobs are rtol and maxiter)
# This is the specialized version assuming a mass matrix form / linear in dx
# M(x, p)*dx = f(x, p)
# Where M is diagonal!

# TODO preallocation

"""
    ISInDy(X, Y, Ψ; weights, f_target, maxiter, rtol, p, t, opt)

Performs an implicit sparse identification of nonlinear dynamics given the data matrices `X` and `Y` via the `AbstractBasis` `basis.`
Keyworded arguments include the parameter (values) of the basis `p` and the timepoints `t`, which are passed in optionally.
`opt` is an `AbstractSubspaceOptimizer` useable for sparse regression within the nullspace, `maxiter` the maximum iterations to perform, and `convergence_error` the
bound which causes the optimizer to stop.

The best vectors of the sparse nullspace are selected via multi-objective optimization.
The best candidate is determined via the `AbstractScalarizationMethod` given in `alg`. The evaluation has two fields, the sparsity of the coefficients and the L2-Norm error at index 1 and 2, respectively.

`rtol` gets directly passed into the computation of the nullspace.

Returns a `SInDyResult`.
"""
function ISInDy(X::AbstractArray, Ẋ::AbstractArray, Ψ::Basis; alg::Optimize.AbstractScalarizationMethod = WeightedSum(), maxiter::Int64 = 10, rtol::Float64 = 0.99, p::AbstractArray = [], t::AbstractVector = [], opt::T = ADM()) where T <: DataDrivenDiffEq.Optimize.AbstractSubspaceOptimizer
    @assert size(X)[end] == size(Ẋ)[end]

    # Compute the library and the corresponding nullspace
    θ = Ψ(X, p, t)

    # Init for sweep over the differential variables
    Ξ = zeros(eltype(θ), length(Ψ)*2, size(Ẋ, 1))

    opt_front = ParetoFront(1, scalarization = alg)
    tmp_front = ParetoFront(1, scalarization = alg)

    @inbounds for i in 1:size(Ẋ, 1)
        dθ = hcat(map((dxi, ti)->dxi.*ti, Ẋ[i, :], eachcol(θ))...)
        Θ = vcat(dθ, θ)
        N = nullspace(Θ', rtol = rtol)
        Q = deepcopy(N) # Deepcopy for inplace

        # Find sparse vectors in nullspace
        # Calls effectively the ADM algorithm with varying initial conditions
        DataDrivenDiffEq.fit!(Q, N', opt, maxiter = maxiter)

        # Compute pareto front
        @inbounds for (i, qi) in enumerate(eachcol(Q))
            if i == 1
                set_candidate!(opt_front, 1, [norm(qi, 0); norm(Θ'*qi, 2)], qi, maxiter, get_threshold(opt))
                set_candidate!(tmp_front, 1, [norm(qi, 0); norm(Θ'*qi, 2)], qi, maxiter, get_threshold(opt))
            else
                set_candidate!(tmp_front, 1, [norm(qi, 0); norm(Θ'*qi, 2)], qi, maxiter, get_threshold(opt))
                conditional_add!(opt_front, tmp_front)
            end
        end
        Ξ[:, i] .= Optimize.parameter(opt_front[1])
        Ξ[:, i] .= Ξ[:, i] ./ maximum(abs.(Ξ[:, i]))
    end

    return ImplicitSparseIdentificationResult(Ξ, Ψ, maxiter, opt, true, Ẋ, X, p = p, t = t)
end

function ImplicitSparseIdentificationResult(coeff::AbstractArray, equations::Basis, iters::Int64, opt::T, convergence::Bool, Y::AbstractVecOrMat, X::AbstractVecOrMat; p::AbstractArray = [], t::AbstractVector = []) where T <: Union{Optimize.AbstractOptimizer, Optimize.AbstractSubspaceOptimizer}

    sparsities = Int64.(norm.(eachcol(coeff), 0))

    b_, p_ = derive_implicit_parameterized_eqs(coeff, equations)
    ps = [p; p_]

    Ŷ = b_(X, ps, t)
    training_error = norm.(eachrow(Y-Ŷ), 2)
    aicc = similar(training_error)

    for i in 1:length(aicc)
        aicc[i] = AICC(sum(sparsities[i]), view(Ŷ, i, :) , view(Y, i, :))
    end
    return SparseIdentificationResult(coeff, [p...;p_...], b_ , opt, iters, convergence,  training_error, aicc,  sparsities)
end

function derive_implicit_parameterized_eqs(Ξ::AbstractArray{T, 2}, b::Basis) where T <: Real

    sparsity = Int64(norm(Ξ, 0))

    @parameters p[1:sparsity]
    p_ = zeros(eltype(Ξ), sparsity)
    cnt = 1

    b_ = Basis(Operation[], variables(b), parameters = [parameters(b)...; p...])

    for i=1:size(Ξ, 2)
        eq_d = nothing
        eq_n = nothing
        # Denominator
        for j = 1:length(b)
            if !iszero(Ξ[j,i])
                if eq_d === nothing
                    eq_d = p[cnt]*b[j]
                else
                    eq_d += p[cnt]*b[j]
                end
                p_[cnt] = Ξ[j,i]
                cnt += 1
            end
        end
        # Numerator
        for j = 1:length(b)
            if !iszero(Ξ[j+length(b),i])
                if eq_n === nothing
                    eq_n = p[cnt]*b[j]
                else
                    eq_n += p[cnt]*b[j]
                end
                p_[cnt] = Ξ[j+length(b),i]
                cnt += 1
            end
        end

        push!(b_, -eq_n ./ eq_d)
    end
    b_, p_
end