text:
# Wrapper for the low level functions defined in https://github.com/oxfordcontrol/osqp/blob/master/include/osqp.h

# Ensure compatibility between Julia versions with @gc_preserve
macro compat_gc_preserve(args...)
    vars = args[1:end - 1]
    body = args[end]
    if VERSION > v"0.7.0-"
        return esc(Expr(:macrocall, Expr(:., :Base, Base.Meta.quot(Symbol("@gc_preserve"))), __source__, args...))
    else
        return esc(body)
    end
end

mutable struct Model
    workspace::Ptr{OSQP.Workspace}
    lcache::Vector{Float64} # to facilitate converting l to use OSQP_INFTY
    ucache::Vector{Float64} # to facilitate converting u to use OSQP_INFTY

    """
        Module()

    Initialize OSQP module
    """
    function Model()
        model = new(C_NULL, Float64[], Float64[])
        @compat finalizer(OSQP.clean!, model)
        return model

    end

end

"""
    setup!(model, P, q, A, l, u, settings)

Perform OSQP solver setup of model `model`, using the inputs `P`, `q`, `A`, `l`, `u`.
"""
function setup!(model::OSQP.Model;
        P::Union{SparseMatrixCSC,Nothing} = nothing,
        q::Union{Vector{Float64},Nothing} = nothing,
        A::Union{SparseMatrixCSC,Nothing} = nothing,
        l::Union{Vector{Float64},Nothing} = nothing,
        u::Union{Vector{Float64},Nothing} = nothing,
        settings...)

    # Check problem dimensions
    if P == nothing
        if q != nothing
            n = length(q)
        elseif A != nothing
            n = size(A, 2)
        else
            error("The problem does not have any variables!")
        end

    else
        n = size(P, 1)
    end

    if A == nothing
        m = 0
    else
        m = size(A, 1)
    end

    # Check if parameters are nothing
    if ((A == nothing) & ( (l != nothing) | (u != nothing))) |
        ((A != nothing) & ((l == nothing) | (u == nothing)))
        error("A must be supplied together with l and u")
    end

    if (A != nothing) & (l == nothing)
        l = -Inf * ones(m)
    end
    if (A != nothing) & (u == nothing)
        u = Inf * ones(m)
    end

    if P == nothing
        P = sparse([], [], [], n, n)
    end
    if q == nothing
        q = zeros(n)
    end
    if A == nothing
        A = sparse([], [], [], m, n)
        l = zeros(m)
        u = zeros(m)
    end

    # Check if dimensions are correct
    if length(q) != n
        error("Incorrect dimension of q")
    end
    if length(l) != m
        error("Incorrect dimensions of l")
    end
    if length(u) != m
        error("Incorrect dimensions of u")
    end

    # Check or sparsify matrices
    if !issparse(P)
        warn("P is not sparse. Sparsifying it now (it might take a while)")
        P = sparse(P)
    end
    if !issparse(A)
        warn("A is not sparse. Sparsifying it now (it might take a while)")
        A = sparse(A)
    end

    # Convert lower and upper bounds from Julia infinity to OSQP infinity
    u = min.(u, OSQP_INFTY)
    l = max.(l, -OSQP_INFTY)

    # Resize caches
    resize!(model.lcache, m)
    resize!(model.ucache, m)

    # Create managed matrices to avoid segfaults (See SCS.jl)
    managedP = OSQP.ManagedCcsc(P)
    managedA = OSQP.ManagedCcsc(A)

    # Get managed pointers (Ref) Pdata and Adata
    Pdata = Ref(OSQP.Ccsc(managedP))
    Adata = Ref(OSQP.Ccsc(managedA))

    # Create OSQP data using the managed matrices pointers
    data = OSQP.Data(n, m,
                     Base.unsafe_convert(Ptr{OSQP.Ccsc}, Pdata),
                     Base.unsafe_convert(Ptr{OSQP.Ccsc}, Adata),
                     pointer(q),
                     pointer(l), pointer(u))

    # Create OSQP settings
    settings_dict = Dict{Symbol,Any}()
    if !isempty(settings)
        for (key, value) in settings
            settings_dict[key] = value
        end
    end

    stgs = OSQP.Settings(settings_dict)

    # Perform setup
    @compat_gc_preserve managedP Pdata managedA Adata q l u begin
        model.workspace = ccall((:osqp_setup, OSQP.osqp), Ptr{OSQP.Workspace},
            (Ptr{OSQP.Data}, Ptr{OSQP.Settings}), Ref(data), Ref(stgs))
    end

    if model.workspace == C_NULL
        error("Error in OSQP setup")
    end

end

function solve!(model::OSQP.Model, results::Results = Results())
    ccall((:osqp_solve, OSQP.osqp), Cc_int,
             (Ptr{OSQP.Workspace}, ),
             model.workspace)
    workspace = unsafe_load(model.workspace)
    info = results.info
    Compat.copyto!(info, unsafe_load(workspace.info))
    solution = unsafe_load(workspace.solution)
    data = unsafe_load(workspace.data)
    n = data.n
    m = data.m
    resize!(results, n, m)
    has_solution = false
    for status in SOLUTION_PRESENT
        info.status == status && (has_solution = true; break)
    end
    if has_solution
        # If solution exists, copy it
        unsafe_copyto!(pointer(results.x), solution.x, n)
        unsafe_copyto!(pointer(results.y), solution.y, m)
        fill!(results.prim_inf_cert, NaN)
        fill!(results.dual_inf_cert, NaN)
    else
        # else fill with NaN and return certificates of infeasibility
        fill!(results.x, NaN)
        fill!(results.y, NaN)
        if info.status == :Primal_infeasible || info.status == :Primal_infeasible_inaccurate
            unsafe_copyto!(pointer(results.prim_inf_cert), workspace.delta_y, m)
            fill!(results.dual_inf_cert, NaN)
        elseif info.status == :Dual_infeasible || info.status == :Dual_infeasible_inaccurate
            fill!(results.prim_inf_cert, NaN)
            unsafe_copyto!(pointer(results.dual_inf_cert), workspace.delta_x, n)
        else
            fill!(results.prim_inf_cert, NaN)
            fill!(results.dual_inf_cert, NaN)
        end
    end

    results
end

function version()
    return unsafe_string(ccall((:osqp_version, OSQP.osqp), Cstring, ()))
end

function clean!(model::OSQP.Model)
    exitflag = ccall((:osqp_cleanup, OSQP.osqp), Cc_int,
             (Ptr{OSQP.Workspace},), model.workspace)
    if exitflag != 0
        error("Error in OSQP cleanup")
    end
end

function update_q!(model::OSQP.Model, q::Vector{Float64})
    (n, m) = OSQP.dimensions(model)
    if length(q) != n
        error("q must have length n = $(n)")
    end
    exitflag = ccall((:osqp_update_lin_cost, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}), model.workspace, q)
    if exitflag != 0 error("Error updating q") end
end

function update_l!(model::OSQP.Model, l::Vector{Float64})
    (n, m) = OSQP.dimensions(model)
    if length(l) != m
        error("l must have length m = $(m)")
    end
    model.lcache .= max.(l, -OSQP_INFTY)
    exitflag = ccall((:osqp_update_lower_bound, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}), model.workspace, model.lcache)
    if exitflag != 0 error("Error updating l") end
end

function update_u!(model::OSQP.Model, u::Vector{Float64})
    (n, m) = OSQP.dimensions(model)
    if length(u) != m
        error("u must have length m = $(m)")
    end
    model.ucache .= min.(u, OSQP_INFTY)
    exitflag = ccall((:osqp_update_upper_bound, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}), model.workspace, model.ucache)
    if exitflag != 0 error("Error updating u") end
end

function update_bounds!(model::OSQP.Model, l::Vector{Float64}, u::Vector{Float64})
    (n, m) = OSQP.dimensions(model)
    if length(l) != m
        error("l must have length m = $(m)")
    end
    if length(u) != m
        error("u must have length m = $(m)")
    end
    model.lcache .= max.(l, -OSQP_INFTY)
    model.ucache .= min.(u, OSQP_INFTY)
    exitflag = ccall((:osqp_update_bounds, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}, Ptr{Cdouble}),
        model.workspace, model.lcache, model.ucache)
    if exitflag != 0 error("Error updating bounds l and u") end
end

prep_idx_vector_for_ccall(idx::Nothing, n::Int, namesym::Symbol) = C_NULL
function prep_idx_vector_for_ccall(idx::Vector{Int}, n::Int, namesym::Symbol)
    if length(idx) != n
        error("$(namesym) and $(namesym)_idx must have the same length")
    end
    idx .-= 1 # Shift indexing to match C
    idx
end

restore_idx_vector_after_ccall!(idx::Nothing) = nothing
function restore_idx_vector_after_ccall!(idx::Vector{Int})
    idx .+= 1 # Unshift indexing
    nothing
end

function update_P!(model::OSQP.Model, Px::Vector{Float64}, Px_idx::Union{Vector{Int}, Nothing})
    Px_idx_prepped = prep_idx_vector_for_ccall(Px_idx, length(Px), :P)
    exitflag = ccall((:osqp_update_P, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}, Ptr{Cc_int}, Cc_int),
        model.workspace, Px, Px_idx_prepped, length(Px))
    restore_idx_vector_after_ccall!(Px_idx)
    if exitflag != 0 error("Error updating P") end
end

function update_A!(model::OSQP.Model, Ax::Vector{Float64}, Ax_idx::Union{Vector{Int}, Nothing})
    Ax_idx_prepped = prep_idx_vector_for_ccall(Ax_idx, length(Ax), :A)
    exitflag = ccall((:osqp_update_A, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}, Ptr{Cc_int}, Cc_int),
        model.workspace, Ax, Ax_idx_prepped, length(Ax))
    restore_idx_vector_after_ccall!(Ax_idx)
    if exitflag != 0 error("Error updating A") end
end

function update_P_A!(model::OSQP.Model, Px::Vector{Float64}, Px_idx::Union{Vector{Int}, Nothing}, Ax::Vector{Float64}, Ax_idx::Union{Vector{Int}, Nothing})
    Px_idx_prepped = prep_idx_vector_for_ccall(Px_idx, length(Px), :P)
    Ax_idx_prepped = prep_idx_vector_for_ccall(Ax_idx, length(Ax), :A)
    exitflag = ccall((:osqp_update_P_A, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble},
        Ptr{Cc_int}, Cc_int, Ptr{Cdouble}, Ptr{Cc_int}, Cc_int),
        model.workspace, Px, Px_idx_prepped, length(Px), Ax, Ax_idx_prepped, length(Ax))
    restore_idx_vector_after_ccall!(Ax_idx)
    restore_idx_vector_after_ccall!(Px_idx)
    if exitflag != 0 error("Error updating P and A") end
end

function update!(model::OSQP.Model; q = nothing, l = nothing, u = nothing, Px = nothing, Px_idx = nothing, Ax = nothing, Ax_idx = nothing)
    # q
    if q != nothing
        update_q!(model, q)
    end

    # l and u
    if l != nothing && u != nothing
        update_bounds!(model, l, u)
    elseif l != nothing
        update_l!(model, l)
    elseif u != nothing
        update_u!(model, u)
    end

    # P and A
    if Px != nothing && Ax != nothing
        update_P_A!(model, Px, Px_idx, Ax, Ax_idx)
    elseif Px != nothing
        update_P!(model, Px, Px_idx)
    elseif Ax != nothing
        update_A!(model, Ax, Ax_idx)
    end
end

function update_settings!(model::OSQP.Model; kwargs...)

    if isempty(kwargs)
        return
    else
        data = Dict{Symbol,Any}()
        for (key, value) in kwargs
            if !(key in UPDATABLE_SETTINGS)
                error("$(key) cannot be updated or is not recognized")
            else
                data[key] = value
            end
        end
    end

    # Get arguments
    max_iter = get(data, :max_iter, nothing)
    eps_abs = get(data, :eps_abs, nothing)
    eps_rel = get(data, :eps_rel, nothing)
    eps_prim_inf = get(data, :eps_prim_inf, nothing)
    eps_dual_inf = get(data, :eps_dual_inf, nothing)
    rho = get(data, :rho, nothing)
    alpha = get(data, :alpha, nothing)
    delta = get(data, :delta, nothing)
    polish = get(data, :polish, nothing)
    polish_refine_iter = get(data, :polish_refine_iter, nothing)
    verbose = get(data, :verbose, nothing)
    scaled_termination = get(data, :early_terminate, nothing)
    check_termination = get(data, :check_termination, nothing)
    warm_start = get(data, :warm_start, nothing)
    time_limit = get(data, :time_limit, nothing)

    # Update individual settings
    if max_iter != nothing
        exitflag = ccall((:osqp_update_max_iter, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cc_int), model.workspace, max_iter)
        if exitflag != 0 error("Error updating max_iter") end
    end

    if eps_abs != nothing
        exitflag = ccall((:osqp_update_eps_abs, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cdouble), model.workspace, eps_abs)
        if exitflag != 0 error("Error updating eps_abs") end
    end

    if eps_rel != nothing
        exitflag = ccall((:osqp_update_eps_rel, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cdouble), model.workspace, eps_rel)
        if exitflag != 0 error("Error updating eps_rel") end
    end

    if eps_prim_inf != nothing
        exitflag = ccall((:osqp_update_eps_prim_inf, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cdouble), model.workspace, eps_prim_inf)
        if exitflag != 0 error("Error updating eps_prim_inf") end
    end

    if eps_dual_inf != nothing
        exitflag = ccall((:osqp_update_eps_dual_inf, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cdouble), model.workspace, eps_dual_inf)
        if exitflag != 0 error("Error updating eps_dual_inf") end
    end

    if rho != nothing
        exitflag = ccall((:osqp_update_rho, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cdouble), model.workspace, rho)
        if exitflag != 0 error("Error updating rho") end
    end

    if alpha != nothing
        exitflag = ccall((:osqp_update_alpha, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cdouble), model.workspace, alpha)
        if exitflag != 0 error("Error updating alpha") end
    end

    if delta != nothing
        exitflag = ccall((:osqp_update_delta, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cdouble), model.workspace, delta)
        if exitflag != 0 error("Error updating delta") end
    end

    if polish != nothing
        exitflag = ccall((:osqp_update_polish, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cc_int), model.workspace, polish)
        if exitflag != 0 error("Error updating polish") end
    end

    if polish_refine_iter != nothing
        exitflag = ccall((:osqp_update_polish_refine_iter, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cc_int), model.workspace, polish_refine_iter)
        if exitflag != 0 error("Error updating polish_refine_iter") end
    end

    if verbose != nothing
        exitflag = ccall((:osqp_update_verbose, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cc_int), model.workspace, verbose)
        if exitflag != 0 error("Error updating verbose") end
    end

    if scaled_termination != nothing
        exitflag = ccall((:osqp_update_scaled_termination, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cc_int), model.workspace, scaled_termination)
        if exitflag != 0 error("Error updating scaled_termination") end
    end

    if check_termination != nothing
        exitflag = ccall((:osqp_update_check_termination, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cc_int), model.workspace, check_termination)
        if exitflag != 0 error("Error updating check_termination") end
    end

    if warm_start != nothing
        exitflag = ccall((:osqp_update_warm_start, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cc_int), model.workspace, warm_start)
        if exitflag != 0 error("Error updating warm_start") end
    end

   if time_limit != nothing
        exitflag = ccall((:osqp_update_time_limit, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Cdouble), model.workspace, time_limit)
        if exitflag != 0 error("Error updating time_limit") end
    end

    return nothing
end

function warm_start_x!(model::OSQP.Model, x::Vector{Float64})
    (n, m) = OSQP.dimensions(model)
    length(x) == n || error("Wrong dimension for variable x")
    exitflag = ccall((:osqp_warm_start_x, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}), model.workspace, x)
    exitflag == 0  || error("Error in warm starting x")
    nothing
end

function warm_start_y!(model::OSQP.Model, y::Vector{Float64})
    (n, m) = OSQP.dimensions(model)
    length(y) == m || error("Wrong dimension for variable y")
    exitflag = ccall((:osqp_warm_start_y, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}), model.workspace, y)
    exitflag == 0 || error("Error in warm starting y")
    nothing
end

function warm_start_x_y!(model::OSQP.Model, x::Vector{Float64}, y::Vector{Float64})
    (n, m) = OSQP.dimensions(model)
    length(x) == n || error("Wrong dimension for variable x")
    length(y) == m || error("Wrong dimension for variable y")
    exitflag = ccall((:osqp_warm_start, OSQP.osqp), Cc_int, (Ptr{OSQP.Workspace}, Ptr{Cdouble}, Ptr{Cdouble}), model.workspace, x, y)
    exitflag == 0 || error("Error in warm starting x and y")
    nothing
end

function warm_start!(model::OSQP.Model; x::Union{Vector{Float64}, Nothing} = nothing, y::Union{Vector{Float64}, Nothing} = nothing)
    if x isa Vector{Float64} && y isa Vector{Float64}
        warm_start_x_y!(model, x, y)
    elseif x isa Vector{Float64}
        warm_start_x!(model, x)
    elseif y isa Vector{Float64}
        warm_start_y!(model, y)
    end
end

# Auxiliary low-level functions
"""
    dimensions(model::OSQP.Model)

Obtain problem dimensions from OSQP model
"""
function dimensions(model::OSQP.Model)

    workspace = unsafe_load(model.workspace)
    if workspace == C_NULL
        error("Workspace has not been setup yet")
    end
    data = unsafe_load(workspace.data)
    return data.n, data.m
end

function linsys_solver_str_to_int!(settings_dict::Dict{Symbol,Any})
         # linsys_str = pop!(settings_dict, :linsys_solver)
    linsys_str = get(settings_dict, :linsys_solver, nothing)

    if linsys_str != nothing
         # Check type
        if !isa(linsys_str, String)
            error("linsys_solver is required to be a string")
        end

         # Convert to lower case
        linsys_str = lowercase(linsys_str)

        if linsys_str == "suitesparse ldl"
            settings_dict[:linsys_solver] = SUITESPARSE_LDL_SOLVER
        elseif linsys_str == "mkl pardiso"
            settings_dict[:linsys_solver] = MKL_PARDISO_SOLVER
        elseif linsys_str == ""
            settings_dict[:linsys_solver] = SUITESPARSE_LDL_SOLVER
        else
            warn("Linear system solver not recognized. Using default SuiteSparse LDL")
            settings_dict[:linsys_solver] = SUITESPARSE_LDL_SOLVER

        end
    end
end


text:
module AustralianWages

using Matte

const title = "Matte Example :: Australian Wages by Occupation"

function ui()
    sidebar_layout(
        side_panel(
            h3("Choose an occupation"),
            br(),
            selector("sub_major", "Occupation major group", "sub_groups"; multiple = true, autocomplete = true),
            show_if("show_unit",
                selector("unit", "Occupation unit group", "unit_groups", multiple = true, autocomplete = true),
                show_if("show_occ",
                    selector("selected_occs", "Occupation", "occupation_list", multiple = true, autocomplete = true)
                )
            ),
            h3("Sex"),
            radio("sex", ["Male", "Female", "Total"], ["Male", "Female", "Total"]),
            button("calculate_plot", "Create plot!"; color = "primary")
        ),
        main_panel(
            h1("Australian Wages by Occupation"),
            br(),
            hide_if("plot_loaded",
                circular_loader()
            ),
            show_if("plot_loaded",
                plots_output("plot_output")
            )
        )
    ),
    footer(p("Powered by Matte.jl; Data from Australian Taxation Office, 2016/17"))
end

module Server

import CSV, Matte
using DataFrames, Plots
using Plots.PlotMeasures

## DATA IMPORT
wage_data = CSV.read("ts17individual15occupationgender.csv")

function show_unit(sub_major)
    length(sub_major) > 0
end

function show_occ(unit)
    length(unit) > 0
end

function sub_groups()
    wage_data[:, Symbol("Sub-major group")]
end

function unit_groups(sub_major)
    if !(typeof(sub_major) <: AbstractArray)
        sub_major = [sub_major] # selector with multiple selector returns element if one selected, array otherwise
    end
    keep = [x in sub_major for x in wage_data[:, Symbol("Sub-major group")]]
    wage_data[keep, Symbol("Unit group")]
end

function occupation_list(unit)
    if !(typeof(unit) <: AbstractArray)
        unit = [unit] # selector with multiple selector returns element if one selected, array otherwise
    end
    keep = [x in unit for x in wage_data[:, Symbol("Unit group")]]
    wage_data[keep, Symbol("Occupation")]
end

function plot_loaded()
    false
end

function plot_output(selected_occs, sex, calculate_plot, session)
    if calculate_plot && sex != nothing && length(selected_occs) > 0
        Matte.update_output("plot_loaded", false, session)
        session.plot_loaded = true
        keep_occ = [x in selected_occs for x in wage_data[:, Symbol("Occupation")]]
        keep_sex = wage_data[:, :Sex] .== sex
        plot_data = wage_data[keep_occ .& keep_sex, Symbol("Average taxable income")]
        plot_labels = wage_data[keep_occ .& keep_sex, Symbol("Occupation")]
        p = plot(plot_labels, plot_data, linetype = :bar, legend = false, left_margin = 100px)
        Matte.update_output("plot_loaded", true, session)
        p
    else
        nothing
    end
end

end

end

ext: jl


text:
export LightAutomaton, AbstractAutomaton, OneStateAutomaton
export states, modes, nstates, nmodes, transitions, ntransitions
export source, event, symbol, target, transitiontype
export add_transition!, has_transition, rem_transition!, rem_state!, add_state!
export in_transitions, out_transitions

abstract type StateProperty{T} end
function typed_map(T::Type, f::Function, sp::StateProperty)
    new_sp = state_property(sp.automaton, T)
    for s in states(sp.automaton)
        new_sp[s] = f(sp[s])
    end
    return new_sp
end

abstract type TransitionProperty{T} end
function typed_map(T::Type, f::Function, tp::TransitionProperty)
    new_tp = transition_property(tp.automaton, T)
    for t in transitions(tp.automaton)
        new_tp[t] = f(tp[t])
    end
    return new_tp
end

"""
    AbstractAutomaton

Abstract type for an automaton.
"""
abstract type AbstractAutomaton end

Base.broadcastable(A::AbstractAutomaton) = Ref(A)

"""
    AbstractTransition

Abstract type for the transition of an automaton.
"""
abstract type AbstractTransition end

Base.broadcastable(t::AbstractTransition) = Ref(t)

"""
    states(A::AbstractAutomaton)

Returns an iterator over the states of the automaton `A`.
It has the alias `modes`.
"""
function states end
const modes = states

"""
    nstates(A::AbstractAutomaton)

Returns the number of states of the automaton `A`.
It has the alias `nmodes`.
"""
function nstates end
const nmodes = nstates

"""
    transitiontype(A::AbstractAutomaton)

Returns type of the transitions of the automaton `A`.
"""
function transitiontype end

"""
    transitions(A::AbstractAutomaton)

Returns an iterator over the transitions of the automaton `A`.
"""
function transitions end

"""
    ntransitions(A::AbstractAutomaton)

Returns the number of transitions of the automaton `A`.
"""
function ntransitions end

"""
    add_transition!(A::AbstractAutomaton, q, r, σ)

Adds a transition between states `q` and `r` with symbol `σ` to the automaton `A`.
"""
function add_transition! end

"""
    has_transition(A::AbstractAutomaton, t::AbstractTransition)::Bool

Returns a `Bool` indicating whether the automaton `A` has the transition `t`.

    has_transition(A::AbstractAutomaton, q, r)::Bool

Returns a `Bool` indicating whether the automaton `A` has a transition from
state `q` to state `r`.
"""
function has_transition end

"""
    rem_transition!(A::AbstractAutomaton, t::AbstractTransition)

Remove the transition `t` from the automaton `A`.
"""
function rem_transition! end

"""
    rem_state!(A::AbstractAutomaton, q)

Remove the state `q` to the automaton `A`.
"""
function rem_state! end

"""
    source(A::AbstractAutomaton, t::AbstractTransition)

Returns the source of the transition `t`.
"""
function source end

"""
    event(A::AbstractAutomaton, t::AbstractTransition)

Returns the event/symbol of the transition `t` in the automaton `A`.
It has the alias `symbol`.
"""
function event end
const symbol = event

"""
    target(A::AbstractAutomaton, t::AbstractTransition)

Returns the target of the transition `t`.
"""
function target end

"""
    in_transitions(A::AbstractAutomaton, s)

Returns an iterator over the transitions with target `s`.
"""
function in_transitions end

"""
    out_transitions(A::AbstractAutomaton, s)

Returns an iterator over the transitions with source `s`.
"""
function out_transitions end

include("onestate.jl")
include("light.jl")


text:
@testset "4.2.9 trig^m (a+b cos^n+c cos^(2 n))^p" begin
    (a, b, c, x, ) = @variables a b c x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*Trig(d+e*x)^m*(a+b*cos(d+e*x)^n+c*cos(d+e*x)^(2*n))^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*sin(d+e*x)^m*(a+b*cos(d+e*x)^n+c*cos(d+e*x)^(2*n))^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(d+e*x)^m*(a+b*cos(d+e*x)+c*cos(d+e*x)^2)^p=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [sin(x)^5/(a + b*cos(x) + c*cos(x)^2), x, 7, ((b^4 + 2*c^2*(a + c)^2 - 2*b^2*c*(2*a + c))*atanh((b + 2*c*cos(x))/sqrt(b^2 - 4*a*c)))/(c^4*sqrt(b^2 - 4*a*c)) - ((b^2 - c*(a + 2*c))*cos(x))/c^3 + (b*cos(x)^2)/(2*c^2) - cos(x)^3/(3*c) + (b*(b^2 - 2*c*(a + c))*log(a + b*cos(x) + c*cos(x)^2))/(2*c^4)]
    @test_int [sin(x)^3/(a + b*cos(x) + c*cos(x)^2), x, 7, -(((b^2 - 2*c*(a + c))*atanh((b + 2*c*cos(x))/sqrt(b^2 - 4*a*c)))/(c^2*sqrt(b^2 - 4*a*c))) + cos(x)/c - (b*log(a + b*cos(x) + c*cos(x)^2))/(2*c^2)]
    @test_int [sin(x)^1/(a + b*cos(x) + c*cos(x)^2), x, 3, (2*atanh((b + 2*c*cos(x))/sqrt(b^2 - 4*a*c)))/sqrt(b^2 - 4*a*c)]
    @test_int [csc(x)^1/(a + b*cos(x) + c*cos(x)^2), x, 9, -(((b^2 - 2*a*c - 2*c^2)*atanh((b + 2*c*cos(x))/sqrt(b^2 - 4*a*c)))/((a - b + c)*(a + b + c)*sqrt(b^2 - 4*a*c))) + log(1 - cos(x))/(2*(a + b + c)) - log(1 + cos(x))/(2*(a - b + c)) + (b*log(a + b*cos(x) + c*cos(x)^2))/(2*(a - b + c)*(a + b + c))]
    @test_int [csc(x)^3/(a + b*cos(x) + c*cos(x)^2), x, 10, ((b^4 + 2*c^2*(a + c)^2 - 2*b^2*c*(2*a + c))*atanh((b + 2*c*cos(x))/sqrt(b^2 - 4*a*c)))/(sqrt(b^2 - 4*a*c)*(a^2 - b^2 + 2*a*c + c^2)^2) + ((b - (a + c)*cos(x))*csc(x)^2)/(2*(a - b + c)*(a + b + c)) + ((a + 2*b + 3*c)*log(1 - cos(x)))/(4*(a + b + c)^2) - ((a - 2*b + 3*c)*log(1 + cos(x)))/(4*(a - b + c)^2) - (b*(b^2 - 2*c*(a + c))*log(a + b*cos(x) + c*cos(x)^2))/(2*(a^2 - b^2 + 2*a*c + c^2)^2)]

    @test_int [sin(x)^4/(a + b*cos(x) + c*cos(x)^2), x, 10, If(13<9, x/(2*c) + ((b^2 - c*(a + 2*c))*x)/c^3 - (2*(b*(b^2 - 2*c*(a + c)) - (b^4 + 2*c^2*(a + c)^2 - 2*b^2*c*(2*a + c))/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(c^3*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) - (2*(b^4 + 2*c^2*(a + c)^2 - 2*b^2*c*(2*a + c) + b^3*sqrt(b^2 - 4*a*c) - 2*b*c*(a + c)*sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(c^3*sqrt(b^2 - 4*a*c)*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c))) - (b*sin(x))/c^2 + (cos(x)*sin(x))/(2*c), x/(2*c) + ((b^2 - c*(a + 2*c))*x)/c^3 + (2*(b^2*(b^2 - 2*c*(a + c)) - b*sqrt(b^2 - 4*a*c)*(b^2 - 2*c*(a + c)) - 2*c*(a*b^2 - c*(a + c)^2))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(c^3*sqrt(b^2 - 4*a*c)*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) - (2*(b^4 + 2*c^2*(a + c)^2 - 2*b^2*c*(2*a + c) + b^3*sqrt(b^2 - 4*a*c) - 2*b*c*(a + c)*sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(c^3*sqrt(b^2 - 4*a*c)*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c))) - (b*sin(x))/c^2 + (cos(x)*sin(x))/(2*c))]
    @test_int [sin(x)^2/(a + b*cos(x) + c*cos(x)^2), x, 7, -(x/c) + (2*(b - (b^2 - 2*c*(a + c))/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(c*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) + (2*(b + (b^2 - 2*c*(a + c))/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(c*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c)))]
    @test_int [csc(x)^2/(a + b*cos(x) + c*cos(x)^2), x, 9, -((2*b*c*(1 + (b^2 - 2*c*(a + c))/(b*sqrt(b^2 - 4*a*c)))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/((a - b + c)*(a + b + c)*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c)))) - (2*b*c*(1 - (b^2 - 2*c*(a + c))/(b*sqrt(b^2 - 4*a*c)))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/((a - b + c)*(a + b + c)*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c))) - sin(x)/(2*(a + b + c)*(1 - cos(x))) + sin(x)/(2*(a - b + c)*(1 + cos(x)))]

    @test_int [sin(x)/(-2 + cos(x) + cos(x)^2), x, 4, (-(1/3))*log(1 - cos(x)) + (1/3)*log(2 + cos(x))]
    @test_int [sin(x)/(4 - 5*cos(x) + cos(x)^2), x, 4, (1/3)*log(1 - cos(x)) - (1/3)*log(4 - cos(x))]
    @test_int [sin(x)/(3 - 2*cos(x) + cos(x)^2), x, 3, atan((1 - cos(x))/sqrt(2))/sqrt(2)]
    @test_int [sin(x)/(13 - 4*cos(x) + cos(x)^2)^2, x, 4, (-(1/54))*atan((1/3)*(-2 + cos(x))) + (2 - cos(x))/(18*(13 - 4*cos(x) + cos(x)^2))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*sin(d+e*x)^m*(a+b*cos(d+e*x)^2+c*cos(d+e*x)^4)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cos(d+e*x)^m*(a+b*cos(d+e*x)^n+c*cos(d+e*x)^(2*n))^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(d+e*x)^m*(a+b*cos(d+e*x)+c*cos(d+e*x)^2)^p=#

    #= ::Subsubsection:: =#
    #=p>0=#

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [cos(x)^4/(a + b*cos(x) + c*cos(x)^2), x, 10, x/(2*c) + ((b^2 - a*c)*x)/c^3 - (2*(b^3 - 2*a*b*c - (b^4 - 4*a*b^2*c + 2*a^2*c^2)/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(c^3*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) - (2*(b^3 - 2*a*b*c + (b^4 - 4*a*b^2*c + 2*a^2*c^2)/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(c^3*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c))) - (b*sin(x))/c^2 + (cos(x)*sin(x))/(2*c)]
    @test_int [cos(x)^3/(a + b*cos(x) + c*cos(x)^2), x, 8, -((b*x)/c^2) + (2*(b^2 - a*c - b^3/sqrt(b^2 - 4*a*c) + (3*a*b*c)/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(c^2*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) + (2*(b^2 - a*c + b^3/sqrt(b^2 - 4*a*c) - (3*a*b*c)/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(c^2*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c))) + sin(x)/c]
    @test_int [cos(x)^2/(a + b*cos(x) + c*cos(x)^2), x, 7, x/c - (2*(b - (b^2 - 2*a*c)/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(c*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) - (2*(b + (b^2 - 2*a*c)/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(c*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c)))]
    @test_int [cos(x)^1/(a + b*cos(x) + c*cos(x)^2), x, 6, (2*(1 - b/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) + (2*(1 + b/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c)))]
    @test_int [cos(x)^0/(a + b*cos(x) + c*cos(x)^2), x, 5, (4*c*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(sqrt(b^2 - 4*a*c)*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) - (4*c*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(sqrt(b^2 - 4*a*c)*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c)))]
    @test_int [sec(x)^1/(a + b*cos(x) + c*cos(x)^2), x, 8, -((2*c*(1 + b/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(a*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c)))) - (2*c*(1 - b/sqrt(b^2 - 4*a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(a*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c))) + atanh(sin(x))/a]
    @test_int [sec(x)^2/(a + b*cos(x) + c*cos(x)^2), x, 10, (2*b*c*(1 + (b^2 - 2*a*c)/(b*sqrt(b^2 - 4*a*c)))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(a^2*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c))) + (2*b*c*(1 - (b^2 - 2*a*c)/(b*sqrt(b^2 - 4*a*c)))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(a^2*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c))) - (b*atanh(sin(x)))/a^2 + tan(x)/a]
    @test_int [sec(x)^3/(a + b*cos(x) + c*cos(x)^2), x, 12, -((2*c*(b^3 - 3*a*b*c + sqrt(b^2 - 4*a*c)*(b^2 - a*c))*atan((sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c - sqrt(b^2 - 4*a*c))))/(a^3*sqrt(b^2 - 4*a*c)*sqrt(b - 2*c - sqrt(b^2 - 4*a*c))*sqrt(b + 2*c - sqrt(b^2 - 4*a*c)))) + (2*c*(b^3 - 3*a*b*c - sqrt(b^2 - 4*a*c)*(b^2 - a*c))*atan((sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*tan(x/2))/sqrt(b + 2*c + sqrt(b^2 - 4*a*c))))/(a^3*sqrt(b^2 - 4*a*c)*sqrt(b - 2*c + sqrt(b^2 - 4*a*c))*sqrt(b + 2*c + sqrt(b^2 - 4*a*c))) + atanh(sin(x))/(2*a) + ((b^2 - a*c)*atanh(sin(x)))/a^3 - (b*tan(x))/a^2 + (sec(x)*tan(x))/(2*a)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*cos(d+e*x)^m*(a+b*cos(d+e*x)^2+c*cos(d+e*x)^4)^p=#

    #= ::Section:: =#
    #=Integrands*of*the*form*tan(d+e*x)^m*(a+b*cos(d+e*x)^n+c*cos(d+e*x)^(2*n))^p=#
    end

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.9 trig^m (a+b cos^n+c cos^(2 n))^p.jl


text:
using Test
using CUDA: CUDA, functional
using AutoGrad: AutoGrad, gcnode, set_gc_function, @diff
using Knet.AutoGrad_gpu: knetgcnode
using Knet.KnetArrays: KnetArray
using Knet.Ops20: RNN

if CUDA.functional(); @testset "gcnode" begin

    # 506: knetgcnode garbage collects rnn fields
    M1 = RNN(2,3; atype=KnetArray{Float32})
    M1.h = M1.c = 0
    M1.dx = M1.dhx = M1.dcx = nothing
    xcpu = randn(Float32,2,4,8)
    xgpu = KnetArray(xcpu)
    M1(xgpu) # sets M1.h,c
    @diff sum(M1(xgpu)) # sets M1.h,c,dx,dhx,dcx
    @test pointer(M1.h) != C_NULL
    @test pointer(M1.c) != C_NULL
    @test pointer(M1.dx) != C_NULL
    @test pointer(M1.dhx) != C_NULL
    @test pointer(M1.dcx) != C_NULL

end; end


text:
"""
    FiniteDifferencesBackend(method::FiniteDifferenceMethod = central_fdm(5, 1))

Differentiation backend based on the FiniteDifferences package.
"""
struct FiniteDifferencesBackend{TM<:FiniteDifferenceMethod} <: AbstractDiffBackend
    method::TM
end

function FiniteDifferencesBackend()
    return FiniteDifferencesBackend(central_fdm(5, 1))
end

function _derivative(f, t, backend::FiniteDifferencesBackend)
    return backend.method(f, t)
end

function _gradient(f, p, backend::FiniteDifferencesBackend)
    return FiniteDifferences.grad(backend.method, f, p)[1]
end

function _jacobian(f, p, backend::FiniteDifferencesBackend)
    return FiniteDifferences.jacobian(backend.method, f, p)[1]
end

if default_differential_backend() === NoneDiffBackend()
    set_default_differential_backend!(FiniteDifferencesBackend())
end


text:
"""
`module  JAC.Radial`
    ... a submodel of JAC that contains all structs and methods to deal with the radial grid, radial orbitals and potentials.
"""
module Radial

    using  QuadGK, Printf, JAC, ..Basics

    export Grid, Potential, Orbital, SingleElecSpectrum

    """
    `abstract type Radial.AbstractMesh`
        ... defines an abstract (singleton) type for the radial mesh.

        + struct MeshNone       ... if an empty instance of Radial.Grid() need to be used.
        + struct MeshGrasp      ... to use a GRASP-like grid with Newton-Cortes integration rules.
        + struct MeshGL         ... to use a Gauss-Legendre type grid where the mesh points between the break points of the t-grid
                                    are chosen to support a Gauss-Legendre integration in this potential.
    """
    abstract type  AbstractMesh                     end
    struct   MeshNone   <:  Radial.AbstractMesh     end
    struct   MeshGrasp  <:  Radial.AbstractMesh     end
    struct   MeshGL     <:  Radial.AbstractMesh     end

    """
    `struct  Radial.Grid`  ... defines a type for the radial grid which contains all information about the grid parameters, the genration
                               of the B-spline basis as well as for performing radial integrations.

        ** Physical grid parameter **
        + rnt            ::Float64           ... smalles grid point > 0.
        + h              ::Float64           ... stepsize in the construction of the exponential grid.
        + hp             ::Float64           ... asymptotic stepsize of the log-lin grid.
        + NoPoints       ::Int64             ... No. of grid points so that r[NoPoints] coincides also with the largest break point of
                                                 the B-spline knot.
        ** B-spline grid parameters and break points **
        + tL             ::Array{Float64,1}  ... radial break points for the B-splines of the large component.
        + tS             ::Array{Float64,1}  ... radial break points for the B-splines of the small component.
        + ntL            ::Int64             ... number of break points in the t-grid of the large component.
        + ntS            ::Int64             ... number of break points in the t-grid of the small component.
        ##x + nth            ::Int64             ... take each nth-point from the 'physical grid' as break point.
        + orderL         ::Int64             ... B-spline order of large components of the large component
        + orderS         ::Int64             ... B-spline order of small components
        + nsL            ::Int64             ... number of B-splines for large components
        + nsS            ::Int64             ... number of B-splines for small components
        + orderGL        ::Int64             ... order of the Gauss-Lengendre integration if mesh == Radial.MeshGL(); this is also
                                                 equal of how the break points are generated as nth-point from the physical grid points.
        ** Radial mesh points **
        + meshType       ::Radial.AbstractMesh
        ##x + nr             ::Int64             ... number of mesh points in the r-grid.
        + r              ::Array{Float64,1}  ... radial grid points
        + rp             ::Array{Float64,1}  ... derivative of the radial grid at the grid points
        + rpor           ::Array{Float64,1}  ... rp over r
        + wr             ::Array{Float64,1}  ... integration weights for all grid points, for instance, GL weights.
    """
    struct  Grid
        rnt              ::Float64
        h                ::Float64
        hp               ::Float64
        NoPoints         ::Int64
        #
        tL               ::Array{Float64,1}
        tS               ::Array{Float64,1}
        ntL              ::Int64
        ntS              ::Int64
        ##x nth              ::Int64
        orderL           ::Int64
        orderS           ::Int64
        nsL              ::Int64
        nsS              ::Int64
        orderGL          ::Int64
        #
        meshType         ::Radial.AbstractMesh
        ##x nr               ::Int64
        r                ::Array{Float64,1}
        rp               ::Array{Float64,1}
        rpor             ::Array{Float64,1}
        wr               ::Array{Float64,1}
    end

    """
    `Radial.Grid()` ... constructor to define an 'empty' grid.
    """
    function Grid()
        Radial.Grid(0., 0., 0., 0,   Float64[], Float64[], 0, 0, 0, 0, 0, 0, 0,   Radial.MeshNone(), Float64[], Float64[], Float64[], Float64[] )
    end

    """
    `Radial.Grid(exponential::Bool; printout::Bool=false)`
        ... constructor to define either a standard 'exponential' mesh (true) or a 'log-lin' mesh (false).
            The standard mesh type is Radial.MeshGL().
    """
    function Grid(exponential::Bool; printout::Bool=false)
        # The standard mesh type is Radial.MeshGL() with some proper parameters
        nth = orderL = 7;    orderS = 8;   orderGL = 7;    meshType = MeshGL()

        if  exponential
            NoPoints = 392;    NoPoints = NoPoints - rem(NoPoints,orderGL)
            grid = Radial.Grid(2.0e-6, 5.0e-2,     0., NoPoints,   Float64[], Float64[], 0, 0, orderL, orderS, 0, 0, orderGL, meshType,
                               Float64[], Float64[], Float64[], Float64[] )
        else
            NoPoints = 600;    NoPoints = NoPoints - rem(NoPoints,orderGL)
            grid = Radial.Grid(2.0e-6, 5.0e-2, 2.0e-2, NoPoints,   Float64[], Float64[], 0, 0, orderL, orderS, 0, 0, orderGL, meshType,
                               Float64[], Float64[], Float64[], Float64[] )
        end
        return( Radial.determineGrid(grid, printout=printout) )
    end

    """
    `Radial.Grid(gr::Radial.Grid;`

            rnt=..,     h=..,       hp=..,      rbox=..,    orderL=..,  orderS=..,  orderGL=..,    meshType..,  printout=..)

        ... constructor for modifying the given Radial.Grid by 'overwriting' the previously selected parameters.
    """
    function Grid(gr::Radial.Grid;
        rnt::Union{Nothing,Float64}=nothing,        h::Union{Nothing,Float64}=nothing,      hp::Union{Nothing,Float64}=nothing,
        rbox::Union{Nothing,Float64}=nothing,       nth::Union{Nothing,Int64}=nothing,      orderL::Union{Nothing,Int64}=nothing,
        orderS::Union{Nothing,Int64}=nothing,       orderGL::Union{Nothing,Int64}=nothing,  meshType::Union{Nothing,Radial.AbstractMesh}=nothing,
        printout::Bool=false)

        if  rnt      == nothing   rntx      = gr.rnt        else    rntx      = rnt       end
        if  h        == nothing   hx        = gr.h          else    hx        = h         end
        if  hp       == nothing   hpx       = gr.hp         else    hpx       = hp        end
        if  rbox     == nothing   rboxx     = nothing       else    rboxx     = rbox      end
        ##x if  nth      == nothing   nthx      = gr.nth        else    nthx      = nth       end
        if  orderL   == nothing   orderLx   = gr.orderL     else    orderLx   = orderL    end
        if  orderS   == nothing   orderSx   = gr.orderS     else    orderSx   = orderS    end
        if  orderGL  == nothing   orderGLx  = gr.orderGL    else    orderGLx  = orderGL   end
        if  meshType == nothing   meshTypex = gr.meshType   else    meshTypex = meshType  end

        if      rboxx == nothing    NoPointsx = gr.NoPoints - rem(gr.NoPoints, orderGLx)
        elseif  rboxx  > 0.         NoPointsx = Radial.determineNoPoints(rntx, hx, hpx, rboxx, orderGLx)
        else    error("stop a")
        end

        grid = Radial.Grid(rntx, hx, hpx, NoPointsx, Float64[], Float64[], 0, 0, orderLx, orderSx, 0, 0, orderGLx, meshTypex,
                           Float64[], Float64[], Float64[], Float64[] )
        return( Radial.determineGrid(grid, printout=printout) )
    end

    # `Base.show(io::IO, grid::Radial.Grid)`  ... prepares a proper printout of the variable grid::Radial.Grid.
    function Base.show(io::IO, grid::Radial.Grid)

        sa = Base.string(grid::Radial.Grid);    print(io, sa * "\n")

        nr = grid.NoPoints;    ntL = length(grid.tL);    ntS = length(grid.tS);    if  nr < 6   return( nothing )    end
        print(io, "r:    ", grid.r[1:3],    "  ...  ", grid.r[nr-2:nr],         "\n")
        if  grid.meshType == Radial.MeshGrasp()
            print(io, "rp:   ", grid.rp[1:3],   "  ...  ", grid.rp[nr-2:nr],    "\n")
            print(io, "rpor: ", grid.rpor[1:3], "  ...  ", grid.rpor[nr-2:nr],  "\n")
        end
        print(io, "wr:   ", grid.wr[1:3],   "  ...  ", grid.wr[nr-2:nr],        "\n")
        ## print(io, "tL:   ", grid.tL[1:3],   "  ...  ", grid.tL[ntL-2:ntL],   "\n")
        print(io, "tS:   ", grid.tS[1:3],   "  ...  ", grid.tS[ntS-2:ntS]           )
    end

    # `Base.string(grid::Radial.Grid)`  ... provides a String notation for the variable grid::Radial.Grid.
    function Base.string(grid::Radial.Grid)
        if   grid.NoPoints == 0
            sa = "Radial grid not defined;  NoPoints = $(grid.NoPoints) ..."
        else
            sa = "Radial grid:  rnt = $(grid.rnt),  h = $(grid.h),  hp = $(grid.hp),  NoPoints = $(grid.NoPoints),  "
            sa = sa * "ntL = $(grid.ntL),  ntS = $(grid.ntS), "
            sa = sa * "orderL = $(grid.orderL),  orderS = $(grid.orderS),  nsL = $(grid.nsL),  nsS = $(grid.nsS),  mesh = $(grid.meshType), ...  "
        end
    end

    """
    `Radial.determineGrid(grid::Radial.Grid; printout::Bool=false)`
        ... determines the detailed radial grid due to the given gridType and parameters; a gr::Radial.Grid is returned
    """
    function determineGrid(grid::Radial.Grid; printout::Bool=false)

        # Read the meshType and general parameters from the given grid
        meshType = grid.meshType
        rnt    = grid.rnt;       h      = grid.h;         hp      = grid.hp;        NoPoints = grid.NoPoints
        nth    = grid.orderGL;   orderL = grid.orderL;    orderS  = grid.orderS;    orderGL  = grid.orderGL
        nr     = grid.NoPoints;  r      = zeros(nr);      rp      = zeros(nr);      rpor     = zeros(nr);      nsL = 0;   nsS = 0
        ##x @show nr, NoPoints

        # Ensure that the largest grid points is always consistent with the largest break point of the B-splines
        if  NoPoints != NoPoints - rem(NoPoints,orderGL)    error("stop a")     end

        # Now define the physical grid due to the given parameters: either an exponential or exponential-linear grid
        if  hp == 0.
            r[1] = rnt;   rp[1] = rnt;   eph = exp(h);    ett = 1.0

            for  i  in 2:nr
                ett     = eph * ett
                ettm1   = ett - 1.0
                r[i]    = rnt * ettm1
                rp[i]   = rnt * ett
                rpor[i] = ett/ettm1
            end

        elseif  hp != 0.
            function rprime(r :: Float64)
                return( 1/ (1/(r + rnt) + h/hp) )
            end

            function f(r :: Float64, i :: Int)
                return( log( r/rnt + 1) + h/hp * r - (i - 1) * h )
            end

            function fprime(r :: Float64)
                return( 1. / (r + rnt) + h/hp )
            end

            r[1]    = 0.
            rp[1]   = rprime(0.)
            rpor[1] = 0.

            rc = 0.
            rn = 0.

            for i = 2:nr
                rn = rn + rnt
                while ((abs((rc - rn)/rn) > 100 * eps(Float64)))
                    rc = rn
                    rn = rc - f(rc, i) / fprime(rc)
                end
                r[i]    = rn
                rp[i]   = rprime(rn)
                rpor[i] = rp[i] / r[i]
            end
        else
            error("stop a")
        end

        # Define the radial break points and the number of B-splines
        ntL = orderL;    ntS = orderS;    tL = zeros( orderL );    tS = zeros( orderS )
        for  i = nth:nth:NoPoints
            ntL = ntL + 1;    push!( tL, r[i]);    ntS = ntS + 1;    push!( tS, r[i])
        end
        for   i = ntL+1:ntL+orderL-1   push!( tL, tL[ntL] )   end;    nsL = ntL - 1;   ntL = length(tL)
        for   i = ntS+1:ntS+orderS-1   push!( tS, tS[ntS] )   end;    nsS = ntS - 1;   ntS = length(tS)
        ##x @show r
        ##x @show " "
        ##x @show t
        ##x ## nsL = nt + orderL - 1 - 8;    nsS = nt + orderS - 1 - 8;    nr = nrmax - 10
        ##x ## nsL = nt + orderL -3;    nsS = nt + orderS - 3;    nr = nrmax - 2
        ##x nsL = nt - 3;    nsS = nt - 2
        ##x for   ip = nt+1:nt+orderL    push!( t, t[nt] )   end
        ##x @show nsL, nsS, nt+orderL, t[nt+orderL]
        ##x @show t

        # Define the radial grid due to the definition
        wr = zeros(nr)
        ##x @show nr, NoPoints
        if  meshType == Radial.MeshGL()
            nr = 0;   r = Float64[];   rp = Float64[];   rpor = Float64[];   wr = Float64[];    rlow = 0.
            wax = gauss(orderGL);   ra = wax[1];   wa = wax[2]
            for  i = 1:length(tL)
                rstep = tL[i] - rlow
                if  rstep > 0
                    bma = tL[i] - rlow;   bpa = tL[i] + rlow
                    for  j = 1:orderGL
                       nr = nr + 1;    push!(r, ra[j] * bma/2. + bpa/2. );    push!(wr, wa[j] * bma/2. )
                    end
                end
                rlow = tL[i]
            end
            #
            if  printout
                println("Define a radial grid of type $meshType with $nr grid points")
                println(" [rnt=" * @sprintf("%.3e",rnt) * ", h=" * @sprintf("%.3e",h) *
                        ", hp=" * @sprintf("%.3e",hp) * ", NoPoints=$NoPoints, r_max=" * @sprintf("%.3e",r[nr]) * ";")
                println("  B-splines with break points at every $(nth)th point, nsL=$nsL, nsS=$nsS, orderL=$orderL, orderS=$orderS, orderGL=$orderGL,  " *
                        "ntL=$ntL, ntS=$ntS] ")
            end
        end

        Grid(rnt, h, hp, NoPoints,  tL, tS, ntL, ntS, orderL, orderS, nsL, nsS, orderGL,   meshType, r, rp, rpor, wr)
    end

    """
    `struct  Radial.GridGL`  ... defines a type for Gauss-Legendre grid of given order but where the interval is hard-coded due to
                                 the given keystring

        + nt           :Int64                 ... number of mesh points in the grid.
        + t            ::Array{Float64,1}     ... mesh points in the variable t.
        + wt           ::Array{Float64,1}     ... weights to the mesh points in t.
    """
    struct  GridGL
        nt             ::Int64
        t              ::Array{Float64,1}
        wt             ::Array{Float64,1}
    end

    """
    `Radial.GridGL("QED", orderGL::Int64; printout::Bool=false)`
        ... constructor to define Gauss-Legendre grid for the typical QED computation in the interval [1.0, infinity].
    """
    function GridGL(sa::String, orderGL::Int64; printout::Bool=false)
        !(sa == "QED")  && error("Unrecognized keystring; sa = $sa")
        txlow = 1.;    t = Float64[];    wt = Float64[];    nt = 0
        for i = 1:100000
            # Define the exponential increase (1.5) and the maximum size (infinity=150.)
            txup = txlow * 1.3;    if  txup > 100.   break   end
            #
            wax = QuadGK.gauss(orderGL);   tx = wax[1];   wtx = wax[2]
            bma  = txup - txlow;   bpa = txup + txlow
            for  j = 1:orderGL
                nt = nt + 1;    push!(t, tx[j] * bma/2. + bpa/2. );    push!(wt, wtx[j] * bma/2. )
            end
            txlow = txup
        end
        if  printout   println("Gauss-Legendre grid with $nt mesh points from t = 1 ... $txlow.")   end

        GridGL(nt, t, wt)
    end

    """
    `Radial.GridGL("Finite", tmin::Float64, tmax::Float64, orderGL::Int64; printout::Bool=false)`
        ... constructor to define Gauss-Legendre grid in the interval [tmin, tmax].
    """
    function GridGL(sa::String, tmin::Float64, tmax::Float64, orderGL::Int64; printout::Bool=false)
        !(sa == "Finite")  && error("Unrecognized keystring; sa = $sa")
        t = Float64[];    wt = Float64[]

        wax = QuadGK.gauss(orderGL);    t = wax[1];     wt = wax[2]
        fac1 = 0.5 * ( tmax - tmin );   fac2 = 0.5 * ( tmax + tmin )

        for  j = 1:orderGL
                t[j]    = fac1 * t[j] + fac2
                wt[j]   = fac1 * wt[j]
        end

        if  printout   println("Gauss-Legendre grid with $orderGL mesh points from t = $tmin ... $tmax.")   end

        GridGL(orderGL, t, wt)
    end

    # `Base.show(io::IO, grid::GridGL)`  ... prepares a proper printout of the variable grid::GridGL.
    function Base.show(io::IO, grid::GridGL)
        print(io, "Gauss-Legendre grid with $(grid.nt) mesh points: \n")
        print(io, "t:   ", grid.t[1:5],    "  ...  ", grid.t[grid.nt-4:grid.nt],     "\n")
        print(io, "w:   ", grid.wt[1:5],   "  ...  ", grid.wt[grid.nt-4:grid.nt],    "\n")
    end

    """
    `struct  Radial.GridGH`  ... defines a type for Gauss-Hermite grid of given order

        + nt           ::Int64                ... number of mesh points in the grid.
        + t            ::Array{Float64,1}     ... mesh points in the variable t.
        + wt           ::Array{Float64,1}     ... weights to the mesh points in t.
    """
    struct  GridGH
        nt             ::Int64
        t              ::Array{Float64,1}
        wt             ::Array{Float64,1}
    end

    """
    `Radial.GridGH(orderGH::Int64; printout::Bool=false)`
        ... constructor to define Gauss-Hermite grid.
    """
    function GridGH(orderGH::Int64; printout::Bool=false)
        t = Float64[];    wt = Float64[]

        wax = gausshermite(orderGH) #this function is included in the package FastGaussQuadrature

        t = wax[1]
        wt = wax[2]

        if  printout   println("Gauss-Hermite grid with $orderGH mesh points.")   end

        GridGH(orderGH, t, wt)
    end

    # `Base.show(io::IO, grid::GridGH)`  ... prepares a proper printout of the variable grid::GridGH.
    function Base.show(io::IO, grid::GridGH)
        print(io, "Gauss-Hermite grid with $(grid.nt) mesh points: \n")
        print(io, "t:   ", grid.t[1:5],    "  ...  ", grid.t[grid.nt-4:grid.nt],    "\n")
        print(io, "w:   ", grid.wt[1:5],   "  ...  ", grid.wt[grid.nt-4:grid.nt],    "\n")
    end

    """
    `struct  Radial.Potential`  ... defines a struct for a local radial potential.

        + name           ::String            ... A name for the potential.
        + Zr             ::Array{Float64,1}  ... radial potential function Z(r) = - r * V(r) as usual in atomic structure theory.
        + grid           ::RadialGrid        ... radial grid on which the potential is defined.
    """
    struct  Potential
        name             ::String
        Zr               ::Array{Float64,1}
        grid             ::Radial.Grid
    end

    """
    `Radial.Potential()`  ... constructor to define an 'empty' instance of the radial potential.
    """
    function Potential()
        Potential("", Float64[], Radial.Grid())
    end

    # `Base.show(io::IO, potential::Radial.Potential)`  ... prepares a proper printout of the variable potential::Radial.Potential.
    function Base.show(io::IO, potential::Radial.Potential)

        sa = Base.string(potential);    print(io, sa * "\n")

        n = length(potential.Zr);                    if  n < 6   return( nothing )    end
        print(io, "Zr:    ", potential.Zr[1:23],    "  ...  ", potential.Zr[n-22:n],    "\n")
        print(io, potential.grid)
    end

    # `Base.string(potential::Radial.Potential)`  ... provides a String notation for the variable potential::Radial.Potential.
    function Base.string(potential::Radial.Potential)
        if  length(potential.Zr)  == 0
            sa = "Radial potential not yet defined; kind = $(potential.name) ..."
        else
            sa = "$(potential.name) (radial) potential ... defined on $(length(potential.Zr)) grid points ..."
        end
    end

    """
    `struct  Radial.Orbital`
        ... defines a type for a single-electron radial orbital function with a large and small component, and which can refer to
            either the standard or an explicitly given grid due to the logical flag useStandardGrid. Bound-state orbitals with energy < 0 are
            distinguished from free-electron orbitals by the flag isBound.

        + subshell        ::Subshell          ... Relativistic subshell.
        + isBound         ::Bool              ... Logical flag to distinguish between bound (true) and free-electron orbitals (false).
        + useStandardGrid ::Bool              ... Logical flag for using the standard grid (true) or an explicitly given grid (false).
        + energy          ::Float64           ... Single-electron energies of bound orbitals are always negative.
        + P               ::Array{Float64,1}  ... Large and ..
        + Q               ::Array{Float64,1}  ... small component of the radial orbital.
        + Pprime          ::Array{Float64,1}  ... dP/dr.
        + Qprime          ::Array{Float64,1}  ... dQ/dr.
        + grid            ::Array{Float64,1}  ... explic. defined radial grid array for P, Q, if StandardGrid = false.
    """
    mutable struct Orbital
        subshell          ::Subshell
        isBound           ::Bool
        useStandardGrid   ::Bool
        energy            ::Float64
        P                 ::Array{Float64,1}
        Q                 ::Array{Float64,1}
        Pprime            ::Array{Float64,1}
        Qprime            ::Array{Float64,1}
        grid              ::Radial.Grid
    end

    """
    `Radial.Orbital(subshell::Subshell, energy::Float64)`
        ... constructor for given subshell and energy, and where useStandardGrid is set to true; the grid must be defined
            explicitly and neither the large and small components nor their derivatives are yet defined in this case.
    """
    function Orbital(subshell::Subshell, energy::Float64)
        if energy < 0    isBound = true    else    isBound = false    end
        useStandardGrid = true
        P = Array{Float64,1}[];    Q = Array{Float64,1}[];    grid = Radial.Grid()
        Pprime = Array{Float64,1}[];    Qprime = Array{Float64,1}[]

        Orbital(subshell, isBound, useStandardGrid, energy, P, Q, Pprime, Qprime, grid)
    end

    """
    `Radial.Orbital(label::String, energy::Float64)`
        ... constructor for given string identifier and energy, and where useStandardGrid is set to true; the grid must be
            defined explicitly and neither the large and small components nor their derivatives are yet defined in this case.
    """
    function Orbital(label::String, energy::Float64)
        if energy < 0    isBound = true    else    isBound = false    end

        subshell = Subshell(label);    useStandardGrid = true
        P = Array{Float64,1}[];        Q = Array{Float64,1}[];    grid = Radial.Grid()
        Pprime = Array{Float64,1}[];    Qprime = Array{Float64,1}[]

        Orbital(subshell, isBound, useStandardGrid, energy, P, Q, Pprime, Qprime, grid)
    end

    # `Base.show(io::IO, orbital::Orbital)`  ... prepares a proper printout of the variable orbital::Orbital.
    function Base.show(io::IO, orbital::Orbital)
        return()
        n = length(orbital.P)

        if   orbital.useStandardGrid
            stdgrid = Defaults.getDefaults("standard grid")
            stdgrid.NoPoints == 0                     &&   return( print("Standard grid has not yet been defined.") )

            n = min(length(orbital.P), stdgrid.NoPoints)
            n > stdgrid.NoPoints                      &&   error("length of P does not match to standard grid; n=$n  NoPoints=$(stdgrid.NoPoints) ")
            length(orbital.P) != length(orbital.Q)    &&   error("P and Q have different length")
        else
            !(n == length(orbital.P) == length(orbital.Q))   &&    error("P, Q, grid have different length")
        end

        sa = Base.string(orbital::Orbital);    print(io, sa * "\n")

        if n <= 6
            print(io, "Large component P: ", orbital.P[1:end], "\n")
            print(io, "Small component Q: ", orbital.Q[1:end], "\n")
            print(io, "Pprime:            ", orbital.Pprime[1:end], "\n")
            print(io, "Qprime:            ", orbital.Qprime[1:end], "\n")
        else
            n = length(orbital.P)
            print(io, "Large component P: ", orbital.P[1:25], "  ...  ", orbital.P[n-10:n], "\n")
            print(io, "Small component Q: ", orbital.Q[1:25], "  ...  ", orbital.Q[n-10:n], "\n")
            print(io, "Pprime:            ", orbital.Pprime[1:25], "  ...  ", orbital.Pprime[n-10:n], "\n")
            print(io, "Qprime:            ", orbital.Qprime[1:25], "  ...  ", orbital.Qprime[n-10:n], "\n")
        end

        if orbital.useStandardGrid
           n <= 6   &&   print(io, "Defined on Grid:   ", stdgrid.r[1:n], "\n")
           n >  6   &&   print(io, "Defined on Grid:   ", stdgrid.r[1:3], "  ...  ", stdgrid.r[n-3:n], "\n")
        else
           n <= 6   &&   print(io, "Defined on Grid:   ", orbital.grid.r[1:n], "\n")
           n > 6    &&   print(io, "Defined on Grid:   ", orbital.grid.r[1:5], "  ...  ", orbital.grid.r[n-4:n], "\n")
        end
    end

    # `Base.string(orbital::Orbital)`  ... provides a String notation for the variable orbital::Orbital with just basic information.
    function Base.string(orbital::Orbital)
        if  orbital.isBound           sa = "Bound-state orbital "     else    sa = "Free-electron orbital "         end
        if  orbital.useStandardGrid   sb = "the standard grid: "      else    sb = "an explicitly-defined grid: "   end
        energy = orbital.energy
        sc = string(orbital.subshell) * " with energy $energy a.u. ";   n = length(orbital.P)

        return( sa * sc * "is defined with $n (grid) points on " * sb )
    end

    """
    `struct  Radial.SingleSymOrbitals`
        ... defines a type for a (relativistic and quasi-complete) single-electron spectrum for symmetry kappa with N positive
            and/or N negative states. All these states are defined with regard to the same grid.

        + kappa        ::Int64                 ... symmetry of the one-electron spectrum
        + NoStates     ::Int64                 ... Number of positive and negative states (if onlyPositive = true)
        + onlyPositive ::Bool                  ... True if only the positive part is kept.
        + pOrbitals    ::Array{Orbital,1}      ... Positive-energy orbitals states, in increasing order.
        + nOrbitals    ::Array{Orbital,1}      ... Negative-energy orbitals states, in increasing order.
        + grid         ::RadialGrid            ... radial grid on which the states are represented.
    """
    struct SingleSymOrbitals
        kappa          ::Int64
        NoStates       ::Int64
        onlyPositive   ::Bool
        pOrbitals      ::Array{Orbital,1}
        nOrbitals      ::Array{Orbital,1}
        grid           ::Radial.Grid
    end

    """
     `Radial.SingleSymOrbitals()`  ... constructor for providing an 'empty' instance of this struct.
    """
    function SingleSymOrbitals()
        SingleSymOrbitals(0, 0, true, Orbital[], Orbital[], Radial.Grid() )
    end

    # `Base.show(io::IO, symOrbitals::SingleSymOrbitals)`  ... prepares a proper printout of the variable symOrbitals::SingleSymOrbitals.
    function Base.show(io::IO, symOrbitals::SingleSymOrbitals)
        sa = "Single-symmetry orbital set for kappa=$(symOrbitals.kappa) with $(symOrbitals.NoStates) states "
        if  symOrbitals.onlyPositive   sa = sa * "from just the positive part of the spectrum."
        else                           sa = sa * "from the positive and negative part of the spectrum."
        end
        println(io, sa)
   end

    """
    `struct  Radial.SingleElecSpectrum`
        ... defines a type for a (relativistic and quasi-complete) single-electron spectrum for different symmetries kappa and either
            (individually) N_kappa positive  or  N_kappa positive and negative states. All these orbitals are defined with regard to
            the same grid.

        + name         ::String                      ... A name for this spectrum, may contain information about the original
                                                         type of basis functions.
        + symOrbitals  ::Array{SingleSymOrbitals,1}  ... Set of orbitals of the same kappa-symmetry.
    """
    struct  SingleElecSpectrum
        name           ::String
        symOrbitals    ::Array{SingleSymOrbitals,1}
    end

    """
     `Radial.SingleElecSpectrum()`  ... constructor for providing an 'empty' instance of this struct.
    """
    function SingleElecSpectrum()
        SingleElecSpectrum("", SingleSymOrbitals[])
    end

    # `Base.show(io::IO, spectrum::SingleElecSpectrum)`  ... prepares a proper printout of the variable spectrum::SingleElecSpectrum.
    function Base.show(io::IO, spectrum::SingleElecSpectrum)
        sa = "Single-electron spectrum $(spectrum.name) for different kappa-symmetries and (individual) number of orbitals " *
             "positive and/or positive & negative energy."
        println(io, sa)
   end

    """
    `Radial.OrbitalBunge1993(subshell::Subshell, Z::Int64)`  ... to calculate the radial orbital on the standard grid.
    """
    function OrbitalBunge1993(subshell::Subshell, Z::Int64, grid::Radial.Grid)
        #x grid = Defaults.getDefaults("standard grid");
        #x qn = SubshellQuantumNumbers( string(subshell) );    n = qn[1];    kappa = qn[2];   l = qn[3]
        n = subshell.n;    kappa = subshell.kappa;   l = Basics.subshell_l(subshell)
        ##x println("n, l = $n, $l")

        wa = Basics.store("orbital functions: NR, Bunge (1993)", Z)
        if      l == 0  &&  n > wa[5][1]+0     error("orbitals with s symmetry for Z = $Z only available up to $(wa[5][1]+0) only.")
        elseif  l == 1  &&  n > wa[6][1]+1     error("orbitals with p symmetry for Z = $Z only available up to $(wa[6][1]+1) only.")
        elseif  l == 2  &&  n > wa[7][1]+2     error("orbitals with d symmetry for Z = $Z only available up to $(wa[7][1]+2) only.")
        elseif  l == 3  &&  n > wa[8][1]+3     error("orbitals with f symmetry for Z = $Z only available up to $(wa[8][1]+3) only.")
        elseif  l == 4  &&  n > wa[9][1]+4     error("orbitals with g symmetry for Z = $Z only available up to $(wa[9][1]+4) only.")
        end

        P = Float64[]
        for  k in 1:grid.NoPoints-5
            Rnl = 0.; r = grid.r[k]
            for i in 4:length(wa[5+l])
                njl = wa[5+l][i][1]
                Zjl = wa[5+l][i][2]
                Cjl = wa[5+l][i][2+n-l]
                Njl = (2*Zjl)^(njl+0.5) / sqrt( factorial(2*njl) )
                Rnl = Rnl + Cjl*Njl*r^(njl-1)*exp(-Zjl*r);
            end
            ##x println("OrbitalBunge1993-aa: k, r, Rnl = $k  $r  $Rnl ")
            push!(P, Rnl*r )
        end

        # Now define small component by apllying kinetic-balance condition
        Q = zeros(size(P, 1))
        #x println("small component not yet defined properly; P = $P ")

        dP(i) = JAC.Math.derivative(P, i)
        for i = 2:size(Q, 1)
            Q[i] = -1/(2 * JAC.Defaults.INVERSE_FINE_STRUCTURE_CONSTANT) * (dP(i) / grid.h / grid.rp[i] + kappa/grid.r[i]) * P[i]
        end
        ##x println("OrbitalBunge1993-ab: ",Q)

        o     = Orbital( subshell, true, true, -0.5, P, Q, grid)
        norma = RadialIntegrals.overlap(o, o, grid)
        o.P = o.P/sqrt(norma)
        o.Q = o.Q/sqrt(norma)

        normb = RadialIntegrals.overlap(o, o, grid)
        ##x println("OrbitalBunge1993-ac:  for subshell $subshell : norm-before = $norma, norm-after = $normb")

        return( o )
    end

    """
    `Radial.OrbitalMcLean1981(subshell::Subshell,Z::Int64)`  ... to calculate the radial orbital on the standard grid.
    """
    function compute_McLean1981(subshell::Subshell,Z::Int64)
        grid = Defaults.getDefaults("standard grid");
        #x qn = Basics.SubshellQuantumNumbers( string(subshell) );    n = qn[1];   l = qn[3]
        n = subshell.n;   l = Basics.subshell_l( subshell )

        wa = Basics.store("orbital functions: NR, McLean (1981)", Z)

        !(55 <= Z <= 92)   &&                     error("Nuclear charge must be 55 <= (Z=$Z) <= 92.")
        wa == Any[]        &&                     error("No data available for neutral Z = $Z.")

        if      l == 0     &&  n > wa[5][1]+0     error("Orbitals with s symmetry for Z = $Z only available up to $(wa[5][1]+0) only.")
        elseif  l == 1     &&  n > wa[6][1]+1     error("orbitals with p symmetry for Z = $Z only available up to $(wa[6][1]+1) only.")
        elseif  l == 2     &&  n > wa[7][1]+2     error("orbitals with d symmetry for Z = $Z only available up to $(wa[7][1]+2) only.")
        elseif  l == 3     &&  n > wa[8][1]+3     error("orbitals with f symmetry for Z = $Z only available up to $(wa[8][1]+3) only.")
        elseif  l == 4     &&  n > wa[9][1]+4     error("orbitals with g symmetry for Z = $Z only available up to $(wa[9][1]+4) only.")
        end

        P = Float64[]
        for  k in 1:grid.NoPoints-5
            Rnl = 0.; r = grid.r[k]
            for i in 4:length(wa[5+l])
                njl = wa[5+l][i][1]
                Zjl = wa[5+l][i][2]
                Cjl = wa[5+l][i][2+n-l]
                Njl = (2*Zjl)^(njl+0.5) / sqrt( factorial(2*njl) )
                Rnl = Rnl + Cjl*Njl*r^(njl-1)*exp(-Zjl*r);
            end
            println("k, r, Rnl = $k  $r  $Rnl ")
            push!(P, Rnl*r )
        end

        println("small component not yet defined properly; P = $P ")

        return( Orbital( subshell, -0.5) )
    end

    """
    `Radial.OrbitalPrimitiveSlater(subshell::Subshell, N::Int64, alpha0::Float64, beta0::Float64, grid::Radial.Grid)`
        ... to calculate a list of (non-relativistic) radial Slater primitives::Array{Orbital,1} for the subshell-symmetry on
            the given grid. All orbitals have the same subshell, isBound = true, useStandardGrid = true and energy = 0.
            The small components are constructed by applying kinetic balance to the large component, and all the primitives are
            properly normalized.
    """
    function OrbitalPrimitiveSlater(subshell::Subshell, N::Int64, alpha0::Float64, beta0::Float64, grid::Radial.Grid)
        function P(i::Int64, l::Int64, eta::Float64)   return( grid.r[i]^(l+1) * exp(-eta*grid.r[i]) )   end
        dP(i) = JAC.Math.derivative(Px, i)

        orbitals = Orbital[];    ll = Basics.subshell_l(subshell)
        for  i = 1:N
            eta = alpha0 * beta0^(i-1)
            # Px  = [for i = 1:grid.NoPoints P(i, ll, eta)]
            Qx  = zeros(size(Px, 1))
            for j = 2:size(Q, 1)
                Qx[j] = -1/(2 * JAC.Defaults.INVERSE_FINE_STRUCTURE_CONSTANT) * (dP(j) / grid.h / grid.rp[j] + kappa/grid.r[j]) * Px[j]
            end

            o     = Orbital( subshell, true, true, 0., Px, Qx, grid() )
            norma = JAC.RadialIntegrals.overlap(o, o, grid)
            o.P = o.P/sqrt(norma)
            o.Q = o.Q/sqrt(norma)

            normb = JAC.RadialIntegrals.overlap(o, o, grid)
            println("OrbitalPrimitiveSlater-aa:  for i = $i : norm-before = $norma, norm-after = $normb")

            push!( orbitals, o)
        end

        return( SingleElecSpectrum() )
    end

    """
    `Radial.determineZbar(pot::Radial.Potential)`
        ... determines the effective charge that is asymptotically seen by the electron in the potential pot.
            A Zbar::Float64 is returned.
    """
    function determineZbar(pot::Radial.Potential)
        nx = 5   # Number of grid points for determining the effective charge Zbar
        mtp = size( pot.Zr, 1);    grid = pot.grid;    meanZbar = zeros(nx);    devsZbar = zeros(nx);   ny = 0
        ##x println("mtp = $mtp,  pot.Zr = $(pot.Zr[mtp-nx+1:mtp])")
        for  i = mtp-nx+1:mtp
            ny  = ny + 1;    meanZbar[ny] = pot.Zr[mtp]
        end
        mZbar   = sum(meanZbar) / nx;   for  i = 1:nx    devsZbar[i] = (meanZbar[i] - mZbar)^2    end
        stdZbar = sqrt( sum(devsZbar) / nx )
        if  true  println(">> Radial potential with effective charge Zbar=" * @sprintf("%.4e",mZbar) *
                          " (Delta-Zbar=" * @sprintf("%.4e",stdZbar) * ") at r=" * @sprintf("%.4e",pot.grid.r[mtp]) * " a.u." )    end

        return( mZbar )
    end

    """
    `Radial.determineNoPoints(rnt::Float64, h::Float64, hp::Float64, rbox::Float64, orderGL::Int64)`
        ... determines the number of points of the physical size if the grid parameters and the box-size is given; moreover,
            it is ensured that this NoPoints is consistent with the largest break point for the B-spline grid.
            An NoPoints::Int64 is returned.
    """
    function determineNoPoints(rnt::Float64, h::Float64, hp::Float64, rbox::Float64, orderGL::Int64)

        if  hp == 0.
            NoPoints = 1;  eph = exp(h);    ett = 1.0
            while true
                NoPoints = NoPoints + 1
                ett      = eph * ett
                ettm1    = ett - 1.0
                rmax     = rnt * ettm1
                if  rmax > rbox  &&  rem(NoPoints, orderGL) == 0    break   end
            end

        elseif  hp != 0.
            NoPoints = 1;  rc = 0.;    rn = 0.

            function f(r :: Float64, i :: Int)
                return( log( r/rnt + 1) + h/hp * r - (i - 1) * h )
            end

            function fprime(r :: Float64)
                return( 1. / (r + rnt) + h/hp )
            end
            while true
                NoPoints = NoPoints + 1
                rn = rn + rnt
                while ((abs((rc - rn)/rn) > 100 * eps(Float64)))
                    rc = rn
                    rn = rc - f(rc, NoPoints) / fprime(rc)
                end
                rmax    = rn
                if  rmax > rbox  &&  rem(NoPoints, orderGL) == 0    break   end
            end
        end

        return( NoPoints )
    end

end # module


text:
"""
```
Module      : FinanceLib.Derivatives.Options
Description : Implement Derivatives modules for the FinanceLib library
Copyright   : (c) 2022 Kishaloy Neogi
License     : MIT
Maintainer  : Kishaloy Neogi
Email       : nkishaloy@yahoo.com
```
The module describes the base modules of Options Derivatives.

You may see the github repository at <https://github.com/n-kishaloy/FinanceLib.jl>
"""
module Options

using Dates: Date
using FinanceLib: yearFrac

"""
```
Call
X     # Exercise price
price # Price of option
```
"""
struct Call
  X     :: Float64  # Exercise price
  price :: Float64  # Price of option
end

"""
```
Put
  X     # Exercise price
  price # Price of option
```
"""
struct Put
  X     :: Float64  # Execise price
  price :: Float64  # Price of option
end

valuation(c :: Call; ST)  = max(0, ST - c.X)
valuation(p :: Put ; ST)  = max(0, p.X - ST)

profit(x; ST) = valuation(x; ST = ST) - x.price

breakeven(c :: Call) = c.X + c.price
breakeven(p :: Put)  = p.X - p.price

"""
```
CoveredCall = Stock - Call
  S0    # Stock value @ S₀
  call  # Call
```
"""
struct CoveredCall
  S0    :: Float64  # Stock value @ S₀
  call  :: Call     # Call
end

"""
```
ProtectivePut = Stock + Put
  S0    # Stock value @ S₀
  put   # Put
```
"""
struct ProtectivePut
  S0    :: Float64  # Stock value @ S₀
  put   :: Put      # Put
end

valuation(cv :: CoveredCall; ST) = min(ST, cv.call.X)
valuation(pv :: ProtectivePut; ST) = max(ST, pv.put.X)

profit(x :: CoveredCall; ST) = valuation(x, ST = ST) - x.S0 + x.call.price
profit(x :: ProtectivePut; ST) = valuation(x, ST = ST) - x.S0 - x.put.price

breakeven(cv :: CoveredCall) = cv.S0 - cv.call.price
breakeven(pt :: ProtectivePut) = pt.S0 + pt.put.price

end


text:
#!/usr/bin/env julia

#load path to qjulia home directory
push!(LOAD_PATH, joinpath(@__DIR__, "..", "matrix"))
push!(LOAD_PATH, joinpath(@__DIR__, "..", "core"))

using MatrixBase
using MatrixMarket
using LinearAlgebra

using QJuliaSolvers
using QJuliaEnums
using MPI
using Random

#initialize MPI
MPI.Init()

Random.seed!(2019)

matrix_path = "./nasa2146.mtx"
#matrix_path = "./nasa2910.mtx"
#matrix_path = "/home/astrel/data/nasa2146/nasa2146.mtx"
#matrix_path = "/home/astrel/data/nasa2910/nasa2910.mtx"
#matrix_path = "/home/astrel/data/nasa4704/nasa4704.mtx"
#matrix_path = "/home/astrel/data/bcsstk03/bcsstk03.mtx"

csrM = MatrixBase.CSRMat{Float64}(matrix_path)

MatrixBase.print_CSRMat_info(csrM)

data_type = typeof(csrM.csrVals[1])

b = zero(Vector{data_type}(undef, csrM.Dims[2]))
x = zero(Vector{data_type}(undef, csrM.Dims[2]))
z = zero(Vector{data_type}(undef, csrM.Dims[2]))

M(b,x) = MatrixBase.csrmv(b, csrM, x)
#or just:
#M=MatrixMarket.mmread(matrix_path)

solv_param = QJuliaSolvers.QJuliaSolverParam_qj()
# Set up parameters
#solv_param.inv_type               = QJuliaEnums.QJULIA_PCG_INVERTER
solv_param.inv_type               = QJuliaEnums.QJULIA_GMRESDR_INVERTER
solv_param.inv_type_precondition  = QJuliaEnums.QJULIA_INVALID_INVERTER
solv_param.tol                    = 1e-7
#
solv_param.maxiter                = 1000
solv_param.Nsteps                 = 2
solv_param.nKrylov                = 128
solv_param.pipeline               = 2

do_unit_source = true

if do_unit_source == true
  z .=@. 1.0 / sqrt(csrM.Dims[2])
else
  z .= rand(data_type, csrM.Dims[2])
end
#
M(b,z)
#
@time QJuliaSolvers.solve(x, b, M, M, solv_param)
#
#compute error:
z .=@. z - x
println("Solution error:  ", norm(z))
#
#compute residual:
M(z, x)
z .= b - z

println("Solution norm: ", norm(x), ", absolute residual: ", norm(z), " relative ", norm(z) / norm(b))

MPI.Finalize()

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
# Copyright (c) 2017: Hendrik Ranocha.

"""
An abstract type representing a balance law in `Dim` space dimensions.
"""
abstract type AbstractBalanceLaw{Dim} end

"""
    Burgers{T,Dim}
Burgers' equation
``
  \\partial_t u + \\partial_x \\frac{u^2}{2} = 0
``
in `Dim` space dimensions using `T` as scalar type.
"""
struct Burgers{T,Dim} <: AbstractBalanceLaw{Dim} end

function Burgers(T=Float64, Dim=1)
  Burgers{T,Dim}()
end

"""
    RiemannProblem{Model<:AbstractBalanceLaw{1}, U, T<:Real}
A Riemann problem in one space dimension for `model` with left and right values
`uₗ`,`uᵣ` at `x₀`,`t₀`.
"""
struct RiemannProblem{Model<:AbstractBalanceLaw{1}, U, T<:Real}
  model::Model
  uₗ::U
  uᵣ::U
  x₀::T
  t₀::T
end

"""
    RiemannProblem(model::AbstractBalanceLaw{1}, uₗ, uᵣ, x₀::Real=0, t₀::Real=0)
Create the `RiemannProblem` for `model` with left and right values  `uₗ`,`uᵣ` at
`x₀`,`t₀`.
"""
function RiemannProblem(model::AbstractBalanceLaw{1}, uₗ, uᵣ, x₀::Real=0, t₀::Real=0)
  @assert (typeof(uₗ) == typeof(uᵣ))
  x₀, t₀ == promote(x₀, t₀)
  RiemannProblem{typeof(model),typeof(uₗ),typeof(x₀)}(model, uₗ, uᵣ, x₀, t₀)
end

"""
    BurgersRiemannSolution{T,T1}
The solution of a Riemann problem `prob` for Burgers' equation.
"""
struct BurgersRiemannSolution{T,T1} <: Function
  prob::RiemannProblem{Burgers{T,1},T,T1}
  σ⁻::T
  σ⁺::T
end

"""
    minmax_speeds(sol::BurgersRiemannSolution)
Return the minimal and maximal speeds `σ⁻, σ⁺` that appear in the solution `sol`.
"""
function minmax_speeds(sol::BurgersRiemannSolution)
  sol.σ⁻, sol.σ⁺
end

"""
    (sol::BurgersRiemannSolution)(ξ::Real)
Evaluate the solution `sol` at the value `ξ` of the self-similarity variable
`ξ = (x - x₀) / (t - t₀)`.
"""
function (sol::BurgersRiemannSolution)(ξ::Real)
  σ⁻ = sol.σ⁻; σ⁺ = sol.σ⁺
  uₗ = sol.prob.uₗ; uᵣ = sol.prob.uᵣ

  if ξ < σ⁻
    uₗ
  elseif ξ < σ⁺
    uₗ + (ξ-σ⁻)/(σ⁺-σ⁻) * (uᵣ-uₗ)
  else
    uᵣ
  end
end

"""
    (sol::BurgersRiemannSolution)(t::Real, x::Real)
Evaluate the solution `sol` at the time and space coordinates `t` and `x`.
"""
function (sol::BurgersRiemannSolution)(x::Real, t::Real)
  x₀ = sol.prob.x₀; t₀ = sol.prob.t₀

  sol((x-x₀)/(t-t₀))
end

"""
    solve{T,T1}(prob::RiemannProblem{Burgers{T,1},T,T1})
Compute the solution of the Riemann prolem `prob`.
"""
function get_solution(prob::RiemannProblem{Burgers{T,1},T,T1}) where {T,T1}
  uₗ = prob.uₗ; uᵣ = prob.uᵣ
  if uₗ > uᵣ
    σ⁻ = σ⁺ = (uₗ + uᵣ) / 2
  else
    σ⁻ = uₗ
    σ⁺ = uᵣ
  end
  BurgersRiemannSolution(prob, σ⁻, σ⁺)
end


text:
## Extrapolations are currently just constant

function sde_extrapolant!(out,Θ,integrator::DEIntegrator,idxs,deriv::Type)
  if idxs === nothing
    recursivecopy!(out,integrator.u)
  else
    out[idxs] .= @view integrator.u[idxs]
  end
end

function sde_extrapolant(Θ,integrator::DEIntegrator,idxs,deriv::Type)
  if idxs === nothing
    return integrator.u
  else
    return integrator.u[idxs]
  end
end

function sde_interpolant!(out,Θ,integrator::DEIntegrator,idxs,deriv::Type)
  sde_interpolant!(out,Θ,integrator.dt,integrator.uprev,integrator.u,idxs,deriv)
end

function sde_interpolant(Θ,integrator::DEIntegrator,idxs,deriv::Type{T}) where T
  sde_interpolant(Θ,integrator.dt,integrator.uprev,integrator.u,idxs,deriv)
end

@muladd function sde_interpolant(Θ,dt,u0,u1,idxs::Nothing,deriv::Type{Val{0}})
  @.. (1-Θ)*u0 + Θ*u1
end

@muladd function sde_interpolant(Θ,dt,u0,u1,idxs,deriv::Type{Val{0}})
  @.. (1-Θ)*u0[idxs] + Θ*u1[idxs]
end

function sde_interpolant(Θ,dt,u0,u1,idxs::Nothing,deriv::Type{Val{1}})
  @.. (u1-u0)/dt
end

function sde_interpolant(Θ,dt,u0,u1,idxs,deriv::Type{Val{1}})
  @.. (u1[idxs]-u0[idxs])/dt
end

@muladd function sde_interpolant!(out,Θ,dt,u0,u1,idxs,deriv::Type{Val{0}})
  Θm1 = (1-Θ)
  if idxs === nothing
    @.. out = Θm1*u0 + Θ*u1
  else
    @views @.. out = Θm1*u0[idxs] + Θ*u1[idxs]
  end
end

function sde_interpolant!(out,Θ,dt,u0,u1,idxs,deriv::Type{Val{1}})
  if idxs === nothing
    @.. out = (u1-u0)/dt
  else
    @views @.. out = (u1[idxs]-u0[idxs])/dt
  end
end

@inline function current_interpolant(t::Number,integrator::DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  sde_interpolant(Θ,integrator,idxs,deriv)
end

@inline function current_interpolant(t,integrator::DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [sde_interpolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_interpolant!(val,t::Number,integrator::DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  sde_interpolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_interpolant!(val,t,integrator::DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [sde_interpolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant(t::Number,integrator::DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  sde_extrapolant(Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant!(val,t::Number,integrator::DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  sde_extrapolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant(t::AbstractArray,integrator::DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [sde_extrapolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant!(val,t,integrator::DEIntegrator,idxs=nothingderiv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [sde_extrapolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

"""
sde_interpolation(tvals,ts,timeseries,ks)

Get the value at tvals where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function sde_interpolation(tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries = id
  @inbounds tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)

  # start the search thinking it's ts[1]-ts[2]
  i₋ = 1
  i₊ = 2
  vals = map(idx) do j
    @inbounds begin
      t = tvals[j]

      if continuity === :left
        # we have i₋ = i₊ = 1 if t = ts[1], i₊ = i₋ + 1 = lastindex(ts) if t > ts[end],
        # and otherwise i₋ and i₊ satisfy ts[i₋] < t ≤ ts[i₊]
        i₊ = min(lastindex(ts), OrdinaryDiffEq._searchsortedfirst(ts,t,i₊,tdir > 0))
        i₋ = i₊ > 1 ? i₊ - 1 : i₊
      else
        # we have i₋ = i₊ - 1 = 1 if t < ts[1], i₊ = i₋ = lastindex(ts) if t = ts[end],
        # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ t < ts[i₊]
        i₋ = max(1, OrdinaryDiffEq._searchsortedlast(ts,t,i₋,tdir > 0))
        i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
      end

      dt = ts[i₊] - ts[i₋]
      Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t-ts[i₋]) / dt

      sde_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
    end
  end

  if tdir<0
    permute!(vals,idx)
  end

  DiffEqArray(vals, tvals)
end

"""
sde_interpolation(tval::Number,ts,timeseries,ks)

Get the value at tval where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function sde_interpolation(tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries = id
  @inbounds tdir = sign(ts[end]-ts[1])

  if continuity === :left
    # we have i₋ = i₊ = 1 if tval = ts[1], i₊ = i₋ + 1 = lastindex(ts) if tval > ts[end],
    # and otherwise i₋ and i₊ satisfy ts[i₋] < tval ≤ ts[i₊]
    i₊ = min(lastindex(ts), OrdinaryDiffEq._searchsortedfirst(ts,tval,2,tdir > 0))
    i₋ = i₊ > 1 ? i₊ - 1 : i₊
  else
    # we have i₋ = i₊ - 1 = 1 if tval < ts[1], i₊ = i₋ = lastindex(ts) if tval = ts[end],
    # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ tval < ts[i₊]
    i₋ = max(1, OrdinaryDiffEq._searchsortedlast(ts,tval,1,tdir > 0))
    i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
  end

  @inbounds begin
    dt = ts[i₊] - ts[i₋]
    Θ = iszero(dt) ? oneunit(tval) / oneunit(dt) : (tval-ts[i₋]) / dt
    val = sde_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
  end

  val
end

function sde_interpolation!(out,tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries = id
  @inbounds tdir = sign(ts[end]-ts[1])

  if continuity === :left
    # we have i₋ = i₊ = 1 if tval = ts[1], i₊ = i₋ + 1 = lastindex(ts) if tval > ts[end],
    # and otherwise i₋ and i₊ satisfy ts[i₋] < tval ≤ ts[i₊]
    i₊ = min(lastindex(ts), OrdinaryDiffEq._searchsortedfirst(ts,tval,2,tdir > 0))
    i₋ = i₊ > 1 ? i₊ - 1 : i₊
  else
    # we have i₋ = i₊ - 1 = 1 if tval < ts[1], i₊ = i₋ = lastindex(ts) if tval = ts[end],
    # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ tval < ts[i₊]
    i₋ = max(1, OrdinaryDiffEq._searchsortedlast(ts,tval,1,tdir > 0))
    i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
  end

  @inbounds begin
    dt = ts[i₊] - ts[i₋]
    Θ = iszero(dt) ? oneunit(tval) / oneunit(dt) : (tval-ts[i₋]) / dt
    sde_interpolant!(out,Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
  end

  out
end

function sde_interpolation!(vals,tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries = id
  @inbounds tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)

  # start the search thinking it's in ts[1]-ts[2]
  i₋ = 1
  i₊ = 2
  @inbounds for j in idx
    t = tvals[j]

    if continuity === :left
      # we have i₋ = i₊ = 1 if t = ts[1], i₊ = i₋ + 1 = lastindex(ts) if t > ts[end],
      # and otherwise i₋ and i₊ satisfy ts[i₋] < t ≤ ts[i₊]
      i₊ = min(lastindex(ts), OrdinaryDiffEq._searchsortedfirst(ts,t,i₊,tdir > 0))
      i₋ = i₊ > 1 ? i₊ - 1 : i₊
    else
      # we have i₋ = i₊ - 1 = 1 if t < ts[1], i₊ = i₋ = lastindex(ts) if t = ts[end],
      # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ t < ts[i₊]
      i₋ = max(1, OrdinaryDiffEq._searchsortedlast(ts,t,i₋,tdir > 0))
      i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
    end

    dt = ts[i₊] - ts[i₋]
    Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t-ts[i₋]) / dt

    if eltype(timeseries) <: AbstractArray
      sde_interpolant!(vals[j],Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
    else
      vals[j] = sde_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
    end
  end

  vals
end


text:
# Some branch in Italy

bus1 = Dict{String, Any}()
bus1["longitude"] = 16.2487678
bus1["latitude"] = 40.358515
bus2 = Dict{String, Any}()
bus2["longitude"] = 14.1482998
bus2["latitude"] = 37.5900782

# strategy = "all_permitted"
# strategy = "OHL_on_existing_corridors"
strategy = "cables_only"

spatial_weights, voltages, resolution, impedances = OTR.define_weights_voltages(strategy)

rgb_values, nodes_lp, boundaries, plot_dictionary = OTR.convert_image_files_to_weights(bus1, bus2)

#define inputs
input_data = Dict{String, Any}()

input_data["resolution_factor"] = 2 # resolution_factor 1,2,3, ... to speed up algorithm
input_data["algorithm_factor"] = 1 # algorithm_factor 1.....1.3 to speed up Astar algorithm, goes at expense of accuracy
input_data["distance"] = 2.5  # do not change: this is the standard resolution of the environmental data
input_data["algorithm"] = "Astar"  # "Astar" or "Dijkstra"
input_data["voltages"] = voltages # transmission voltages
input_data["spatial_weights"] = spatial_weights # spatial weights
input_data["rgb_values"] = rgb_values # Spatial data as RGB values
input_data["boundaries"] = boundaries # VBoundaries of the area (to avoid using full European range)
input_data["overlapping_area_weight"] = "average" # "average" = average weight of the spatial weights; "minimum" = minimum of the overlapping weights
input_data["strategy"] = strategy # or "OHL_on_existing_corridors" or "cables_only"
input_data["losses"] = 0.01 # proxy for losses
input_data["lifetime"] = 30 # lifetime: NOT USED in FLEXPLAN
input_data["interest"] = 0.02 # Interest: NOT USED in FLEXPLAN
input_data["technology"] = "dc" # or "dc"
input_data["power_rating"] = 2000 # power rating
input_data["start_node"] = Dict{String, Any}()
input_data["start_node"]["x"] = nodes_lp["x1"]
input_data["start_node"]["y"] = nodes_lp["y1"]
input_data["end_node"] = Dict{String, Any}()
input_data["end_node"]["x"] = nodes_lp["x2"]
input_data["end_node"]["y"] = nodes_lp["y2"]
input_data["impedances"] = impedances # Provide look-up table for OHL & OGC impedances

@testset "DC cables only" begin
    strategy = "cables_only"
    input_data["technology"] = "dc" # or "dc"
    input_data["strategy"] = strategy
    @time spatial_data, spatial_data_matrices, cost_data, equipment_data, c_tot, optimal_path, ac_dc, ac_cab, dc_cab, route_impedance, route_length  = OTR.do_optimal_routing(input_data)

    @testset "cost data" begin
        @test isapprox(cost_data["ohl"]["ohl_km"][1], 1.24362898072976e7; atol = 1e0)
        @test isapprox(cost_data["ohl"]["ohl_km"][17], 1.5231374779172337e7; atol = 1e0)
        @test isapprox(cost_data["ohl"]["ohl_km"][9639], 2.1540416786293026e7; atol = 1e0)
    end
    @testset "total cost" begin
        @test isapprox(c_tot, 9.807293638249868e8; atol = 1e0)
    end
    @testset "optimal path" begin
        @test isapprox(optimal_path[52], 13019)
        @test isapprox(optimal_path[15], 10790)
    end
    @testset "impedance_data" begin
        @test isapprox(route_impedance["r_pu"], 0.001619311421391319, atol = 1e-5)
        @test isapprox(route_impedance["x_pu"], 0.0, atol = 1e-3)
        @test isapprox(route_impedance["bc_pu"], 0.0, atol = 1e-3)
    end

    @testset "route_length" begin
        @test isapprox(route_length["total_length"], 174.544725, atol = 1e-1)
        @test isapprox(route_length["ohl_length"], 0.0, atol = 1e-3)
        @test isapprox(route_length["ugc_length"], 174.54472, atol = 1e-3)
    end
end
@testset "DC all permitted" begin
    strategy = "all_permitted"
    input_data["technology"] = "dc" # or "dc"
    input_data["strategy"] = strategy
    @time spatial_data, spatial_data_matrices, cost_data, equipment_data, c_tot, optimal_path, ac_dc, ac_cab, dc_cab, route_impedance, route_length  = OTR.do_optimal_routing(input_data)

    @testset "cost data" begin
        @test isapprox(cost_data["ohl"]["ohl_km"][1], 1.24362898072976e7; atol = 1e0)
        @test isapprox(cost_data["ohl"]["ohl_km"][17], 1.5231374779172337e7; atol = 1e0)
        @test isapprox(cost_data["ohl"]["ohl_km"][9639], 2.1540416786293026e7; atol = 1e0)
    end
    @testset "total cost" begin
        @test isapprox(c_tot, 9.807293638249868e8; atol = 1e0)
    end
    @testset "optimal path" begin
        @test isapprox(optimal_path[52], 13019)
        @test isapprox(optimal_path[15], 10790)
    end

    @testset "impedance_data" begin
        @test isapprox(route_impedance["r"], 1.6581748955047106, atol = 1e-3)
        @test isapprox(route_impedance["x_pu"], 0.0, atol = 1e-3)
        @test isapprox(route_impedance["bc_pu"], 0.0, atol = 1e-3)
    end

    @testset "route_length" begin
        @test isapprox(route_length["total_length"], 174.544725, atol = 1e-1)
        @test isapprox(route_length["ohl_length"], 0.0, atol = 1e-3)
        @test isapprox(route_length["ugc_length"], 174.54472, atol = 1e-3)
    end
end

@testset "AC cables only" begin
    strategy = "cables_only"
    input_data["technology"] = "ac" # or "dc"
    input_data["strategy"] = strategy
    @time spatial_data, spatial_data_matrices, cost_data, equipment_data, c_tot, optimal_path, ac_dc, ac_cab, dc_cab, route_impedance, route_length  = OTR.do_optimal_routing(input_data)

    @testset "cost data" begin
        @test isapprox(cost_data["ohl"]["ohl_km"][1], 4.072540970760065e7; atol = 1e0)
        @test isapprox(cost_data["ohl"]["ohl_km"][17], 4.072540970760065e7; atol = 1e0)
        @test isapprox(cost_data["ohl"]["ohl_km"][9639], 5.759442674168975e7; atol = 1e0)
    end
    @testset "total cost" begin
        @test isapprox(c_tot, 2.0405224791992671e9; atol = 1e0)
    end
    @testset "optimal path" begin
        @test isapprox(optimal_path[52], 13027)
        @test isapprox(optimal_path[15], 10790)
    end

    @testset "impedance_data" begin
        @test isapprox(route_impedance["r_pu"], 0.0005509232966856801, atol = 1e-6)
        @test isapprox(route_impedance["x"], 9.09313399161207, atol = 1e-3)
        @test isapprox(route_impedance["bc_pu"], 485.0546916571823, atol = 1e-3)
    end

    @testset "route_length" begin
        @test isapprox(route_length["total_length"], 185.57416, atol = 1e-1)
        @test isapprox(route_length["ohl_length"], 0.0, atol = 1e-3)
        @test isapprox(route_length["ugc_length"], 185.57416, atol = 1e-3)
    end
end

@testset "AC all permitted" begin
    strategy = "all_permitted"
    input_data["technology"] = "ac" # or "dc"
    input_data["strategy"] = strategy
    @time spatial_data, spatial_data_matrices, cost_data, equipment_data, c_tot, optimal_path, ac_dc, ac_cab, dc_cab, route_impedance, route_length  = OTR.do_optimal_routing(input_data)

    @testset "cost data" begin
        @test isapprox(cost_data["ohl"]["ohl_km"][1], 4.072540970760065e7; atol = 1e0)
        @test isapprox(cost_data["ohl"]["ohl_km"][17], 4.072540970760065e7; atol = 1e0)
        @test isapprox(cost_data["ohl"]["ohl_km"][9639], 5.759442674168975e7; atol = 1e0)
    end
    @testset "total cost" begin
        @test isapprox(c_tot, 2.0405224791992671e9; atol = 1e0)
    end
    @testset "optimal path" begin
        @test isapprox(optimal_path[52], 13027)
        @test isapprox(optimal_path[15], 10790)
    end

    @testset "impedance_data" begin
        @test isapprox(route_impedance["r_pu"], 0.0005509232966856801, atol = 1e-6)
        @test isapprox(route_impedance["x"], 9.09313399161207, atol = 1e-3)
        @test isapprox(route_impedance["bc_pu"], 485.0546916571823, atol = 1e-3)
    end

    @testset "route_length" begin
        @test isapprox(route_length["total_length"], 185.57416, atol = 1e-1)
        @test isapprox(route_length["ohl_length"], 0.0, atol = 1e-3)
        @test isapprox(route_length["ugc_length"], 185.57416, atol = 1e-3)
    end
end

# For plotting the pdf of the result
# OTR.plot_result(plot_dictionary, input_data, spatial_data, spatial_data_matrices, optimal_path)
#


text:
# Install
using Pkg

Pkg.add("Plots")
Pkg.add("PlotlyJS")
Pkg.add("ORCA")

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
abstract type BridgeAlgorithm <: DiffEqBase.DEAlgorithm end
struct BridgeEuler <: BridgeAlgorithm end
struct BridgeHeun <: BridgeAlgorithm end
struct BridgeMdb <: BridgeAlgorithm end
struct BridgeSRK <: BridgeAlgorithm end
struct BridgeR3 <: BridgeAlgorithm end
struct BridgeBS3 <: BridgeAlgorithm end