text:
# GFF3 Reader
# ===========

"""
    GFF3Reader(input::IO; save_directives::Bool=false)

Create a reader for data in GFF3 format.

# Arguments:
* `input`: data source
* `save_directives=false`: if true, store directive lines, which can be accessed
  with the `directives` function
"""
type Reader <: Bio.IO.AbstractReader
    state::Bio.Ragel.State
    save_directives::Bool
    targets::Vector{Symbol}
    found_fasta::Bool
    directives::Vector{Record}
    directive_count::Int
    preceding_directive_count::Int

    function Reader(input::BufferedStreams.BufferedInputStream,
                    save_directives::Bool=false,
                    skip_features::Bool=false, skip_directives::Bool=true, skip_comments::Bool=true)
        targets = Symbol[]
        if !skip_features
            push!(targets, :feature)
        end
        if !skip_directives
            push!(targets, :directive)
        end
        if !skip_comments
            push!(targets, :comment)
        end
        return new(Bio.Ragel.State(body_machine.start_state, input), save_directives, targets, false, Record[], 0, 0)
    end
end

function Reader(input::IO;
                save_directives::Bool=false,
                skip_features::Bool=false, skip_directives::Bool=true, skip_comments::Bool=true)
    return Reader(BufferedStreams.BufferedInputStream(input), save_directives, skip_features, skip_directives, skip_comments)
end

function Base.eltype(::Type{Reader})
    return Record
end

function Bio.IO.stream(reader::Reader)
    return reader.state.stream
end

function Base.eof(reader::Reader)
    return reader.state.finished || eof(reader.state.stream)
end

function Base.close(reader::Reader)
    # make trailing directives accessable
    reader.directive_count = reader.preceding_directive_count
    reader.preceding_directive_count = 0
    close(Bio.IO.stream(reader))
end

function Bio.Intervals.IntervalCollection(reader::Reader)
    intervals = collect(Interval{Record}, reader)
    return IntervalCollection(intervals, true)
end

"""
Return all directives that preceded the last GFF entry parsed as an array of
strings.

Directives at the end of the file can be accessed by calling `close(reader)`
and then `directives(reader)`.
"""
function directives(reader::Reader)
    ret = String[]
    for i in endof(reader.directives)-reader.directive_count+1:endof(reader.directives)
        push!(ret, content(reader.directives[i]))
    end
    return ret
end

"""
Return true if the GFF3 stream is at its end and there is trailing FASTA data.
"""
function hasfasta(reader::Reader)
    if eof(reader)
        return reader.found_fasta
    else
        error("GFF3 file must be read until the end before any FASTA sequences can be accessed")
    end
end

"""
Return a Bio.Seq.FASTA.Reader initialized to parse trailing FASTA data.

Throws an exception if there is no trailing FASTA, which can be checked using
`hasfasta`.
"""
function getfasta(reader::Reader)
    if !hasfasta(reader)
        error("GFF3 file has no FASTA data")
    end
    return Bio.Seq.FASTA.Reader(reader.state.stream)
end

info("compiling GFF3")
const record_machine, body_machine = (function ()
    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    rep1 = Automa.RegExp.rep1
    alt = Automa.RegExp.alt
    opt = Automa.RegExp.opt

    feature = let
        seqid = re"[a-zA-Z0-9.:^*$@!+_?\-|%]*"
        seqid.actions[:enter] = [:mark]
        seqid.actions[:exit]  = [:feature_seqid]

        source = re"[ -~]*"
        source.actions[:enter] = [:mark]
        source.actions[:exit]  = [:feature_source]

        type_ = re"[ -~]*"
        type_.actions[:enter] = [:mark]
        type_.actions[:exit]  = [:feature_type_]

        start = re"[0-9]+|\."
        start.actions[:enter] = [:mark]
        start.actions[:exit]  = [:feature_start]

        end_ = re"[0-9]+|\."
        end_.actions[:enter] = [:mark]
        end_.actions[:exit]  = [:feature_end_]

        score = re"[ -~]*[0-9][ -~]*|\."
        score.actions[:enter] = [:mark]
        score.actions[:exit]  = [:feature_score]

        strand = re"[+\-?]|\."
        strand.actions[:enter] = [:feature_strand]

        phase = re"[0-2]|\."
        phase.actions[:enter] = [:feature_phase]

        attributes = let
            char = re"[^=;,\t\r\n]"
            key = rep1(char)
            key.actions[:enter] = [:mark]
            key.actions[:exit]  = [:feature_attribute_key]
            val = rep(char)
            attr = cat(key, '=', val, rep(cat(',', val)))

            cat(rep(cat(attr, ';')), opt(attr))
        end

        cat(seqid,  '\t',
            source, '\t',
            type_,  '\t',
            start,  '\t',
            end_,   '\t',
            score,  '\t',
            strand, '\t',
            phase,  '\t',
            attributes)
    end
    feature.actions[:exit] = [:feature]

    directive = re"##[^\r\n]*"
    directive.actions[:exit] = [:directive]

    comment = re"#([^#\r\n][^\r\n]*)?"
    comment.actions[:exit] = [:comment]

    record = alt(feature, directive, comment)
    record.actions[:enter] = [:anchor]
    record.actions[:exit]  = [:record]

    blank = re"[ \t]*"

    newline = let
        lf = re"\n"
        lf.actions[:enter] = [:countline]

        cat(opt('\r'), lf)
    end

    body = rep(cat(alt(record, blank), newline))
    body.actions[:exit] = [:body]

    # look-ahead of the beginning of FASTA
    body′ = cat(body, opt('>'))

    return map(Automa.compile, (record, body′))
end)()

const record_actions = Dict(
    :feature_seqid   => :(record.seqid  = (mark:p-1) - offset),
    :feature_source  => :(record.source = (mark:p-1) - offset),
    :feature_type_   => :(record.type_  = (mark:p-1) - offset),
    :feature_start   => :(record.start  = (mark:p-1) - offset),
    :feature_end_    => :(record.end_   = (mark:p-1) - offset),
    :feature_score   => :(record.score  = (mark:p-1) - offset),
    :feature_strand  => :(record.strand = p - offset),
    :feature_phase   => :(record.phase  = p - offset),
    :feature_attribute_key => :(push!(record.attribute_keys, (mark:p-1) - offset)),
    :feature         => :(record.kind = :feature),
    :directive       => :(record.kind = :directive),
    :comment         => :(record.kind = :comment),
    :record          => quote
        Bio.ReaderHelper.resize_and_copy!(record.data, data, 1:p-1)
        record.filled = (offset+1:p-1) - offset
    end,
    :anchor          => :(),
    :mark            => :(mark = p))
eval(
    Bio.ReaderHelper.generate_index_function(
        Record,
        record_machine,
        quote
            mark = offset = 0
        end,
        record_actions))
eval(
    Bio.ReaderHelper.generate_read_function(
        Reader,
        body_machine,
        quote
            mark = offset = 0
        end,
        merge(record_actions, Dict(
            :record => quote
                Bio.ReaderHelper.resize_and_copy!(record.data, data, Bio.ReaderHelper.upanchor!(stream):p-1)
                record.filled = (offset+1:p-1) - offset
                if isfeature(record)
                    reader.directive_count = reader.preceding_directive_count
                    reader.preceding_directive_count = 0
                elseif isdirective(record)
                    reader.preceding_directive_count += 1
                    if reader.save_directives
                        push!(reader.directives, copy(record))
                    end
                    if is_fasta_directive(record)
                        reader.found_fasta = true
                        reader.state.finished = true
                    end
                end
                if record.kind ∈ reader.targets
                    found_record = true
                    @escape
                end
            end,
            :body => quote
                if data[p] == UInt8('>')
                    reader.found_fasta = true
                    reader.state.finished = true
                    # HACK: any better way?
                    cs = 0
                    @goto exit
                end
            end,
            :countline => :(linenum += 1),
            :anchor    => :(Bio.ReaderHelper.anchor!(stream, p); offset = p - 1)))))


text:
struct Knnda
    X::Array{Float64}
    y::AbstractMatrix
    fm
    nlvdis::Int
    metric::String
    h::Real
    k::Int
    tol::Real
    lev::AbstractVector
    ni::AbstractVector
end

"""
    knnda(X, y; nlvdis = 0, metric = "eucl", h = Inf, k = 1, tol = 1e-4)
k-Nearest-Neighbours weighted discrimination (kNN-DA).
* `X` : X-data.
* `y` : y-data (class membership).
* `nlvdis` : Number of latent variables (LVs) to consider in the
    global PLS used for the dimension reduction before
    calculating the dissimilarities. If `nlvdis = 0`, there is no dimension reduction.
* `metric` : Type of dissimilarity used to select the neighbors.
    Possible values are "eucl" (default; Euclidean distance)
    and "mahal" (Mahalanobis distance).
* `h` : A scalar defining the shape of the weight function. Lower is h,
    sharper is the function. See function `wdist`.
* `k` : The number of nearest neighbors to select for each observation to predict.
* `tol` : For stabilization when very close neighbors.

For each new observation to predict, i) a number of
`k` nearest neighbors is selected (= "weighting 1") and ii) a weigthed vote (= "weighting 2")
is then computed to select the most frequent class in this neighborhood.

Weightings 1 and 2 are computed from the dissimilarities between the observation
to predict and the training observations. Depending on argument `nlvdis`,
the computation is done from the raw X-data or after a dimension reduction.
In the last case, global PLS2 scores (LVs) are
computed from {`X`, Y-dummy} (where Y-dummy is the dummy table build from `y`),
and the dissimilarities are computed over these scores.

In general, for high dimensional X-data, using the Mahalanobis distance requires
preliminary dimensionality reduction of the data.
"""
function knnda(X, y; nlvdis = 0, metric = "eucl", h = Inf, k = 1, tol = 1e-4)
    X = ensure_mat(X)
    y = ensure_mat(y)
    ztab = tab(y)
    if nlvdis == 0
        fm = nothing
    else
        fm = plskern(X, dummy(y).Y; nlv = nlvdis)
    end
    return Knnda(X, y, fm, nlvdis, metric, h, k, tol,
        ztab.keys, ztab.vals)
end

"""
    predict(object::Knnda1, X)
Compute the y-predictions from the fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::Knnda, X)
    X = ensure_mat(X)
    m = size(X, 1)
    # Getknn
    if isnothing(object.fm)
        res = getknn(object.X, X; k = object.k, metric = object.metric)
    else
        res = getknn(object.fm.T, transform(object.fm, X); k = object.k, metric = object.metric)
    end
    listw = copy(res.d)
    for i = 1:m
        w = wdist(res.d[i]; h = object.h)
        w[w .< object.tol] .= object.tol
        listw[i] = w
    end
    # End
    pred = similar(object.y, m, 1)
    @inbounds for i = 1:m
        s = res.ind[i]
        pred[i, :] .= findmax_cla(object.y[s], listw[i])
    end
    (pred = pred, listnn = res.ind, listd = res.d, listw = listw)
end

ext: jl


text:
using Pkg

dependencies = [
    "IJulia",
    "Genie"
]

Pkg.add(dependencies)


text:
# arnoldi.jl
"""
    mutable struct ArnoldiFactorization{T,S} <: KrylovFactorization{T,S}

Structure to store an Arnoldi factorization of a linear map `A` of the form

```julia
A * V = V * B + r * b'
```

For a given Arnoldi factorization `fact` of length `k = length(fact)`, the basis `V` is
obtained via [`basis(fact)`](@ref basis) and is an instance of [`OrthonormalBasis{T}`](@ref
Basis), with also `length(V) == k` and where `T` denotes the type of vector like objects
used in the problem. The Rayleigh quotient `B` is obtained as
[`rayleighquotient(fact)`](@ref) and is of type [`B::PackedHessenberg{S<:Number}`](@ref
PackedHessenberg) with `size(B) == (k,k)`. The residual `r` is obtained as
[`residual(fact)`](@ref) and is of type `T`. One can also query [`normres(fact)`](@ref) to
obtain `norm(r)`, the norm of the residual. The vector `b` has no dedicated name but can be
obtained via [`rayleighextension(fact)`](@ref). It takes the default value ``e_k``, i.e. the
unit vector of all zeros and a one in the last entry, which is represented using
[`SimpleBasisVector`](@ref).

An Arnoldi factorization `fact` can be destructured as `V, B, r, nr, b = fact` with
`nr = norm(r)`.

`ArnoldiFactorization` is mutable because it can [`expand!`](@ref) or [`shrink!`](@ref).
See also [`ArnoldiIterator`](@ref) for an iterator that constructs a progressively expanding
Arnoldi factorizations of a given linear map and a starting vector. See
[`LanczosFactorization`](@ref) and [`LanczosIterator`](@ref) for a Krylov factorization that
is optimized for real symmetric or complex hermitian linear maps.
"""
mutable struct ArnoldiFactorization{T,S} <: KrylovFactorization{T,S}
    k::Int # current Krylov dimension
    V::OrthonormalBasis{T} # basis of length k
    H::Vector{S} # stores the Hessenberg matrix in packed form
    r::T # residual
end

Base.length(F::ArnoldiFactorization) = F.k
Base.sizehint!(F::ArnoldiFactorization, n) = begin
    sizehint!(F.V, n)
    sizehint!(F.H, (n * n + 3 * n) >> 1)
    return F
end
Base.eltype(F::ArnoldiFactorization) = eltype(typeof(F))
Base.eltype(::Type{<:ArnoldiFactorization{<:Any,S}}) where {S} = S

basis(F::ArnoldiFactorization) = F.V
rayleighquotient(F::ArnoldiFactorization) = PackedHessenberg(F.H, F.k)
residual(F::ArnoldiFactorization) = F.r
@inbounds normres(F::ArnoldiFactorization) = abs(F.H[end])
rayleighextension(F::ArnoldiFactorization) = SimpleBasisVector(F.k, F.k)

# Arnoldi iteration for constructing the orthonormal basis of a Krylov subspace.
"""
    struct ArnoldiIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
    ArnoldiIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth])

Iterator that takes a general linear map `f::F` and an initial vector `v₀::T` and generates
an expanding `ArnoldiFactorization` thereof. In particular, `ArnoldiIterator` iterates over
progressively expanding Arnoldi factorizations using the
[Arnoldi iteration](https://en.wikipedia.org/wiki/Arnoldi_iteration).

The argument `f` can be a matrix, or a function accepting a single argument `v`, so that
`f(v)` implements the action of the linear map on the vector `v`.

The optional argument `orth` specifies which [`Orthogonalizer`](@ref) to be used. The
default value in [`KrylovDefaults`](@ref) is to use [`ModifiedGramSchmidtIR`](@ref), which
possibly uses reorthogonalization steps.

When iterating over an instance of `ArnoldiIterator`, the values being generated are
instances of [`ArnoldiFactorization`](@ref), which can be immediately destructured into a
[`basis`](@ref), [`rayleighquotient`](@ref), [`residual`](@ref), [`normres`](@ref) and
[`rayleighextension`](@ref), for example as

```julia
for (V, B, r, nr, b) in ArnoldiIterator(f, v₀)
    # do something
    nr < tol && break # a typical stopping criterion
end
```

Since the iterator does not know the dimension of the underlying vector space of
objects of type `T`, it keeps expanding the Krylov subspace until the residual norm `nr`
falls below machine precision `eps(typeof(nr))`.

The internal state of `ArnoldiIterator` is the same as the return value, i.e. the
corresponding `ArnoldiFactorization`. However, as Julia's Base iteration interface (using
`Base.iterate`) requires that the state is not mutated, a `deepcopy` is produced upon every
next iteration step.

Instead, you can also mutate the `ArnoldiFactorization` in place, using the following
interface, e.g. for the same example above

```julia
iterator = ArnoldiIterator(f, v₀)
factorization = initialize(iterator)
while normres(factorization) > tol
    expand!(iterator, factorization)
    V, B, r, nr, b = factorization
    # do something
end
```

Here, [`initialize(::KrylovIterator)`](@ref) produces the first Krylov factorization of
length 1, and `expand!(::KrylovIterator, ::KrylovFactorization)`(@ref) expands the
factorization in place. See also [`initialize!(::KrylovIterator,
::KrylovFactorization)`](@ref) to initialize in an already existing factorization (most
information will be discarded) and [`shrink!(::KrylovFactorization, k)`](@ref) to shrink an
existing factorization down to length `k`.
"""
struct ArnoldiIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
    operator::F
    x₀::T
    orth::O
end
ArnoldiIterator(A, x₀) = ArnoldiIterator(A, x₀, KrylovDefaults.orth)

Base.IteratorSize(::Type{<:ArnoldiIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:ArnoldiIterator}) = Base.EltypeUnknown()

function Base.iterate(iter::ArnoldiIterator)
    state = initialize(iter)
    return state, state
end
function Base.iterate(iter::ArnoldiIterator, state)
    nr = normres(state)
    if nr < eps(typeof(nr))
        return nothing
    else
        state = expand!(iter, deepcopy(state))
        return state, state
    end
end

function initialize(iter::ArnoldiIterator; verbosity::Int = 0)
    # initialize without using eltype
    x₀ = iter.x₀
    β₀ = norm(x₀)
    iszero(β₀) && throw(ArgumentError("initial vector should not have norm zero"))
    Ax₀ = apply(iter.operator, x₀)
    α = dot(x₀, Ax₀) / (β₀ * β₀)
    T = typeof(α)
    # this line determines the vector type that we will henceforth use
    v = mul!(zero(T)*Ax₀, x₀, 1 / β₀) # (one(T) / β₀) * x₀ # mul!(similar(x₀, T), x₀, 1/β₀)
    if typeof(Ax₀) != typeof(v)
        r = mul!(similar(v), Ax₀, 1 / β₀)
    else
        r = rmul!(Ax₀, 1 / β₀)
    end
    βold = norm(r)
    r = axpy!(-α, v, r)
    β = norm(r)
    # possibly reorthogonalize
    if iter.orth isa Union{ClassicalGramSchmidt2,ModifiedGramSchmidt2}
        dα = dot(v, r)
        α += dα
        r = axpy!(-dα, v, r)
        β = norm(r)
    elseif iter.orth isa Union{ClassicalGramSchmidtIR,ModifiedGramSchmidtIR}
        while eps(one(β)) < β < iter.orth.η * βold
            βold = β
            dα = dot(v, r)
            α += dα
            r = axpy!(-dα, v, r)
            β = norm(r)
        end
    end
    V = OrthonormalBasis([v])
    H = T[α, β]
    if verbosity > 0
        @info "Arnoldi iteration step 1: normres = $β"
    end
    return state = ArnoldiFactorization(1, V, H, r)
end
function initialize!(iter::ArnoldiIterator, state::ArnoldiFactorization; verbosity::Int = 0)
    x₀ = iter.x₀
    V = state.V
    while length(V) > 1
        pop!(V)
    end
    H = empty!(state.H)

    v = mul!(V[1], x₀, 1 / norm(x₀))
    w = apply(iter.operator, v)
    r, α = orthogonalize!(w, v, iter.orth)
    β = norm(r)
    state.k = 1
    push!(H, α, β)
    state.r = r
    if verbosity > 0
        @info "Arnoldi iteration step 1: normres = $β"
    end
    return state
end
function expand!(iter::ArnoldiIterator, state::ArnoldiFactorization; verbosity::Int = 0)
    state.k += 1
    k = state.k
    V = state.V
    H = state.H
    r = state.r
    β = normres(state)
    push!(V, rmul!(r, 1 / β))
    m = length(H)
    resize!(H, m + k + 1)
    r, β = arnoldirecurrence!(iter.operator, V, view(H, (m+1):(m+k)), iter.orth)
    H[m+k+1] = β
    state.r = r
    if verbosity > 0
        @info "Arnoldi iteration step $k: normres = $β"
    end
    return state
end
function shrink!(state::ArnoldiFactorization, k)
    length(state) <= k && return state
    V = state.V
    H = state.H
    while length(V) > k + 1
        pop!(V)
    end
    r = pop!(V)
    resize!(H, (k * k + 3 * k) >> 1)
    state.k = k
    state.r = rmul!(r, normres(state))
    return state
end

# Arnoldi recurrence: simply use provided orthonormalization routines
function arnoldirecurrence!(
    operator,
    V::OrthonormalBasis,
    h::AbstractVector,
    orth::Orthogonalizer
)
    w = apply(operator, last(V))
    r, h = orthogonalize!(w, V, h, orth)
    return r, norm(r)
end


text:
using StatsModels: collect_matrix_terms, MatrixTerm

poly(x, n) = x^n

abstract type PolyModel end
struct PolyTerm <: AbstractTerm
    term::Symbol
    deg::Int
end
PolyTerm(t::Term, deg::ConstantTerm) = PolyTerm(t.sym, deg.n)

StatsModels.apply_schema(t::FunctionTerm{typeof(poly)}, sch, ::Type{<:PolyModel}) =
    PolyTerm(t.args_parsed...)

StatsModels.modelcols(p::PolyTerm, d::NamedTuple) =
    reduce(hcat, [d[p.term].^n for n in 1:p.deg])

struct NonMatrixTerm{T} <: AbstractTerm
    term::T
end

StatsModels.is_matrix_term(::Type{<:NonMatrixTerm}) = false
StatsModels.apply_schema(t::NonMatrixTerm, sch, Mod::Type) =
    NonMatrixTerm(apply_schema(t.term, sch, Mod))
StatsModels.modelcols(t::NonMatrixTerm, d) = modelcols(t.term, d)

struct DummyTerm <: AbstractTerm
end

@testset "Extended formula/models" begin
    d = (z = rand(10), y = rand(10), x = collect(1:10))
    sch = schema(d)

    @testset "Poly term" begin
        f = @formula(y ~ poly(x, 3))

        f_plain = apply_schema(f, sch)
        @test f_plain.rhs.terms[1] isa FunctionTerm
        @test f_plain == apply_schema(f, sch, Nothing)
        @test last(modelcols(f_plain, d)) == hcat(d[:x].^3)

        f_special = apply_schema(f, sch, PolyModel)
        @test f_special.rhs.terms[1] isa PolyTerm
        @test last(modelcols(f_special, d)) == hcat(d[:x], d[:x].^2, d[:x].^3)
    end

    @testset "Non-matrix term" begin
        f = @formula(z ~ x + y)
        f2 = term(:z) ~ term(:x) + NonMatrixTerm(term(:y))
        f3 = term(:z) ~ NonMatrixTerm(term(:x)) + term(:y)
        f4 = term(:z) ~ NonMatrixTerm.(f.rhs)
        f5 = term(:z) ~ term(:x) + NonMatrixTerm(term(:y)) + term(:y)

        @test collect_matrix_terms(f.rhs) == MatrixTerm((term(:x) + term(:y)))
        @test collect_matrix_terms(f2.rhs) ==
            (MatrixTerm((term(:x), )), NonMatrixTerm(term(:y)))
        @test collect_matrix_terms(f3.rhs) ==
            (MatrixTerm((term(:y), )), NonMatrixTerm(term(:x)))
        @test collect_matrix_terms(f4.rhs) == f4.rhs
        @test collect_matrix_terms(f5.rhs) ==
            (MatrixTerm((term(:x), term(:y))), NonMatrixTerm(term(:y)))

        f = apply_schema(f, sch)
        @test f.rhs isa MatrixTerm
        @test f.rhs == collect_matrix_terms(f.rhs)
        @test modelcols(f.rhs, d) == hcat(d.x, d.y)

        f2 = apply_schema(f2, sch)
        @test f2.rhs isa Tuple{MatrixTerm, NonMatrixTerm}
        @test f2.rhs == apply_schema((MatrixTerm(term(:x)), NonMatrixTerm(term(:y))), sch)
        @test modelcols(f2.rhs, d) == (hcat(d.x), d.y)

        # matrix term goes first
        f3 = apply_schema(f3, sch)
        @test f3.rhs isa Tuple{MatrixTerm, NonMatrixTerm}
        @test f3.rhs == apply_schema((MatrixTerm(term(:y)), NonMatrixTerm(term(:x))), sch)
        @test modelcols(f3.rhs, d) == (hcat(d.y), d.x)

        f4 = apply_schema(f4, sch)
        @test f4.rhs isa Tuple{NonMatrixTerm, NonMatrixTerm}
        @test f4.rhs == apply_schema((NonMatrixTerm(term(:x)), NonMatrixTerm(term(:y))), sch)
        @test modelcols(f4.rhs, d) == (d.x, d.y)

        # matrix terms are gathered
        f5 = apply_schema(f5, sch)
        @test f5.rhs isa Tuple{MatrixTerm, NonMatrixTerm}
        @test f5.rhs ==
            apply_schema((MatrixTerm((term.((:x, :y)))), NonMatrixTerm(term(:y))), sch)
        @test modelcols(f5.rhs, d) == (hcat(d.x, d.y), d.y)

    end

    @testset "Fallback" begin
        @test_throws ArgumentError modelcols(DummyTerm(), (a=[1], ))
    end

    @testset "Ambiguity detection" begin
        # ambiguities are introduced by adding additional methods here
        @test_broken isempty(Test.detect_ambiguities(StatsModels))
    end
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

    """
    `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.
        ** B-spline grid parameters and break points **
        + t              ::Array{Float64,1}  ... radial break points for the B-splines
        + nt             ::Int64             ... number of break points in the t-grid.
        + nth            ::Int64             ... take each nth-point from the 'physical grid' as break point.
        + orderL         ::Int64             ... B-spline order of large components
        + 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 == MeshGL
        ** Radial mesh points **
        + mesh           ::RadialMesh
        + 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
        #
        t                ::Array{Float64,1}
        nt               ::Int64
        nth              ::Int64
        orderL           ::Int64
        orderS           ::Int64
        nsL              ::Int64
        nsS              ::Int64
        orderGL          ::Int64
        #
        mesh             ::RadialMesh
        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' instance of the radial grid.

    `Radial.Grid("grid: exponential")`  ... constructor to define an instance of the (Grasp) standard exponential grid.

    `Radial.Grid("grid: by given parameters"; rnt = rntx, h = hx, hp = hpx, NoPoints = NoPointx)`
        ... constructor to define an instance of the radial grid with the given parameters.
    """
    function Grid()
        Grid(0., 0., 0., 0,   Float64[], 0, 0, 0, 0, 0, 0, 0,   MeshGrasp, 0, Float64[], Float64[], Float64[], Float64[] )
    end

    function Grid(sa::String; rnt::Float64 = 2.0e-6, h::Float64 = 5.0e-2, hp::Float64 = 0., NoPoints::Int64 = 390, printout::Bool=true)

        # Define the general parameters for the B-spline grid and radial mesh
        mesh = MeshGL  ## MeshGL   MeshGrasp
        nth = 7;    orderL = 7;    orderS = 8;   orderGL = 7;    nsL = 0;   nsS = 0
        nr = NoPoints + 10;    r = zeros(nr);    rp = zeros(nr);    rpor = zeros(nr)

        # Now define the physical grid due to the given parameters: either an exponential or exponential-linear grid
        if  hp == 0.    &&    sa in [ "grid: exponential", "grid: by given parameters"]
            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.    &&    sa in [ "grid: by given parameters"]
            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 += 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("Unsupported keystring = $sa for given parameters.")
        end

        # Define the radial break points and the number of B-splines
        nt = 0;    nrmax = 0;    t = zeros(orderL)
        for  i = 1:nth:nr
            nrmax = i;    nt = nt + 1;    push!( t, r[i])
        end
        nsL = nt + orderL - 1 - 8;    nsS = nt + orderS - 1 - 8;    nr = nrmax - 10
        for   ip = nt+1:nt+max(orderL, orderS)-1    push!( t, r[nt])   end
        ##x ns = nv + k - 1       # No. of splines
        ##x for  ip = ii+1:ii+k-1    push!( tlist, grid.r[i])   end

        # Define the radial grid due to the definition
        wr = zeros(nr)
        if  mesh == MeshGL
            nr = 1;   r = Float64[ 0.];   rp = Float64[];   rpor = Float64[];   wr = Float64[0.];    rlow = 0.
            wax = gauss(orderGL);   ra = wax[1];   wa = wax[2]
            for  i = 1:length(t)
                rstep = t[i] - rlow
                if  rstep > 0
                    bma = t[i] - rlow;   bpa = t[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 = t[i]
            end
            #
            if  printout
                println("Define a radial grid of type $mesh 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 wit break points at every $(nth)th point, nsL=$nsL, nsS=$nsS, orderL=$orderL, orderS=$orderS, orderGL=$orderGL] ")
            end
            ## for j = 1:nr   println("j = $j,  r = $(r[j]),  wr = $(wr[j])")   end
        end

        Grid(rnt, h, hp, NoPoints,  t, nt, nth, orderL, orderS, nsL, nsS, orderGL,   mesh, nr, r, rp, rpor, wr)
    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.nr;    nt = grid.nt;    if  nr < 6   return( nothing )    end
        print(io, "r:    ", grid.r[1:3],    "  ...  ", grid.r[nr-2:nr],     "\n")
        if  grid.mesh == 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, "t:    ", grid.t[1:3],    "  ...  ", grid.t[nt-2:nt]          )
    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 * "orderL = $(grid.orderL),  orderS = $(grid.orderS),  nsL = $(grid.nsL),  nsS = $(grid.nsS),  mesh = $(grid.mesh), ...  "
        end
    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.5;    if  txup > 150.   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

    # `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.Potential`  ... defines a struct for the radial potential which contains all information about its physical content.

        + 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 isStandardGrid. Bound-state orbitals with energy < 0 are
            distinguished from free-electron orbitals by the flag isBound. -- Note that the arrays P, Q and grid cannot be defined only by the
            standard constructor but are typically set explicitly after an instance of this type has been created.

        + 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 isStandardGrid is set to true; the grid must be defined
            explicitly and the large and small components are not yet defined in this case.
    """
    function Orbital(subshell::Subshell, energy::Float64)
        if energy < 0    isBound = true    else    isBound = false    end
        isStandardGrid = true
        P = Array{Float64,1}[];    Q = Array{Float64,1}[];    grid = RadialGrid()

        Orbital(subshell, isBound, isStandardGrid, energy, P, Q, grid)
    end

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

        subshell = Subshell(label);    isStandardGrid = true
        P = Array{Float64,1}[];        Q = Array{Float64,1}[];    grid = RadialGrid()

        Orbital(subshell, isBound, isStandardGrid, energy, P, Q, grid)
    end

    # `Base.show(io::IO, orbital::Orbital)`  ... prepares a proper printout of the variable orbital::Orbital.
    function Base.show(io::IO, orbital::Orbital)
        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 > stdgrid.NoPoints      &&   error("length of P does not match to standard grid")
            n != length(orbital.Q)    &&   error("P and Q have different length")
        else
            !(n == length(orbital.Q) == 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")
        else
            n = length(orbital.P)
            print(io, "Large component P: ", orbital.P[1:25], "  ...  ", orbital.P[n-24:n], "\n")
            print(io, "Small component Q: ", orbital.Q[1:25], "  ...  ", orbital.P[n-24: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()
        SingleSymOrbitals("", 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) * ")." )    end

        return( mZbar )
    end

end # module


text:
# plot recipe for plotting a Measurement object
@recipe function f(
    measurement::Measurement;
    uncertainties = true
)
    total_unc = if uncertainties == false
        0.0
    elseif uncertainties == true
        total_uncertainty(measurement, keys(measurement.uncertainties))
    else
        total_uncertainty(measurement, uncertainties)
    end

    @series begin
        seriestype --> :hline
        ribbon --> total_unc
        [measurement.value]
    end

end

# plot recipe for plotting a MeasurementDistribution object
@recipe function f(
    measurement::MeasurementDistribution;
    uncertainties = true,
    bin_names = String.(measurement.bin_names)
)
    n = length(measurement.value)
    xs = 1:n+1
    ys = [measurement.value[1], measurement.value...]

    total_uncs = if uncertainties == false
        zeros(n)
    elseif uncertainties == true
        total_uncertainty(measurement, keys(measurement.uncertainties))
    else
        total_uncertainty(measurement, uncertainties)
    end

    seriestype = get(plotattributes, :seriestype, :steppre)

    #TODO adapt for arbitrary bin width
    if seriestype == :steppre
        @series begin
            seriestype := :steppre
            ribbon --> [1, total_uncs...]
            label --> "measurements"
            xticks -->  (collect(xs).+0.5, bin_names)
            xs, ys
        end

    elseif seriestype == :scatter
        @series begin
            seriestype := :scatter
            xerror --> 0.5
            yerror --> total_uncs
            seriescolor --> :orange
            label --> "measurements"
            xticks --> (collect(xs).+0.5, bin_names)
            xs[2:end].-0.5, ys[2:end]
        end
    end
end

# plot recipe for plotting a particular Measurement from a NamedTuple
@recipe function f(
    measurements::NamedTuple{<:Any, <:Tuple{Vararg{AbstractMeasurement}}},
    name::Symbol;
    uncertainties = true
)
    @series begin
        label --> string(name)
        uncertainties --> uncertainties
        measurements[name]
    end
end

# plot recipe for plotting a particular MeasurementDistribution from a NamedTuple
@recipe function f(
    measurements::NamedTuple{<:Any, <:Tuple{Vararg{AbstractMeasurement}}},
    name::Symbol;
    uncertainties = true
)
    @series begin
        label --> string(name)
        uncertainties --> uncertainties
        measurements[name]
    end
end

function total_uncertainty(
    measurement::Union{Measurement, MeasurementDistribution},
    uncertainty::Symbol
)
    return measurement.uncertainties[uncertainty]
end

# calculate the total uncertainty of a measurement by adding the quadratures
function total_uncertainty(
    measurement::Union{Measurement, MeasurementDistribution},
    uncertainty::Union{Tuple{Vararg{Symbol}}, Vector{Symbol}}
)
    total = zeros(length(measurement.uncertainties[1]))
    for u in uncertainty
        total .+= measurement.uncertainties[u].^2
    end
    return sqrt.(total)
end


text: ccall((:entry_point, "/home/qq/MyGitHub/BacASable/Julia/JuliaCpp/libentrypoint.so"), Cvoid, ())

max_stars_repo_name: Lecrapouille/BacASable

max_issues_repo_name: Lecrapouille/BacASable

max_forks_repo_name: Lecrapouille/BacASable


text:
using YAXArrayBase, NetCDF, Zarr, Test

@testset "Reading NetCDF" begin
import Downloads
p = Downloads.download("https://www.unidata.ucar.edu/software/netcdf/examples/sresa1b_ncar_ccsm3-example.nc")

p2 = mv(p,string(tempname(),".nc"))

ds_nc = YAXArrayBase.to_dataset(p2)
vn = get_varnames(ds_nc)
@test sort(vn) == ["area", "lat", "lat_bnds", "lon", "lon_bnds", "msk_rgn",
 "plev", "pr", "tas", "time", "time_bnds", "ua"]
@test get_var_dims(ds_nc, "tas") == ["lon", "lat", "time"]
@test get_var_dims(ds_nc, "area") == ["lon", "lat"]
@test get_var_dims(ds_nc, "time") == ["time"]
@test get_var_dims(ds_nc, "time_bnds") == ["bnds", "time"]
@test get_var_attrs(ds_nc,"tas")["long_name"] == "air_temperature"
h = get_var_handle(ds_nc, "tas")
@test all(isapprox.(h[1:2,1:2], [215.893 217.168; 215.805 217.03]))
@test allow_parallel_write(ds_nc) == false
@test allow_missings(ds_nc) == false
end

@testset "Reading Zarr" begin
p = "gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/highresSST-present/r1i1p1f1/6hrPlev/psl/gn/v20170706/"
ds_zarr = to_dataset(p)
vn = get_varnames(ds_zarr)
@test sort(vn) == ["lat", "lat_bnds", "lon", "lon_bnds", "psl", "time", "time_bnds"]
@test get_var_dims(ds_zarr, "psl") == ["lon", "lat", "time"]
@test get_var_dims(ds_zarr, "time") == ["time"]
@test get_var_dims(ds_zarr, "time_bnds") == ["bnds", "time"]
@test get_var_attrs(ds_zarr,"psl")["long_name"] == "Sea Level Pressure"
h = get_var_handle(ds_zarr, "psl")
@test all(isapprox.(h[1:2,1:2,1], [99360.8  99334.9; 99360.8  99335.4]))
@test allow_parallel_write(ds_zarr) == true
@test allow_missings(ds_zarr) == true
end

function test_write(T)
    p = tempname()
    ds = create_empty(T, p)
add_var(ds, 0.5:1:9.5, "lon", ("lon",), Dict{String,Any}("units"=>"degrees_east"))
add_var(ds, 20:-1.0:1, "lat", ("lat",), Dict{String,Any}("units"=>"degrees_north"))
v = add_var(ds, Float32, "tas", (10,20), ("lon", "lat"), Dict{String,Any}("units"=>"Celsius"))

v[:,:] = collect(reshape(1:200, 10, 20))

@test sort(get_varnames(ds)) == ["lat","lon","tas"]
@test get_var_dims(ds, "tas") == ["lon", "lat"]
@test get_var_dims(ds, "lon") == ["lon"]
@test get_var_attrs(ds,"tas")["units"] == "Celsius"
h = get_var_handle(ds, "lon")
@test h[:] == 0.5:1:9.5
v = get_var_handle(ds, "tas")
@test v[1:2,1:2] == [1 11; 2 12]
end

@testset "Writing NetCDF" begin
  test_write(YAXArrayBase.NetCDFDataset)
end

@testset "Writing Zarr" begin
  test_write(YAXArrayBase.ZarrDataset)
end


text:
# for _atm, _bas, _env
const CHARGE_OF  = 1
const PTR_COORD  = 2
const NUC_MOD_OF = 3
const PTR_ZETA   = PTR_FRAC_CHARGE = 4
const ATM_SLOTS  = 6
const ATOM_OF    = 1
const ANG_OF     = 2
const NPRIM_OF   = 3
const NCTR_OF    = 4
const RADI_POWER = 4 # for ECP
const KAPPA_OF   = 5
const SO_TYPE_OF = 5 # for ECP
const PTR_EXP    = 6
const PTR_COEFF  = 7
const BAS_SLOTS  = 8
# pointer to env
const PTR_LIGHT_SPEED = 1
const PTR_COMMON_ORIG = 2
const PTR_RINV_ORIG   = 5
const PTR_RINV_ZETA   = 8
const PTR_RANGE_OMEGA = 9
const PTR_F12_ZETA    = 10
const PTR_GTG_ZETA    = 11
const AS_RINV_ORIG_ATOM = 18
const AS_ECPBAS_OFFSET = 19
const AS_NECPBAS      = 20
const PTR_ENV_START   = 21
# parameters from libcint
const NUC_POINT = 1
const NUC_GAUSS = 2
# nucleus with fractional charges. It can be used to mimic MM particles
const NUC_FRAC_CHARGE = 3
const NUC_ECP = 4  # atoms with pseudo potential

abstract type AbstractIntor end

struct CIntor{N,M} <: AbstractIntor
    atm::SVector{N, SVector{ATM_SLOTS, Cint}}
    bas::SVector{M, SVector{BAS_SLOTS, Cint}}
    env::Vector{Cdouble}
    shl::SVector{M, Int}
end

function CIntor(nuclei::Nuclei, basis::Basis)
    env = zeros(Cdouble, PTR_ENV_START-1)
    atm = make_atm!(nuclei, env)
    bas = make_bas!(nuclei, basis, env)
    shl = shell_size(nuclei, basis)
    N = length(atm)
    M = length(bas)
    return CIntor{N,M}(atm, bas, env, shl)
end

CIntor(mole::Mole) = CIntor(mole.nuclei, mole.basis)

function make_atm!(nuclei::Nuclei, env::Vector{Float64})
    N = length(nuclei)
    zeta = 0.0
    atm = Vector{SVector{ATM_SLOTS, Cint}}(undef, N)
    for i = 1:N
        nuc_charge = nuclei.numbers[i]
        ptr = length(env)
        _atm = zeros(Cint, ATM_SLOTS)
        _atm[CHARGE_OF] = nuc_charge
        _atm[PTR_COORD] = ptr
        _atm[NUC_MOD_OF] = NUC_POINT
        _atm[PTR_ZETA] = ptr + 3
        atm[i] = _atm
        _env = [nuclei.positions[i]; zeta]
        append!(env, _env)
    end
    return atm
end

function make_bas!(nuclei::Nuclei, basis::Basis, env::Vector{Float64})
    ptr_exp = []
    ptr_coeff = []
    for i = 1:length(basis.shells)
        _ptr_exp = []
        _ptr_coeff = []
        for shell in basis.shells[i]
            ptr = length(env)
            append!(env, shell.exponents)
            push!(_ptr_exp, ptr)
            for coeffs in shell.coefficients
                coeffs = normalizaton_sph(shell.angular_momentum, shell.exponents, coeffs)
                append!(env, coeffs)
            end
            push!(_ptr_coeff, ptr + length(shell.exponents))
        end
        push!(ptr_exp, _ptr_exp)
        push!(ptr_coeff, _ptr_coeff)
    end

    bas = SVector{BAS_SLOTS, Cint}[]
    for i = 1:length(nuclei)
        number = nuclei.numbers[i]
        shells_idx = findfirst(==(number), basis.numbers)
        shells = basis.shells[shells_idx]
        for (j, shell) in enumerate(shells)
            _bas = zeros(Cint, BAS_SLOTS)
            _bas[ATOM_OF] = i - 1
            _bas[ANG_OF] = shell.angular_momentum
            _bas[NPRIM_OF] = length(shell.exponents)
            _bas[NCTR_OF] = length(shell.coefficients)
            _bas[KAPPA_OF] = 0
            _bas[PTR_EXP] = ptr_exp[shells_idx][j]
            _bas[PTR_COEFF] = ptr_coeff[shells_idx][j]
            push!(bas, _bas)
        end
    end
    return bas
end

function shell_size(nuclei::Nuclei, basis::Basis)
    size = Int[]
    for i = 1:length(nuclei)
        number = nuclei.numbers[i]
        shells_idx = findfirst(==(number), basis.numbers)
        shells = basis.shells[shells_idx]
        for (j, shell) in enumerate(shells)
            push!(size, 2*shell.angular_momentum + 1)
        end
    end
    return size
end

function (intor::CIntor)(intor_name::String)
    cint_name = 'c'*intor_name*"_sph"
    func = getfield(Libcint, Symbol(cint_name))
    atm = reinterpret(Cint, Array(intor.atm))
    natm = Cint(length(intor.atm))
    bas = reinterpret(Cint, Array(intor.bas))
    nbas = Cint(length(intor.bas))
    if startswith(cint_name, "cint1e")
        return getints2c(func, atm, natm, bas, nbas, intor.env, intor.shl)
    elseif startswith(cint_name, "cint2e")
        return getints4c(func, atm, natm, bas, nbas, intor.env, intor.shl)
    else
        error("Integral name not found")
    end
end

function getints2c(
    func::Function,
    atm::AbstractVector{Cint},
    natm::Cint,
    bas::AbstractVector{Cint},
    nbas::Cint,
    env::AbstractVector{Cdouble},
    shl::AbstractVector{Int},
    )
    n = sum(shl)
    M = Array{Cdouble}(undef, n, n)
    MB = PseudoBlockArray(M, shl, shl)
    buf = Matrix{Cdouble}(undef, 64, 64)
    shls = Vector{Cint}(undef, 2)
    for i = 1:nbas, j = i:nbas
        shls .= (i-1,j-1)
        func(buf, shls, atm, natm, bas, nbas, env)
        mb = MB[Block(i,j)]
        for idx=1:length(mb)
            @inbounds mb[idx] = buf[idx]
        end
    end
    symmetrize!(M)
    return M
end

function getints4c(
    func::Function,
    atm::AbstractVector{Cint},
    natm::Cint,
    bas::AbstractVector{Cint},
    nbas::Cint,
    env::AbstractVector{Cdouble},
    shl::AbstractVector{Int},
    )
    n = sum(shl)
    T = Array{Cdouble}(undef, n, n, n, n)
    TB = PseudoBlockArray(T, shl, shl, shl, shl)
    buf = Matrix{Cdouble}(undef, 64, 64)
    shls = Vector{Cint}(undef, 4)
    for i = 1:nbas, j = i:nbas
        for k = i:nbas, l = k:nbas
            shls .= (i-1,j-1,k-1,l-1)
            func(buf, shls, atm, natm, bas, nbas, env)
            tb = TB[Block(i,j,k,l)]
            for idx=1:length(tb)
                @inbounds tb[idx] = buf[idx]
            end
        end
    end
    symmetrize!(T)
    return T
end

function symmetrize!(M::AbstractMatrix{Float64})
    n = LinearAlgebra.checksquare(M)
    for i = 1:n, j = i+1:n
        M[j, i] = M[i, j]
    end
end

function symmetrize!(T::AbstractArray{Float64,4})
    n = LinearAlgebra.checksquare(T)
    for i = 1:n, j = i:n, k = i:n, l = k:n
        T[i, j, l, k] = T[i, j, k, l]
        T[j, i, k, l] = T[i, j, k, l]
        T[j, i, l, k] = T[i, j, k, l]
        T[k, l, i, j] = T[i, j, k, l]
        T[k, l, j, i] = T[i, j, k, l]
        T[l, k, i, j] = T[i, j, k, l]
        T[l, k, j, i] = T[i, j, k, l]
    end
end

function normalizaton_sph(n::Int, es::AbstractVector{Float64}, d::AbstractVector{Float64})
    N = length(d)
    cs = Vector{Float64}(undef, N)
    for i=1:N
        cs[i] = d[i] * normalizaton_sph(n, es[i])
    end

    val = 0.0
    for i = 1:N
        val += (cs[i] / normalizaton_sph(n, es[i]))^2
        for j = i+1:N
            val += 2 * cs[i] * cs[j] / normalizaton_sph(n, (es[i] + es[j])/2)^2
        end
    end
    cs ./= val
    return cs
end

"""Normalized factor for GTO radial part"""
function normalizaton_sph(n::Int, a::Float64)
    return sqrt(2^(2*n+3) * factorial(n+1) * (2*a)^(n+1.5) /
                (factorial(2*n+2) * sqrt(π)))
end


text:
include("../roblib.jl")

function rect_border(r::Robot)

    while isborder(r, Left) == false || isborder(r, Down) == false # Перемещаем робота в юго-западный угол
        movements!(r, Down)
        movements!(r, Left)
    end

    num_steps = get_num_movements!(r, Right) # определяем размер поля

    max_S = 0
    rect_S = 0

    while isborder(r, Up) == false || isborder(r, Right) == false # пока не дойдёт до верхней границы поля
        for side in(Left, Right)
            if isborder(r, Up) == true && isborder(r, Right) == true # дошёл до верхней границы  поля
                break
            end
            rect_S = search_border(r, side, num_steps)
            if rect_S > max_S
                max_S = rect_S
            end
            if isborder(r, Up)==false
                move!(r, Up)
            end
        end
    end

    print(max_S)
end

function search_border(r::Robot,side::HorizonSide,numsteps::Int)

    num_borders = 0
    rect_S = 0

    num_steps = 0
    length = 0

    state=false

    while numsteps != 0 # сделано так , а не через for , т.к. при обходе прямоугольной нужно уменьшать количество шагов
        if isborder(r,side) == false
            move!(r,side)
            numsteps-=1
            if isborder(r, Up) == true
                state = true
                num_borders += 1
            else
                state=false
                if num_borders > 0
                    length = num_borders
                    rect_S = check(r,side, length)
                end
                num_borders = 0
            end
        else # обходит вертикальную перегородку
            while isborder(r, side) == true
                move!(r, Down)
                num_steps+=1
            end
            move!(r, side)
            numsteps -= 1
            while isborder(r, Up) == true # на случай, если мы обходим прямоугольную
                move!(r, side)
                numsteps -= 1
            end
            for _ in 1:num_steps
                move!(r, Up)
            end
            num_steps = 0
        end
    end

    return(rect_S)
end

function check(r::Robot,side::HorizonSide, length::Int) # проверяем , имеется ли перегородка , совмещённая с горизонтальной , если да, то значит мы нашли прямоугольную
    width = 0

    move!(r, Up)

    if isborder(r, inverse(side)) == true
        width = find_width(r,side)
        for _ in 1: width + 1
            move!(r, Down)
        end
    else
        move!(r, Down)
    end

    return(width * length)
end

function find_width(r::Robot, side::HorizonSide)
    width = 0

    while isborder(r,inverse(side))==true
        width += 1
        move!(r, Up)
    end

    return(width)
end


text:
convert(::Type{EventTraceData}, C::EventChannel) = EventTraceData(C)

function convert(::Type{EventTraceData}, S::SeisData)
  TD = EventTraceData(getfield(S, :n))
  for f in datafields
    if (f in unindexed_fields) == false
      setfield!(TD, f, deepcopy(getfield(S, f)))
    end
  end
  return TD
end

function convert(::Type{SeisData}, TD::EventTraceData)
  S = SeisData(getfield(TD, :n))
  for f in datafields
    if (f in unindexed_fields) == false
      setfield!(S, f, deepcopy(getfield(TD, f)))
    end
  end
  return S
end

function convert(::Type{SeisChannel}, D::EventChannel)
  C = SeisChannel()
  for f in datafields
    setfield!(C, f, deepcopy(getfield(D, f)))
  end
  return C
end

function convert(::Type{EventChannel}, C::SeisChannel)
  D = EventChannel()
  for f in datafields
    setfield!(D, f, deepcopy(getfield(C, f)))
  end
  return D
end

function unsafe_convert(::Type{SeisData}, TD::EventTraceData)
  S = SeisData(getfield(TD, :n))
  for f in datafields
    if (f in unindexed_fields) == false
      setfield!(S, f, getfield(TD, f))
    end
  end
  return S
end

function unsafe_convert(::Type{EventTraceData}, S::SeisData)
  TD = EventTraceData(getfield(S, :n))
  for f in datafields
    if (f in unindexed_fields) == false
      setfield!(TD, f, getfield(S, f))
    end
  end
  return TD
end

push!(TD::EventTraceData, C::SeisChannel) = push!(TD, convert(EventChannel, C))


text:
using CSV, DataFrames, Plots, StatsPlots

export plot_trade_union_membership, plot_GDP_per_capita, plot_maddison_GDP_data
export plot_GDP
export north_vs_south_countries, anglo_germanic_countries
export plot_relative_GDP

ENV["GKS_ENCODING"]="utf8"  # Fix some annoying warning when plotting

# Find names of colors here: http://juliagraphics.github.io/Colors.jl/stable/namedcolors/

function plot_trade_union_membership()
    path = joinpath(datadir, "trade-union-membership.csv")
    unions = CSV.File(path) |> DataFrame
    unions = select(unions, :Country, :Time, :Value)
    country_unions = groupby(unions, :Country)

    country_names = ["Norway", "Sweden", "Denmark", "Finland",
                 "United Kingdom", "United States", "Canada"]
    colors = [:purple, :blue, :cyan, :green, :orange, :red, :brown]
    plot()
    for (i, cname) in enumerate(country_names)
        country = country_unions[(Country=cname,)]
        @df country plot!(:Time, :Value, label=cname, lw=2, seriescolor=colors[i])
    end
    plot!(minorgrid=true,      # turn on minor grid lines
          gridalpha=0.3,       # Make major grid lines a bit more visible
          legend=:outertopleft # Place box with legends outside of plot on top left
    )
end

function plot_GDP_per_capita()
    gdp = load_GDP_per_capita()
    p = @df gdp plot(
           :year,
           [:norway, :sweden, :denmark, :finland, :uk, :usa],
           lw=2,
           minorgrid=true,
           gridalpha=0.3,
           legend=:bottomright,
           label=["norway"  "sweden"  "denmark"  "finland"  "uk"  "usa"],
           seriescolor=[:magenta4  :mediumblue  :deeppink2 :purple4 :lime :green]
           # seriescolor=[:aqua  :teal  :dodgerblue :blue :gold :orange :orangered3]
        )
    # p = plot(gdp.year,   # x-axis
    #         [gdp.norway, gdp.sweden, gdp.finland, gdp.denmark, gdp.uk, gdp.usa, gdp.canada],
    #         lw=2,        # line width
    #         legend=:bottomright, # where to put legend
    #         label=["norway" "sweden" "finland" "denmark" "uk" "usa" "canada"] # labels in legend
    #     )
    # savefig(p, "nordic-GDP-growth-from-60s.svg")
end

function plot_maddison_GDP_data()
    gdp = load_maddison_GDP_data()
    country_gdp = groupby(gdp, :Country)
    plot_GDP(country_gdp)
end

# This is a smarter way of doing the plotting
function simpler_maddison_GDP_plot()
    gdp = load_maddison_GDP_data()
    countries = DataFrame(Country=["Norway", "Sweden", "Denmark", "Finland", "United Kingdom", "United States"])

    # Like a set intersection. Match on :Country column
    sample = innerjoin(gdp, countries, on=:Country)

    # Filter out irrelevant years
    sample = filter(row -> 1950 < row.Year < 2020, sample)

    # Creates a separate series for every unique value in value provided to
    # group property. Every row with that group value gets assigned to the same series
    plot(sample.Year, sample.GDP, group=sample.Country)

    # You can also write this as:
    @df sample plot(:Year, :GDP, group=:Country)
end

const north_vs_south_countries = [
        "Norway"    => :purple,
        "Sweden"    => :cyan,
        "Finland"   => :blue,
        "Portugal"  => :goldenrod2,
        "Argentina" => :orange,
        "Spain"     => :brown,
        "Italy"     => :red,
        "Venezuela" => :orangered3,
      ]

const anglo_germanic_countries = [
      "Denmark"    => :purple,
      "Germany"    => :cyan,
      "Austria"   => :blue,
      "United Kingdom"  => :goldenrod2,
      "United States" => :orange,
      "Canada"     => :brown,
      "Italy"     => :red,
      "France" => :orangered3,
    ]

"""
    plot_GDP(country_gdp::GroupedDataFrame, range)

Plots GDP data for the year `range` given.
"""
function plot_GDP(country_gdp::GroupedDataFrame, range = 1800:2020, countries = north_vs_south_countries)
    plot()
    for (cname, color) in countries
        country = country_gdp[(Country=cname,)]
        country = filter(row -> row.Year in range, country)
        @df country plot!(:Year, :GDP, label=cname, lw=2, seriescolor=color)
        # @df country plot!(:Year, :GDP, label=cname, lw=1)
    end
    plot!(minorgrid=true,      # turn on minor grid lines
          gridalpha=0.3,       # Make major grid lines a bit more visible
          # legend=:outertopleft # Place box with legends outside of plot on top left
          legend=:topleft

    )
end

function plot_relative_GDP(base_country::DataFrame, country_gdp::GroupedDataFrame, range = 1900:2020, countries = anglo_germanic_countries)
    bcountry = base_country = filter(row -> row.Year in range, base_country)
    plot(ylims = (40, 130))
    for (cname, color) in countries
        country = country_gdp[(Country=cname,)]
        country = filter(row -> row.Year in range, country)
        country.rGDP = 100*(country.GDP ./ bcountry.GDP)
        @df country plot!(:Year, :rGDP, label=cname, lw=2, seriescolor=color)
        # @df country plot!(:Year, :GDP, label=cname, lw=1)
    end
    plot!(minorgrid=true,      # turn on minor grid lines
          gridalpha=0.3,       # Make major grid lines a bit more visible
          # legend=:outertopleft # Place box with legends outside of plot on top left
          legend=:topleft

    )
end

function plot_cultural_diff()
    path = joinpath(datadir, "geerthofstede-2015-08-16.csv")
    cultural = CSV.File(path) |> DataFrame
    cultures = select(cultural, :country,
        :pdi => :hierarchy,
        :idv => :individualism,
        :mas => :masculinity,
        :uai => :uncertainty,
        :ltowvs => :long_term,
        :ivr => :indulgence)

    # Countries we want to filter on by performing an innerjoin
    # countries = DataFrame(country=["Norway", "Sweden", "Denmark", "Finland", "U.S.A."])
    countries = DataFrame(country=["Norway", "Sweden", "Denmark", "Japan", "China"])

    # Pick cultures where the :country column matches in both DataFrame objects
    sample = innerjoin(cultures, countries, on = :country)

    # Stack turns column names into row values. E.g. if you got the row:
    # country hierarchy individualism
    # China   80        20
    #
    # then this becomes:
    # country variable       value
    # China   hierarchy      80
    # China   individualism  20
    #
    # So every column entry gets a separate row with the name of that column
    data = stack(sample, 2:7)

    # Why is this useful? Because the `group` attribute allows you to bundle values
    # having identical column value such as "hierarchy" into a seprate series.
    # Here each series is one country. The x-axis values for each country is the
    # cultural dimensions under the `variable` column
    groupedbar(data.variable, data.value, group=data.country)

    # Turn into CSV viewable in Numbers
    # clipboard(repr(MIME("text/csv"), sample))
end

max_stars_repo_name: ordovician/NordicEcon

max_issues_repo_name: ordovician/NordicEcon

max_forks_repo_name: ordovician/NordicEcon


text:
#addprocs(2)
#@show nprocs()

using RoME
using Test

# might be unnecessary
using JLD2  #, HDF5

println("[TEST] numeric root...")
include("testhigherdimroots.jl")
println("[SUCCESS]")

include("testManifoldsPose2Equivalent.jl")

# Requires standardized testing strategy
# println("[TEST] Camera function evaluations...")
# include("testCameraFunctions.jl")
# println("[SUCCESS]")

@warn "Skipping multiple feature constraint test for the time being"
# println("[TEST] MultipleFeatures constraints")
# include("testmultiplefeatures.jl")
# println("[SUCCESS]")

println("[TEST] Linear array function evaluations...")
include("testDidsonFunctions.jl")
println("[SUCCESS]")

include("testPoint2Point2.jl")

include("testPoint2Point2Init.jl")

include("testBasicPose2Stationary.jl")

println("[TEST] Pose2 evaluations...")
include("TestPoseAndPoint2Constraints.jl")
println("[SUCCESS]")

include("testBeehive2D_CliqByCliq.jl")

# include("HexagonalLightGraphs.jl")

@testset "[TEST] Pose2 evaluations..." begin
  include("testDynPoint2D.jl")
end

include("testPoint2Point2WorldBearing.jl")

include("testBearingRange2D.jl")

include("testDeltaOdo.jl")

include("testFixedLagFG.jl")

include("testMultimodalRangeBearing.jl")

include("testDynPose2D.jl")

println("[TEST] Pose3 evaluations...")
include("threeDimLinearProductTest.jl")
println("[SUCCESS]")

println("[TEST] ensure Pose3Pose3NH evaluations...")
include("testPose3Pose3NH.jl")
println("[SUCCESS]")

@warn "Skipping JLD tests in RoME"
# println("[TEST] saving to and loading from .jld2 file")
# savejld(fg, file="tempfg.jld2" )
# fgu = loadjld( file="tempfg.jld2" )
# Base.rm("tempfg.jld2")
# println("Success")

println("[TEST] PartialPose3XYYaw evaluations...")
include("testPartialXYH.jl")
println("[SUCCESS]")

println("[TEST] partial pose3 evaluations...")
include("testpartialpose3.jl")
println("[SUCCESS]")

println("[TEST] packing converters...")
include("testpackingconverters.jl")
println("[SUCCESS]")

include("TestDefaultFGInitialization.jl")

include("testAccumulateFactors.jl")

include("testDeadReckoningTether.jl")


text:
using Thrift
using Test

import Thrift.process, Thrift.meta

# include the generated module, which in-turn includes our implementation code in `proto_tests_impl.jl`
# prepend with @everywhere when running a TProcessPoolServer on multiple Julia processors
@isdefined(srvcctrl) || include("gen-jl/srvcctrl/srvcctrl.jl");
@isdefined(proto_tests) || include("gen-jl/proto_tests/proto_tests.jl");

import .proto_tests: ProtoTestsClient, InvalidOperation, AllTypes, AllTypesDefault, TestEnum
import .proto_tests: test_hello, test_exception, test_oneway, ping, test_enum, test_types, test_types_default
import .srvcctrl: start_service, stop_service

# create a client instance with our choice of protocol and transport
clnt_transport = TSocket(19999)
#proto = TCompactProtocol(clnt_transport)
proto = TBinaryProtocol(clnt_transport)
clnt = ProtoTestsClient(proto)

function run_client()
    println("Transport: $(typeof(clnt_transport))")
    println("Protocol : $(typeof(proto))")
    println("Starting requests...")

    # open a connection
    open(clnt_transport)

    # invoke service and print the result
    println("\nCalling test_hello...")
    ret = test_hello(clnt, "Julia")
    println(ret)
    @test endswith(ret, "Julia")

    println("\nCalling test_exception...")
    try
        test_exception(clnt)
    catch ex
        println(ex)
        @test isa(ex, InvalidOperation)
        @test ex.oper == "test_exception"
    end

    println("\nCalling ping...")
    ping(clnt)

    println("\nCalling test_oneway...")
    test_oneway(clnt)

    println("\nCalling test_types_default...")
    ret = test_types_default(clnt, AllTypesDefault())
    println(ret)
    @test get_field(ret, :bool_val) == false
    @test get_field(ret, :byte_val) == 0x02
    @test get_field(ret, :i16_val) == 11
    @test get_field(ret, :i32_val) == 21
    @test get_field(ret, :i64_val) == 31
    @test get_field(ret, :double_val) == -10.1
    @test get_field(ret, :string_val) == "HELLO WORLD"

    v = get_field(ret, :map_val)
    @test length(v) == 2
    @test v[Int32(2)] == 40
    @test v[Int32(1)] == 20

    @test get_field(ret, :list_val) == Int16[11,12,13]
    v = get_field(ret, :set_val)
    @test length(v) == 3
    @test UInt8(13) in v
    @test UInt8(14) in v
    @test UInt8(15) in v

    println("\nCalling stop_service...")
    stop_service(clnt)

    # close connection
    close(clnt_transport)
end

println("Waiting for server to start...")
sleep(2)
run_client()
sleep(2)


text:
using Test
import Pkg

using InteractiveCodeSearch

macro test_nothrow(ex)
    quote
        @test begin
            $(esc(ex))
            true
        end
    end
end


text:
using DFControl, Test

testdir = joinpath(dirname(dirname(pathof(DFControl))), "test")
testjobpath = joinpath(testdir, "testassets", "test_job")

name = "Pt"
local_dir = testjobpath
server_dir = testjobpath
bin_dir = joinpath(homedir(), "bin")
excs = [Exec("mpirun", bin_dir, :np => 24), Exec("pw.x", bin_dir, Dict{Symbol, Any}(:nk => 2))]

pseudoset = :test

header = ["#SBATCH -N 1", "#SBATCH --ntasks-per-node=24",
          "#SBATCH --time=24:00:00", "#SBATCH -p defpart",
          "module load open-mpi/gcc/1.10.4-hfi", "module load mkl/2016.1.056"
         ]

scf_data = Dict(:kpoints => [6, 6, 6, 1, 1, 1], :flags => [:verbosity => "low"])
bands_data = Dict(:kpoints => [[0.5, 0.5, 0.5, 100.],
                                [0.0, 0.0, 0.0, 100.],
                                [0.0, 0.5, 0.0, 1.]],
                  :flags => [:verbosity => "high", :nbnd => 8])

nscf_data = merge(bands_data, Dict(:kpoints => [10, 10, 10]))
calculations = [:vc_relax => (excs, scf_data), :scf => (excs, scf_data), :bands => (excs, bands_data), :nscf =>(excs, nscf_data)]

job = DFJob(name, local_dir, joinpath(testjobpath, "Pt.cif"), calculations,
      :prefix       => "$name",
      :restart_mode => "from_scratch",
      :mixing_mode  => "plain",
      :mixing_beta  => 0.7,
      :conv_thr     => 1.0e-8,
      #kwargs
      server_dir  = server_dir,
      bin_dir     = bin_dir,
      header      = header,
      pseudoset  = :test,
     )

save(job)
show(job)
@test data(job["scf"], :k_points).data == [6, 6, 6, 1, 1, 1]
@test data(job, "nscf", :k_points).data == DFControl.kgrid(10, 10, 10, :nscf)
@test all(values(job[:ecutwfc]) .== 32.0)
@test job["scf"][:prefix] == job["nscf"][:prefix] == "$name"
@test job["bands"][:verbosity] == "high"

setflags!(job, :prefix => "blabla", print=false)
@test job["scf"][:prefix] == job["nscf"][:prefix] == "blabla"
setflags!(job, :Hubbard_U => [4], print=false)
setflags!(job, :Hubbard_J => [4 4 5], print=false)
@test job["scf"][:Hubbard_U] == job["nscf"][:Hubbard_U] == [4.0]
@test job["scf"][:Hubbard_J] == job["nscf"][:Hubbard_J] == [4.0  4.0  5.0]

struct2 = DFControl.create_supercell(job.structure, 1, 2, 1)
newpositions = [DFControl.position_cart(at) for at in atoms(struct2)]
oldposition = atoms(job.structure)[1].position_cart
cell_ = DFControl.cell(job.structure)

@test atoms(job.structure)[1].position_cart == job.structure.cell' * atoms(job.structure)[1].position_cryst
@test oldposition == newpositions[1]
@test oldposition + cell_[:, 1] ∈ newpositions
@test oldposition + 2*cell_[:, 2] ∈ newpositions
@test oldposition + 1*cell_[:, 3] ∈ newpositions
@test oldposition + (1*cell_[:, 3] + 1*cell_[:, 2]) ∈ newpositions


text:
module SmuUtilities

include("utilities.jl")

end  #  module SmuUtilities

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
module RunOut

using EzXML
using HTTP
using JSON
using Parameters
using TimeZones

export TrackArtist, Contributor, Rating, Community, Company
export Identifier, Format, Image, Label, Track, Video, Release
export Client, fetch_release, fetch_artist
export generate_xml, save_xml

include("models.jl")
include("parse.jl")
include("discogs.jl")
include("xml.jl")

end


text:
#GroupContribution Struct.
#it statically stores the id and values.
struct GroupContribution{N}
    id::NTuple{N,Int}
    value::NTuple{N,Int}
end

#==
Property Struct:
it holds a value and its metadata.
for now, the metadata is the long name and the units.
the property is parametrically typed on the value, to recognize if is a
string, Int, Float, VariableProperty, or GroupContribution.
it can replace IntVal, StringVal,
TODO:in the future, using Unitful to hold the value, something like:
Property{T,UNITS} where UNITS <: Unitful.Unit
==#

struct Property{T}
    name::String
    units::String
    value::T
end

value(prop::Property) = prop.value
name(prop::Property) = prop.name
units(prop::Property) = prop.units

#==
#those structs are no longer used in this branch.

struct ScalarProperty
    units::String
    value::Float64
end

struct RangeBoundary
    units::String
    value::Float64
end

struct VectorProperty
    units::String
    EqNo::Int
    A::Float64
    B::Float64
    C::Float64
    D::Float64
    E::Float64
    Tmin::RangeBoundary
    Tmax::RangeBoundary
end

struct PureComp
    LibraryIndex::Int
    CompoundID::String
    StructureFormula::String
    Family::Int
    CriticalTemperature::ScalarProperty
    CriticalPressure::ScalarProperty
    CriticalVolume::ScalarProperty
    CriticalComperssibility::ScalarProperty
    NormalBoilingPointTemperature::ScalarProperty
    NormalMeltingPointTemperature::ScalarProperty
    TripplePointTemperature::ScalarProperty
    TripplePointPressure::ScalarProperty
    MolecularWeight::ScalarProperty
    LiquidVolumeAtNormalBoilingPoint::ScalarProperty
    AcentricityFactor::ScalarProperty
    SolubilityParameter::ScalarProperty
    DipoleMoment::ScalarProperty
    HeatOfFormation::ScalarProperty
    GibbsEnergyOfFormation::ScalarProperty
    AbsEntropy::ScalarProperty
    HeatOfFusionAtMeltingPoint::ScalarProperty
    HeatOfCombustion::ScalarProperty
    LiquidDensity::VectorProperty
    VaporPressure::VectorProperty
    HeatOfVaporization::VectorProperty
    LiquidHeatCapacityCp::VectorProperty
    IdealHeatCapacityCp::VectorProperty
    SecondVirialCoefficient::VectorProperty
    LiquidViscosity::VectorProperty
    VaporViscosity::VectorProperty
    LiquidThermalConductivity::VectorProperty
    VaporThermalConductivity::VectorProperty
    RPPHeatCapacityCp::VectorProperty
    RelativeStaticPermittivity::VectorProperty
    AntoineVaporPressure::VectorProperty
    LiquidViscosityRPS::VectorProperty
    COSTALDVolume::ScalarProperty
    DiameterLJ::ScalarProperty
    EnergyLJ::ScalarProperty
    RacketParameter::ScalarProperty
    FullerVolume::ScalarProperty
    Parachor::ScalarProperty
    SpecificGravity::ScalarProperty
    Charge::ScalarProperty
    CostaldAcentricFactor::ScalarProperty
    WilsonVolume::ScalarProperty
    ChaoSeaderAcentricFactor::ScalarProperty
    ChaoSeaderSolubilityParameter::ScalarProperty
    ChaoSeaderLiquidVolume::ScalarProperty
    CAS::String
end
==#