text:
using PyPlot, Revise, HiQGA.transD_GP, Random, SpecialFunctions
## model
zfixed   = [-1e5,   0,    ]
rho      = [1e20,   1,  ]
nmax = 200
##  geometry
rTx = 15
rRx = 0.001
zRx = -0.001
zTx = -0.002
doconvramp = false
modelprimary = true
nkᵣeval = 200
ntimesperdecade = 15
nfreqsperdecade = 15
times = 10 .^LinRange(-6,1,40)
## build operator
F = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      ntimesperdecade = ntimesperdecade,
                      nfreqsperdecade = nfreqsperdecade,
                      times  = times,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx,
                      doconvramp = doconvramp,
                      modelprimary = modelprimary,
                      nkᵣeval = nkᵣeval)
## get TD field
transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho)
## Analytic loop
μ = transD_GP.AEM_VMD_HMD.μ
sig0 = 1/rho[2]
theta = sqrt.(μ*sig0./4F.interptimes)
a     = rTx
dbzdt = -1 ./(sig0*a^3)*(3*erf.(theta*a) - 2/sqrt(pi)*theta.*a.*(3 .+ 2*theta.^2*a^2).*exp.(-theta.^2*a^2))
dbzdt = dbzdt/(pi*a^2)
##
figure()
s1 = subplot(211)
loglog(F.interptimes,μ*abs.(F.HTD_z_interp), label="loop W&H")
loglog(F.interptimes,abs.(dbzdt), "--", label="loop W&H transient half space ")
xlim(extrema(times))
ylabel("dBzdt V/Am^4")
grid(which="both")
legend()
s2 = subplot(212, sharex=s1)
# loglog(F.interptimes, 100*abs.(-AEM_VMD_HMD.μ*F.HTDinterp-dbzdt)./abs.(dbzdt))
# signs are supposed to be the same at loop centre ...
loglog(F.interptimes, 100*abs.(abs.(μ*F.HTD_z_interp)-abs.(dbzdt))./abs.(dbzdt))
grid(which="both")
ylabel("% difference")
xlabel("time s")

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
#-----------------------------------------------------------------------# NBClassifier
"""
    NBClassifier(p::Int, T::Type; stat = Hist(15))

Calculate a naive bayes classifier for classes of type `T` and `p` predictors.  For each
class `K`, predictor variables are summarized by the `stat`.

# Example

    x, y = randn(10^4, 10), rand(Bool, 10^4)

    o = fit!(NBClassifier(10, Bool), (x,y))
    collect(keys(o))
    probs(o)

    xi = randn(10)
    predict(o, xi)
    classify(o, xi)
"""
mutable struct NBClassifier{T, G<:Group} <: OnlineStat{XY}
    d::OrderedDict{T, G}
    init::G
    # For trees
    id::Int
    j::Int
    at::Union{Number, String, Symbol, Char}
    ig::Float64
    children::Vector{Int}
end
function NBClassifier{T}(g::G; id=1) where {T,G}
    NBClassifier{T,G}(OrderedDict{T,G}(), g, id, 0, -Inf, -Inf, Int[])
end
function NBClassifier(g::G, T::Type; id=1) where {G<:Group}
    NBClassifier(OrderedDict{T, G}(), g, id, 0, -Inf, -Inf, Int[])
end
function NBClassifier(p::Int, T::Type; id=1, stat=KHist(10))
    NBClassifier(p * stat, T; id=id)
end
function Base.show(io::IO, o::NBClassifier)
    print(io, "NBClassifier")
    for (k, p) in zip(keys(o), probs(o))
        print(io, "\n    > $k ($(round(p, digits=4)))")
    end
end
function _fit!(o::NBClassifier, xy)
    x, y = xy
    if haskey(o.d, y)
        _fit!(o.d[y], x)
    else
        stat = copy(o.init)
        _fit!(stat, x)
        o.d[y] = stat
    end
end
_merge!(o::NBClassifier, o2::NBClassifier) = merge!(merge!, o.d, o2.d)
Base.getindex(o::NBClassifier, j) = [stat[j] for stat in values(o.d)]

Base.keys(o::NBClassifier) = keys(o.d)
Base.values(o::NBClassifier) = values(o.d)
nkeys(o::NBClassifier) = length(o.d)
nvars(o::NBClassifier) = length(o.init)
nobs(o::NBClassifier) = isempty(o.d) ? 0 : sum(nobs, values(o))
probs(o::NBClassifier) = isempty(o.d) ? zeros(0) : map(nobs, values(o)) ./ nobs(o)

function _predict(o::NBClassifier, x::VectorOb, p = zeros(nkeys(o)), n = nobs(o))
    for (k, gk) in enumerate(values(o))
        # P(Ck)
        p[k] = log(nobs(gk) / n + ϵ)
        # P(xj | Ck)
        for j in 1:length(x)
            p[k] += log(pdf(gk[j], x[j]) + ϵ)
        end
        p[k] = exp(p[k])
    end
    sp = sum(p)
    sp == 0.0 ? p : rmul!(p, inv(sp))
end
function _classify(o::NBClassifier, x::VectorOb, p = zeros(nkeys(o)), n = nobs(o))
    _, k = findmax(_predict(o, x, p, n))
    index_to_key(o, k)
end
function index_to_key(d, i)
    for (k, ky) in enumerate(keys(d))
        k == i && return ky
    end
end

predict(o::NBClassifier, x::VectorOb) = _predict(o, x)
predict(o::NBClassifier, x) = [predict(o, xi) for xi in x]
predict(o::NBClassifier, x::AbstractMatrix) = predict(o, OnlineStatsBase.eachrow(x))

classify(o::NBClassifier, x::VectorOb) = _classify(o, x)
classify(o::NBClassifier, x) = [classify(o, xi) for xi in x]
classify(o::NBClassifier, x::AbstractMatrix) = classify(o, OnlineStatsBase.eachrow(x))

# Tree stuff

function classify_node(o::NBClassifier)
    _, k = findmax([nobs(g) for g in values(o)])
    index_to_key(o, k)
end

function split!(o::NBClassifier)
    nroot = [nobs(g) for g in values(o)]
    nleft, nright = copy(nroot), copy(nroot)
    entropy_root = entropy(o)
end

entropy(o::NBClassifier) = entropy(probs(o), 2)

#  function split(o::NBClassifier)
#     nroot = [nobs(g) for g in values(o)]
#     nleft = copy(nroot)
#     nright = copy(nroot)
#     split = NBSplit(length(nroot))
#     entropy_root = entropy(o)
#     for j in 1:nvars(o)
#         ss = o[j]
#         stat = merge(ss)
#         for loc in split_candidates(stat)
#             for k in 1:nkeys(o)
#                 nleft[k] = round(Int, n_sent_left(ss[k], loc))
#             end
#             entropy_left = entropy(nleft ./ sum(nleft))
#             @. nright = nroot - nleft
#             entropy_right = entropy(nright ./ sum(nright))
#             entropy_after = smooth(entropy_right, entropy_left, sum(nleft) / sum(nroot))
#             ig = entropy_root - entropy_after
#             if ig > split.ig
#                 split.j = j
#                 split.at = loc
#                 split.ig = ig
#                 split.nleft .= nleft
#             end
#         end
#     end
#     left = NBClassifier(collect(keys(o)), o.init)
#     right = NBClassifier(collect(keys(o)), o.init)
#     for (i, g) in enumerate(values(left.d))
#         g.nobs = split.nleft[i]
#     end
#     for (i, g) in enumerate(values(right.d))
#         g.nobs = nroot[i] - split.nleft[i]
#     end
#     o, split, left, right
# end

# #-----------------------------------------------------------------------# NBClassifier
# """
#     NBClassifier(group::Group, labeltype::Type)

# Create a naive bayes classifier, using the stats in `group` to approximate the
# distributions of each predictor variable conditioned on label.

# - For continuous variables, use [`Hist(nbin)`](@ref).
# - For categorical variables, use [`CountMap(T)`](@ref).

# # Example

#     x = randn(10^5, 10)
#     y = rand(1:5, 10^5)
#     o = NBClassifier(10Hist(20), Float64)
#     series((x, y), o)
#     predict(o, x)
#     classify(o, x)
# """
# struct NBClassifier{T, G <: Group} <: ExactStat{(1, 0)}
#     d::OrderedDict{T, G}  # class => group
#     init::G        # empty group
# end
# NBClassifier(T::Type, g::G) where {G<:Group} = NBClassifier(OrderedDict{T,G}(), g)
# NBClassifier(g::Group, T::Type) = NBClassifier(T, g)
# function NBClassifier(labels::Vector{T}, g::G) where {T, G<:Group}
#     NBClassifier(OrderedDict{T, G}(lab=>copy(g) for lab in labels), g)
# end
# NBClassifier(p::Int, T::Type, b=20) = NBClassifier(T, p * Hist(b))

# function Base.show(io::IO, o::NBClassifier)
#     print(io, name(o))
#     sd = sort(o.d)
#     for di in sd
#         print(io, "\n    > ", first(di), " (", round(nobs(last(di)) / nobs(o), 4), ")")
#     end
# end

# Base.keys(o::NBClassifier) = keys(o.d)
# Base.values(o::NBClassifier) = values(o.d)
# Base.haskey(o::NBClassifier, y) = haskey(o.d, y)
# nvars(o::NBClassifier) = length(o.init)
# nkeys(o::NBClassifier) = length(o.d)
# nobs(o::NBClassifier) = sum(nobs, values(o))
# probs(o::NBClassifier) = [nobs(g) for g in values(o)] ./ nobs(o)
# Base.getindex(o::NBClassifier, j) = [stat[j] for stat in values(o)]

# # d is an object that iterates keys in known order
# function index_to_key(d, i)
#     for (k, ky) in enumerate(keys(d))
#         k == i && return ky
#     end
# end

# function fit!(o::NBClassifier, xy, γ)
#     x, y = xy
#     if haskey(o, y)
#         g = o.d[y]
#         fit!(g, x, 1 / (nobs(g) + 1))
#     else
#         o.d[y] = fit!(copy(o.init), x, 1.0)
#     end
# end
# entropy(o::NBClassifier) = entropy(probs(o), 2)

# function _predict(o::NBClassifier, x::VectorOb, p = zeros(nkeys(o)), n = nobs(o))
#     for (k, gk) in enumerate(values(o))
#         # P(Ck)
#         p[k] = log(nobs(gk) / n + ϵ)
#         # P(xj | Ck)
#         for j in 1:length(x)
#             p[k] += log(pdf(gk[j], x[j]) + ϵ)
#         end
#         p[k] = exp(p[k])
#     end
#     sp = sum(p)
#     sp == 0.0 ? p : p ./= sum(p)
# end
# function _classify(o::NBClassifier, x::VectorOb, p = zeros(nkeys(o)), n = nobs(o))
#     _, k = findmax(_predict(o, x, p, n))
#     index_to_key(o, k)
# end
# predict(o::NBClassifier, x::VectorOb) = _predict(o, x)
# classify(o::NBClassifier, x::VectorOb) = _classify(o, x)
# function classify_node(o::NBClassifier)
#     _, k = findmax([nobs(g) for g in values(o)])
#     index_to_key(o, k)
# end
# for f in [:(_predict), :(_classify)]
#     @eval begin
#         function $f(o::NBClassifier, x::AbstractMatrix, ::Rows = Rows())
#             n = nobs(o)
#             p = zeros(nkeys(o))
#             mapslices(xi -> $f(o, xi, p, n), x, 2)
#         end
#         function $f(o::NBClassifier, x::AbstractMatrix, ::Cols)
#             n = nobs(o)
#             p = zeros(nkeys(o))
#             mapslices(xi -> $f(o, xi, p, n), x, 1)
#         end
#     end
# end

# function split(o::NBClassifier)
#     nroot = [nobs(g) for g in values(o)]
#     nleft = copy(nroot)
#     nright = copy(nroot)
#     split = NBSplit(length(nroot))
#     entropy_root = entropy(o)
#     for j in 1:nvars(o)
#         ss = o[j]
#         stat = merge(ss)
#         for loc in split_candidates(stat)
#             for k in 1:nkeys(o)
#                 nleft[k] = round(Int, n_sent_left(ss[k], loc))
#             end
#             entropy_left = entropy(nleft ./ sum(nleft))
#             @. nright = nroot - nleft
#             entropy_right = entropy(nright ./ sum(nright))
#             entropy_after = smooth(entropy_right, entropy_left, sum(nleft) / sum(nroot))
#             ig = entropy_root - entropy_after
#             if ig > split.ig
#                 split.j = j
#                 split.at = loc
#                 split.ig = ig
#                 split.nleft .= nleft
#             end
#         end
#     end
#     left = NBClassifier(collect(keys(o)), o.init)
#     right = NBClassifier(collect(keys(o)), o.init)
#     for (i, g) in enumerate(values(left.d))
#         g.nobs = split.nleft[i]
#     end
#     for (i, g) in enumerate(values(right.d))
#         g.nobs = nroot[i] - split.nleft[i]
#     end
#     o, split, left, right
# end

# n_sent_left(o::Union{OrderStats, Hist}, loc) = sum(o, loc)
# n_sent_left(o::CountMap, label) = o[label]

# #-----------------------------------------------------------------------# NBSplit
# # Continuous:  x[j] < at
# # Categorical: x[j] == at
# mutable struct NBSplit{}
#     j::Int
#     at::Any
#     ig::Float64
#     nleft::Vector{Int}
# end
# NBSplit(n=0) = NBSplit(0, -Inf, -Inf, zeros(Int, n))

# whichchild(o::NBSplit, x) = x[o.j] < o.at ? 1 : 2

# #-----------------------------------------------------------------------# NBNode
# mutable struct NBNode{T <: NBClassifier} <: ExactStat{(1, 0)}
#     nbc::T
#     id::Int
#     parent::Int
#     children::Vector{Int}
#     split::NBSplit
# end
# function NBNode(o::NBClassifier; id = 1, parent = 0, children = Int[], split = NBSplit())
#     NBNode(o, id, parent, children, split)
# end
# function Base.show(io::IO, o::NBNode)
#     print(io, "NBNode ", o.id)
#     if o.split.j > 0
#         print(io, " (split on $(o.split.j)")
#     end
# end

# #-----------------------------------------------------------------------# NBTree
# """
#     NBTree(o::NBClassifier; maxsize=5000, splitsize=1000)

# Create a decision tree where each node is a naive bayes classifier.  A node will split
# when it reaches `splitsize` observations and no more splits will occur once `maxsize`
# nodes are in the tree.

# # Example

#     x = randn(10^5, 10)
#     y = rand(Bool, 10^5)
#     o = NBTree(NBClassifier(10Hist(20), Bool))
#     series((x,y), o)
#     classify(o, x)
# """
# mutable struct NBTree{T<:NBNode} <: ExactStat{(1, 0)}
#     tree::Vector{T}
#     maxsize::Int
#     splitsize::Int
# end
# function NBTree(o::NBClassifier; maxsize = 5000, splitsize = 1000)
#     NBTree([NBNode(o)], maxsize, splitsize)
# end
# NBTree(args...; kw...) = NBTree(NBClassifier(args...); kw...)
# function Base.show(io::IO, o::NBTree)
#     print(io, "NBTree(size = $(length(o.tree)), splitsize=$(o.splitsize))")
# end

# function fit!(o::NBTree, xy, γ)
#     x, y = xy
#     i, node = whichleaf(o, x)
#     fit!(node.nbc, xy, γ)
#     if length(o.tree) < o.maxsize && nobs(node.nbc) >= o.splitsize
#         nbc, spl, left_nbc, right_nbc = split(node.nbc)
#         # if spl.ig > o.cp
#             node.split = spl
#             node.children = [length(o.tree) + 1, length(o.tree) + 2]
#             t = length(o.tree)
#             left =  NBNode(left_nbc,  id = t + 1, parent = i)
#             right = NBNode(right_nbc, id = t + 2, parent = i)
#             push!(o.tree, left)
#             push!(o.tree, right)
#         # end
#     end
# end

# function whichleaf(o::NBTree, x::VectorOb)
#     i = 1
#     node = o.tree[i]
#     while length(node.children) > 0
#         i = node.children[whichchild(node.split, x)]
#         node = o.tree[i]
#     end
#     i, node
# end

# function classify(o::NBTree, x::VectorOb)
#     i, node = whichleaf(o, x)
#     classify_node(node.nbc)
# end
# function classify(o::NBTree, x::AbstractMatrix)
#     mapslices(xi -> classify(o, xi), x, 2)
# end

ext: jl


text:
using Test

@testset "0.1.0: Lib provides high-level APIs for implementing BLS Signatures scheme perfomantly" begin
    # https://medium.com/cryptoadvance/bls-signatures-better-than-schnorr-5a7fe30ea716
    # e(P, H(m)) = e(G, S)
    m = "test"
    pk = mod(rand(BN), curve_order(EP))
    G = curve_gen(EP)
    e(P::EP, Q::EP2) = field_final_exp(curve_miller(P, Q))
    H(m::String) = curve_map(EP2, Vector{UInt8}(m))
    P = pk * G
    S = pk * H(m)
    @test !isinf(G) && isvalid(G)
    @test !isinf(P) && isvalid(P)
    @test !isinf(S) && isvalid(S)
    @test e(P, H(m)) == e(G, S)
    # print something to the log
    @time e(P, H(m)) == e(G, S)
end


text:
module TestTSXCore

using Test
using TSXPlayground.TSX

const NTRIES = 1000

function first_transacted(f)
    local result
    for _ in 1:NTRIES
        result = f()
        result.transacted && break
    end
    return result
end

function trivial()
    transacted = false
    if (code = xbegin()) == XBEGIN_STARTED
        transacted = true
        xend()
    end
    return (; transacted, code)
end

function simple_abort()
    transacted = false
    if (code = xbegin()) == XBEGIN_STARTED
        transacted = true
        @xabort 42
        xend()
    end
    return (; transacted, code)
end

function check_xtest()
    transacted = false
    outside = xtest()
    inside = 0
    if (code = xbegin()) == XBEGIN_STARTED
        transacted = true
        inside = xtest()
        xend()
    end
    return (; transacted, code, outside, inside)
end

@inline function tryinc!(counter)
    transacted = false
    if (code = xbegin()) == XBEGIN_STARTED
        counter[] += 1
        transacted = true
        xend()
    end
    return (; transacted, code)
end

function check_conflict!(counter, me, other, limit)
    local result
    n = 0
    t = 0
    c = 0
    allcode = Int32(0)
    ok = false
    for _ in 1:limit
        result = tryinc!(counter)
        n += 1
        t += result.transacted
        c += TSX.is_conflict(result.code)
        allcode |= result.code
        if c > 0
            me[] = ok = true
        end
        ok && other[] && break
    end
    return (; n, t, c, allcode, result...)
end

function check_conflict(limit = 1000_0000)
    padded = Vector{Int}(undef, 128)
    counter = view(padded, 64)
    counter[] = 0
    done1 = Threads.Atomic{Bool}(false)
    done2 = Threads.Atomic{Bool}(false)

    local r1
    tid = mod1(Threads.threadid() + 1, Threads.nthreads())
    task = @task begin
        r1 = check_conflict!(counter, done1, done2, limit)
    end
    task.sticky = false
    ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, tid - 1)
    schedule(task)

    r2 = check_conflict!(counter, done2, done1, limit)
    wait(task)

    return (; limit, r1, r2, counter = counter[])
end

function test_conflict()
    result = check_conflict()
    @debug "`check_conflict()`" result
    @test TSX.is_conflict(result.r1.allcode) && TSX.is_conflict(result.r2.allcode)
    @test result.r1.t + result.r2.t == result.counter
    return result
end

@testset "trivial" begin
    result = first_transacted(trivial)
    @test result.transacted == true
end

@testset "simple_abort" begin
    local result
    for _ in 1:NTRIES
        result = simple_abort()
        TSX.is_explicit(result.code) && break
    end
    @test TSX.is_explicit(result.code)
    @test result.transacted == false
    @test xabort_code(result.code) == 42
end

@testset "check_xtest" begin
    result = first_transacted(check_xtest)
    @test result.transacted == true
    @test result.outside == 0
    @test result.inside != 0
end

@testset "check_conflict" begin
    if Threads.nthreads() > 1
        test_conflict()
    end
end

end  # module


text:
# TODO add double dog leg and subspace dogleg
#===============================================================================
  Dogleg is a trust region sub-problem solver used to generate a cheap and crude
  approximation to the solution. If the Cauchy step is outside of the trust re-
  gion it will scale it down to have the length of the radius. If Cauchy step
  is in the interior, it will find the intersection between the Newton step and
  the Cauchy step. Since the trust region is a Euclidean Ball, it is simple to
  find this point on the trust region boundary.

  The Dogleg solver is only appropriate for positive definite Hessians.
===============================================================================#
"""
    Dogleg()

A trust region sub-problem solver that assumes positive definite hessians (exact
or quasi-Newton approximations such as BFGS or variants).
"""
struct Dogleg{T} <: TRSPSolver
    γ::T # unused, for double-dogleg
end
Dogleg() = Dogleg(nothing)

function (dogleg::Dogleg)(∇f, H, Δ, p, scheme; abstol = 1e-10, maxiter = 50)
    T = eltype(p)
    n = length(∇f)

    # find the Cauchy point; assumes ∇f is not ≈ 0
    d_cauchy = -∇f * norm(∇f)^2 / (∇f' * H * ∇f)

    # If it lies outside of the trust region, accept the Cauchy point and
    # move on
    norm_d_cauchy = norm(d_cauchy)
    if norm_d_cauchy ≥ Δ
        shrink = Δ / norm_d_cauchy # inv(Δ/norm_d_cauchy) puts it on the border
        p .= d_cauchy .* shrink
        interior = false
    else
        # Else, calculate (Quasi-)Newton step. If this is interior, then take the
        # step. Otherwise find where the dog-leg path crosses the trust region

        # find the (quasi-)Newton step
        p = find_direction!(p, H, nothing, ∇f, scheme)
        norm_p = norm(p)
        if norm_p ≤ Δ # fixme really need to add the 20% slack here (see TR book and NTR)
            p .= p
            if norm_p < Δ
                interior = true
            else
                interior = false
            end
        else
            # we now solve a quadratic to find the step size t from d_cauchy
            # towards p. We use a numerically stable way of doing this (see
            # any numerical analysis text book) See [NW, p. 75] for the expression
            # giving rise to the quadratic equation
            dot_cachy_p = dot(d_cauchy, p)

            # a is ||d_c - d_n||^2 expanded into scalar operations
            a = norm_d_cauchy^2 + norm_p^2 - 2 * dot_cachy_p
            b = -dot_cachy_p * norm_d_cauchy^2
            c = norm_d_cauchy^2 - Δ^2 # move the rhs over
            q = -(b + sign(b) * √(b^2 - 4 * a * c)) / 2

            # since we know that c is necessarily negative (since d_cauchy was
            # not at the border) the discriminant is positive, and there are two
            # roots - pick the positive one. There has to be one positive and one
            # negative.
            if b > 0
                # if b is positive, q is negative. Then if c is negative, we must
                # have that c / q is the positive root.
                t = c / q
            elseif b < 0
                # else, the other root must be positive
                t = q / a
            else
                t = T(0)
            end
            p .= d_cauchy .+ t .* p
            interior = false
        end
    end
    m = dot(∇f, p) + dot(p, H * p) / 2

    return (
        p = p,
        mz = m,
        interior = interior,
        λ = nothing,
        hard_case = nothing,
        solved = true,
        Δ = Δ,
    )
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
module Kalman

using Distributions, GaussianDistributions
export LinearHomogSystem, kalmanfilter, kalmanfilter!, kalmanrts, kalmanrts!, kalmanEM

# generate
export sample, randmvn

# iterator
export KalmanFilter, MappedKalmanFilter

# track
export track

include("ellipse.jl")

abstract type StateSpaceModel
end

doc"""
```
LinearHomogSystem <: StateSpaceModel

LinearHomogSystem(x0, P0, Phi, b, Q, y, H, R)
```

Time homogeneous linear dynamical model with linear observation scheme corresponding to

$ x[k] = Φx[k−1] + b + w[k],    w[k] ∼ N(0, Q) $

$ y[k] = Hx[k] + v[k],    v[k] ∼ N(0, R) $

The argument `y` to the constructor is a dummy variable with the same type and `size`
as the observations.

# Example

```
x0 = [1., 0.]
P0 = eye(2)

Phi = [0.8 0.2; 0.0 0.8]
b = zeros(2)
Q = [0.1 0.0; 0.0 1.0]

y = [NaN]
H = [1.0 0.0]
R = eye(1)
M = LinearHomogSystem(x0, P0, Phi, b, Q, y, H, R)
```

"""
mutable struct LinearHomogSystem{Tx,TP,Ty,TPhi,T3,T4} <: StateSpaceModel
        # Initial x1 ~ N(x0, P0)
        x0::Tx  # d
        P0::TP # dxd

        # Evolution x -> Phi x + b + w, w ~ N(0, Q)
        Phi::TPhi # dxd
        b::Tx # d
        Q::TP # dxd

        # Observational model y = H*x + v, v ~ N(0, R)
        y::Ty # d₂ "shadow" value
        H::T3 # d₂xd
        R::T4 # d₂xd₂
end
prior(M) = Gaussian(M.x0, M.P0)
dims(SSM) = size(SSM.H)
llikelihood(yres, S, SSM) = logpdf(Gaussian(zero(yres), S), yres)

doc"""
    predict!(s, x, P, t, M) -> x, Ppred, Phi

Computes predicted state `x` and predicted covariance of model `M` according
the prediction equations

$P_{pred} = \Phi P \Phi' + Q$

Also returns evolution operator / Jacobian `Phi`.
"""
function predict!(s, x, P, t, M::LinearHomogSystem)
    x = evolve(s, x, P, t, M)
    Ppred = M.Phi*P*M.Phi' + M.Q
    x, Ppred, M.Phi
end

doc"""
    evolve(s, x, P, t, M) -> x

Evolve state `x` according to the state dynamics of model `M`
without noise, usually

$ x \mapsto \Phi x + b. $
"""
function evolve(s, x, P, t, M::LinearHomogSystem)
    M.Phi*x + M.b
end

"""
    observe!(s, x, P, t, y, M) -> t, y, H, R

Obtain observation time `t`, observation `y`, observation matrix `H`
and observation covariance `H` using model `M`. Except for certain use cases this is
a no-op on arguments `t` and `y` (but needs to return `t` and `y` nevertheless.)
"""
function observe!(s, x, P, t, Y, M::LinearHomogSystem)
    t, Y, M.H, M.R
end

"""
    correct!(x, Ppred, y, H, R, M) -> x, P, yres, S, K

Perform the Kalman correction step using state space model `M`.
Inputs are predicted state `x` and predicted covariance, as well as observation `y`
and observation matrix `H`, observation covariance `R`.

Outputs filtered state `x`, filtered state covariance `P`, residual `y`
as well as innovation covariance `S` and Kalman gain `K`.
"""
function correct!(x, Ppred, y, H, R, SSM)
    yres = y - H*x # innovation residual

    S = H*Ppred*H' + R # innovation covariance

    K = Ppred*H'/S # Kalman gain
    x = x + K*yres
    P = (I - K*H)*Ppred # Ppred - Ppred*H' * inv(H*Ppred*H' + R) * H* Ppred
    x, P, yres, S, K
end

"""
    kalman_kernel(s, x, P, t, Y, SSM) -> t, x, P, Ppred, ll, K

Single Kalman filter step consisting of a prediction step `predict!`, an observation step `observe!`
and a correction step `correct!`. Return filtered covariance `P` and predicted `Ppred`

Computes and returns as well the log likelihood of the residual and the Kalman gain.
"""
function kalman_kernel(s, x, P, t, Y, SSM)

    x, Ppred, Phi = predict!(s, x, P, t, SSM)

    t, y, H, R = observe!(s, x, P, t, Y, SSM)

    x, P, yres, S, K = correct!(x, Ppred, y, H, R, SSM)

    ll = llikelihood(yres, S, SSM)

    t, x, P, Ppred, ll, K
end

"""
    ll, xxf = kalmanfilter(yy::Matrix, M::LinearHomogSystem)

Kalman filter

    yy -- d₂xn array
    M -- linear model
    ll -- marginal likelihood
    xxf -- filtered process
"""
function kalmanfilter(yy::AbstractMatrix, M::LinearHomogSystem)
    d = size(M.Phi, 1)
    n = size(yy, 2)
    kalmanfilter!(0:n, yy, zeros(d, n), zeros(d, d, n), zeros(d, d, n), M)
end

"""
    ll, X = kalmanfilter{T}(Y::Array{T,3}, M::LinearHomogSystem)

Stack Kalman filter

    Y -- array of m independent processes (dxnxm array)

    X -- filtered processes (dxnxm array)
    ll -- product marginal log likelihood
"""
function kalmanfilter(Y::AbstractArray{T,3}, M::LinearHomogSystem) where {T}
    d₂, n, m= size(Y)
    d = length(M.x0)
    X = zeros(d,n,m)
    PP = zeros(d, d, n)
    PPpred = zeros(d, d, n)
    ll = 0.
    for j in 1:m
        ll += kalmanfilter!(0:n, view(Y, :, :, j), view(X, :, :, j),  PP, PPpred, M)[4]
    end
    ll, X
end

# Matrix-of-vectors-version
function kalmanfilter(Y::AbstractMatrix{Ty}, M::LinearHomogSystem{Tx,TP,Ty}) where {Tx,TP,Ty}
    n, m = size(Y)
    X = zeros(Tx, n, m)
    PP = zeros(TP, n)
    PPpred = zeros(TP, n)
    ll = 0.
    for j in 1:m
        ll += kalmanfilter!(0:n, view(Y, :, j), view(X, :, j),  PP, PPpred, M)[4]
    end
    ll, X
end

function kalmanfilter!(tt, yy, xxf, PP, PPpred, M::LinearHomogSystem{Vector{T}}) where {T}

    d₂, d = dims(M)
    assert(ndims(yy) == 2)
    n = size(yy,2)

    xf = M.x0
    P = M.P0

    ll = 0.0
    t = tt[1]

    for i in 1:n # t, x, P, Ppred, ll, K = kalman_kernel(s, x, P, t, Y, SSM)
        t, xf, P, Ppred, l, _ = kalman_kernel(t, xf, P, tt[i+1], yy[.., i], M)
        xxf[.., i], PP[.., i], PPpred[.., i] = xf, P, Ppred
        ll += l
    end
    xxf, PP, PPpred, ll
end

function kalmanfilter!(tt, yy, xxf, PP, PPpred, M::LinearHomogSystem)

    n = size(yy, 1)

    xf = M.x0
    P = M.P0

    ll = 0.0
    t = tt[1]

    for i in 1:n # t, x, P, Ppred, ll, K = kalman_kernel(s, x, P, t, Y, SSM)
        t, xf, P, Ppred, l, _ = kalman_kernel(t, xf, P, tt[i+1], yy[i], M)
        xxf[i], PP[i], PPpred[i] = xf, P, Ppred
        ll += l
    end
    xxf, PP, PPpred, ll
end

function smoother_kernel(xs, Ps, xf, Pf, Ppred, Phi, b)
    # xs -- previous #h[i+1]
    # Ps -- previous #H[i+1]
    # xf -- xxf[:, i] #m[i]
    # Pf -- PPf[:, :, i] # Cn
    # Ppred -- PPpred[:, :, i+1] # R[i+1] = C[i] + W[i+1]

    J = Pf*Phi'/Ppred # C/(C+w)
    xs = xf +  J*(xs - (Phi*xf  + b))
    Ps = Pf + J*(Ps - Ppred)*J'

    xs, Ps, J
end

function backward_kernel(xs, Ps, xf, Pf, Ppred, Phi, b)
    # xs -- previous
    # Ps -- previous
    # xf -- xxf[:, i]
    # Pf -- PPf[:, :, i]
    # Ppred -- PPpred[:, :, i+1]

    J = Pf*Phi'/Ppred
    xs = xf +  J*(xs - (Phi*xf  + b))
    Ps = Pf + J*Ppred*J' # as Ppred = M.Phi*P*M.Phi' + M.Q

    xs, Ps, J
end

"""
    xxs, PP, PPpred, ll = kalmanrts(yy, xxs, M::LinearHomogSystem)

Rauch-Tung-Striebel smoother

    yy -- d₂xn array or
    xxs -- empty dxn array or view
    M --  linear dynamic system and observation model

    xxs -- smoothed process
    PP -- smoother variance
    PPpred -- smoother prediction
    ll -- filter likelihood
"""
function kalmanrts(yy::Matrix, M::LinearHomogSystem{Vector{T}}) where {T}
    d = size(M.Phi, 1)
    n = size(yy, 2)
    kalmanrts!(yy, zeros(d, n), zeros(d, d, n), zeros(d, d, n), M)
end

function kalmanrts!(yy, xx, PP, PPpred, M::LinearHomogSystem{Vector{T}}) where {T}  # using Rauch-Tung-Striebel
    d₂, d = size(M.H)
    assert(ndims(yy) == 2)
    n = size(yy,2)

    # forward pass

    # in place! xx = xxf = xxs
    # PP = PPf = PPs

    xx, PP, PPpred, ll = kalmanfilter!(0:n, yy, xx, PP, PPpred, M)

    # backwards pass
    Phi, b = M.Phi, M.b

    # start with xf, P from forward pass
    xs = xx[:, n]
    Ps = PP[:, :, n]

    for i in n-1:-1:1
        xs, Ps, J = smoother_kernel(xs, Ps, xx[:, i],  PP[:, :, i], PPpred[:, :, i+1], Phi, b)
        xx[:, i] = xs
        PP[:, :, i] = Ps
    end

    xx, PP, PPpred, ll
end

"""
    X = kalmanrts{T}(Y::Array{T,3}, M::LinearHomogSystem)

Stack Rauch-Tung-Striebel smoother, computes the marginal smoothed states, that is it computes the
law ``p(x_i | y_{1:n})``.

    Y -- m independent processes (d₂xnxm array)

    X -- smoothed processes (dxnxm array)
"""
function kalmanrts(Y::Array{T,3}, M::LinearHomogSystem{Vector{T}}) where {T}
    d₂, n, m= size(Y)
    d = length(M.x0)
    X = zeros(d,n,m)
    PP = zeros(d, d, n)
    PPpred = zeros(d, d, n)
    ll = 0.
    for j in 1:m
        ll += kalmanrts!(view(Y, :, :, j), view(X, :, :, j), PP, PPpred, M)[4]
    end
    ll, X
end

include("generate.jl")

include("iterator.jl")

include("track.jl")

#include("kalmanem.jl")

end # module


text:
using TypeDBClient
using Test

@testset "TypeDBClient" begin
include("unit/test_common.jl")
include("unit/test_request_builder.jl")
include("unit/test_concepts.jl")

end


text:
using FractionalTimeDG
using Test
using OffsetArrays
import FractionalTimeDG, SpecialFunctions
using QuadGK: quadgk

Γ = SpecialFunctions.gamma

include("common.jl")

@testset "FractionalTimeDG.jl" begin
    include("test_coef_G_K.jl")
    include("test_legendre.jl")
    include("test_ABC.jl")
    include("test_H_uniform.jl")
    include("test_H.jl")
    include("test_chebyshev.jl")
end


text:
function solve_heat(problem::FDHeatProblem{TBCL, TBCR, TS}) where {TBCL<:PeriodicBC, TBCR<:PeriodicBC, TS<:BEuler}
    u = similar(problem.u0);
    @. u = problem.u0;
    assemble_system!(problem);
    t = 0.;

    u_trajectory = [deepcopy(u)];
    t_trajectory = [t];

    for j in 1:problem.nΔt
        @. problem.rhs = u;
        u .= problem.A\problem.rhs;
        t+=problem.Δt;
        push!(u_trajectory, deepcopy(u));
        push!(t_trajectory, t);
    end

    return t_trajectory, u_trajectory
end


text:
baremodule Tuples
using Base

# Since we're defining our own getindex and length methods to operate on types
# separately, it's difficult to then use Base's methods at the same time.
# Instead, we define them as t* within an Implementation module, and reassign
# them to the API names we want to use in this baremodule.
Base.include("implementation.jl")
const getindex = Implementation.tgetindex
const length = Implementation.tlength
const collect = Implementation.tcollect
const concatenate = Implementation.concatenate

end


text:
const FIXUP_PR_TITLE = "TagBot: Use issue comment triggers"
const FIXUP_COMMIT_MESSAGE = "[ci skip] $FIXUP_PR_TITLE"
const FIXUP_PR_BODY = """
As requested, I've updated your TagBot configuration to use issue comment triggers.

Please note that this PR does not take into account your existing configuration.
If you had any custom configuration, you'll need to add it back yourself.
"""

const TAGBOT_YML = raw"""
name: TagBot
on:
  issue_comment:
    types:
      - created
  workflow_dispatch:
jobs:
  TagBot:
    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
    runs-on: ubuntu-latest
    steps:
      - uses: JuliaRegistries/TagBot@v1
        with:
          token: ${{ secrets.GITHUB_TOKEN }}
          ssh: ${{ secrets.DOCUMENTER_KEY }}
"""

function should_fixup(repo, issue)
    return fixup_comment_exists(repo, issue) &&
           !fixup_done(repo) &&
           tagbot_file(repo; issue_comments=true) === nothing
end

function get_fork(repo)
    fork = GH.create_fork(repo; auth=AUTH[])
    # Make sure the fork is new, otherwise it might be outdated.
    return if now(UTC) - fork.created_at > Minute(1)
        GH.delete_repo(fork; auth=AUTH[])
        get_fork(repo)
    else
        fork
    end
end

function open_fixup_pr(repo; branch="tagbot/$(randstring())", fork=get_fork(repo))
    head = GH.commit(fork, "HEAD"; auth=AUTH[])
    GH.create_reference(
        fork; auth=AUTH[], params=(; sha=head.sha, ref="refs/heads/$branch")
    )
    path, contents = tagbot_file(fork)
    GH.update_file(
        fork,
        path;
        auth=AUTH[],
        params=(;
            branch=branch,
            content=base64encode(TAGBOT_YML),
            message=FIXUP_COMMIT_MESSAGE,
            sha=bytes2hex(sha1("blob $(length(contents))\0$contents")),
        ),
    )
    result = try
        GH.create_pull_request(
            repo;
            auth=AUTH[],
            params=(;
                title=FIXUP_PR_TITLE,
                body=FIXUP_PR_BODY,
                head="$(TAGBOT_USER[]):$branch",
                base=fork.default_branch,
            ),
        )
    catch ex
        @error "Encountered an error while trying to create the fixup PR" exception = (
            ex, catch_backtrace()
        )
        nothing
    end
    return result
end

function fixup_comment_exists(repo, issue)
    kwargs = Dict(:auth => AUTH[], :page_limit => 1, :params => (; per_page=100))
    while true
        comments, pages = GH.comments(repo, issue; kwargs...)
        for comment in comments
            if is_fixup_trigger(comment)
                return true
            end
        end
        if haskey(pages, "next")
            delete!(kwargs, :params)
            kwargs[:start_page] = pages["next"]
        else
            return false
        end
    end
end

function is_fixup_trigger(comment)
    return comment.user.login != TAGBOT_USER[] && occursin(r"TagBot fix"i, comment.body)
end

function fixup_done(repo)
    pulls, _ = GH.pull_requests(
        repo; auth=AUTH[], params=(; creator=TAGBOT_USER[], state="all")
    )
    for pull in pulls
        if pull.title == FIXUP_PR_TITLE
            return true
        end
    end
    return false
end


text:
const GeometryLike = Union{AbstractGeometry, AbstractMesh, MeshFileGeometry}

"""
Represents a three.js Object, consisting of a geometry and a material.
"""
struct Object{G <: GeometryLike, M <: AbstractMaterial} <: AbstractObject
    geometry::G
    material::M
    _type::String
end

geometry(o::Object) = o.geometry
material(o::Object) = o.material
threejs_type(o::Object) = o._type

defaultmaterial(args...; kw...) = MeshLambertMaterial(args...; kw...)

# Default object types for geometries, point clouds, and triads
Object(g::GeometryLike) = MeshObject(g)
Object(g::GeometryLike, m::AbstractMaterial) = MeshObject(g, m)
Object(c::PointCloud) = Points(c)
Object(c::PointCloud, m::AbstractMaterial) = Points(c, m)
Object(t::Triad) = LineSegments(t, LineBasicMaterial(vertexColors=2))

@deprecate Mesh(args...) MeshObject(args...)

MeshObject(g, m) = Object(g, m, "Mesh")
MeshObject(geometry::GeometryLike) = MeshObject(geometry, defaultmaterial())

Points(g, m) = Object(g, m, "Points")
Points(geometry::GeometryLike; kw...) = Points(geometry, PointsMaterial(kw...))

for line_type in [:LineSegments, :Line, :LineLoop]
    @eval $line_type(g::AbstractGeometry, m::AbstractMaterial=LineBasicMaterial()) = Object(g, m, $(string(line_type)))
    @eval $line_type(points::AbstractVector{<:Point}, m::AbstractMaterial=LineBasicMaterial()) = $line_type(PointCloud(points), m)
end

struct PngImage
    data::Vector{UInt8}
end

PngImage(fname::AbstractString) = PngImage(open(read, fname))

@with_kw struct Texture
    image::PngImage
    wrap::Tuple{Int, Int} = (1001, 1001)  # TODO: replace with enum
    repeat::Tuple{Int, Int} = (1, 1)      # TODO: what does this mean?
end

@with_kw mutable struct GenericMaterial <: AbstractMaterial
    _type::String
    color::RGBA{Float32} = RGB(1., 1., 1.)
    map::Union{Texture, Nothing} = nothing
    depthFunc::Int = 3
    depthTest::Bool = true
    depthWrite::Bool = true
    linewidth::Float64 = 1.
    vertexColors::Int = 0    # TODO: make an enum
    side::Int = 2            # TODO: make an enum https://github.com/mrdoob/three.js/blob/d55897b8e9b2632896d8ac146a05b3b4be3668f8/src/constants.js#L14
    wireframe::Bool = false
    wireframeLinewidth::Float64 = 1
end

threejs_type(m::GenericMaterial) = m._type

MeshBasicMaterial(;kw...) = GenericMaterial(_type="MeshBasicMaterial"; kw...)
MeshLambertMaterial(;kw...) = GenericMaterial(_type="MeshLambertMaterial"; kw...)
MeshPhongMaterial(;kw...) = GenericMaterial(_type="MeshPhongMaterial"; kw...)
LineBasicMaterial(;kw...) = GenericMaterial(_type="LineBasicMaterial"; kw...)

@with_kw struct PointsMaterial <: AbstractMaterial
    color::RGBA{Float32}=RGB(1., 1., 1.)
    size::Float32 = 0.002
    vertexColors::Int = 2
end


text:
function iterate_instructions(f, session)
    regs = Gallium.getregs(session)
    addr = UInt64(Gallium.ip(regs))
    ctx = DisAsmContext()
    cont = true
    while cont
        insts = Gallium.load(session, RemotePtr{UInt8}(addr), 15)
        (Inst, InstSize) = getInstruction(insts, 0; ctx = ctx)
        cont = f(addr, Inst, ctx)
        addr += InstSize
        free(Inst)
        Inst = nothing
    end
end

function NextBranchAddr(session)
    addr = 0
    iterate_instructions(session) do where, Inst, ctx
        mayBranch(Inst, ctx) && return false
        addr = where
        return true
    end
    return addr
end

function LastInstrInRange(session, range)
    addr = 0
    iterate_instructions(session) do where, Inst, ctx
        where > last(range) && return false
        addr = where
        return true
    end
    return addr
end

task_step_until_bkpt!(session) = Gallium.step_until_bkpt!(session)

function step_to_address(session, addr)
    bp = Gallium.breakpoint(session, addr)
    task_step_until_bkpt!(session)
    Gallium.disable(bp)
    nothing
end

function step_over(session, range)
    while true
        nba = NextBranchAddr(session)
        isbranch = true
        if !(nba ∈ range)
            nba = LastInstrInRange(session, range)
            isbranch = false
        end
        bp = Gallium.breakpoint(session, nba)
        # get_frame
        while true
            task_step_until_bkpt!(session)
            #comp = compare_frames(frame, timeline)
            #if older
            #    Gallium.disable(bp)
            #    return
            #elseif younger
            #    continue
            #else
            #    break
            #end
            break
        end
        # Ok, we've arrived at our breakpoint in the same frame
        Gallium.disable(bp)
        # We're at an instruction that may step out of the range. Single
        # step and see where we are
        Gallium.single_step!(session)
        (UInt64(Gallium.ip(session)) ∈ range) && continue
        break
    end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/mixed_convexity.jl
# Defines the operators: cos, sin, sinh, atanh, tanh, atan, asin, tan, acos,
#                        asinh, cosh, rad2deg, deg2rad, sind, cosd..., sec,
#                        csc..., sech, csch, ...asech...
#############################################################################

# convex relaxation (envelope) of cos function
@inline function cv_cos(x::Float64, xL::Float64, xU::Float64, tp1::Float64, tp2::Float64)
     r = 0.0
     kL = Base.ceil(-0.5 - xL/(2.0*pi))
     if x <= (pi - 2.0*pi*kL)
         xL1 = xL + 2.0*pi*kL
         if xL1 >= pi/2.0
             return cos(x), -sin(x), tp1, tp2
         end
         xU1 = min(xU+2.0*pi*kL,pi)
         if (xL1 >= -pi/2) && (xU1 <= pi/2)
             r = (abs(xL - xU) < MC_ENV_TOL) ? 0.0 : (cos(xU) - cos(xL))/(xU - xL)
             return cos(xL) + r*(x - xL), r, tp1, tp2
         end
         val, dval, tp1 = cv_cosin(x + (2.0*pi)*kL, xL1, xU1, tp1)
         return val, dval, tp1, tp2
     end
     kU = Base.floor(0.5 - xU/(2.0*pi))
     if (x >= -pi - 2.0*pi*kU)
         xU2 = xU + 2.0*pi*kU
         (xU2 <= -pi/2.0) && (return cos(x), -sin(x), tp1, tp2)
         val, dval, tp2 = cv_cosin(x + 2.0*pi*kU, max(xL + 2.0*pi*kU, -pi), xU2, tp2)
         return val, dval, tp1, tp2
     end
     return -1.0, 0.0, tp1, tp2
end
# function for computing convex relaxation over nonconvex and nonconcave regions
@inline function cv_cosin(x::Float64, xL::Float64, xU::Float64, xj::Float64)
    if abs(xL) <= abs(xU)
        left = false
        x0 = xU
        xm = xL
    else
        left = true
        x0 = xL
        xm = xU
    end
    if xj === Inf
        xj, flag = newton(x0, xL, xU, cv_cosenv, dcv_cosenv, xm, 0.0)
        flag && (xj = golden_section(xL, xU, cv_cosenv, xm, 0.0))
    end
    if (left && x <= xj) || (~left && x >= xj)
        return cos(x), -sin(x), xj
    else
        r = abs(xm - xj) < MC_ENV_TOL ? 0.0 : (cos(xm) - cos(xj))/(xm - xj)
    end
    return cos(xm) + r*(x - xm), r, xj
end
@inline cv_cosenv(x::Float64, y::Float64, z::Float64) = (x - y)*sin(x) + cos(x) - cos(y)
@inline dcv_cosenv(x::Float64, y::Float64, z::Float64) = (x - y)*cos(x)
# concave relaxation (envelope) of cos function
@inline function cc_cos(x::Float64, xL::Float64, xU::Float64, tp1::Float64, tp2::Float64)
    temp = cv_cos(x - pi, xL - pi, xU - pi, tp1, tp2)
    return -temp[1], -temp[2], temp[3], temp[4]
end
@inline function cos_arg(xL::Float64, xU::Float64)
    kL = Base.ceil(-0.5 - xL/(2.0*pi))
    xL1 = xL + 2.0*pi*kL
    xU1 = xU + 2.0*pi*kL
    ~((xL1 >= -pi) && (xL1 <= pi)) && (return NaN, NaN)
    if xL1 <= 0.0
        if xU1 <= 0.0
            arg1 = xL
            arg2 = xU
        elseif xU1 >= pi
            arg1 = pi-2.0*pi*kL
            arg2 = -2.0*pi*kL
        else
            arg1 = (cos(xL1) <= cos(xU1)) ? xL : xU
            arg2 = -2.0*pi*kL
        end
        return arg1, arg2
    end
    if xU1 <= pi
        arg1 = xU
        arg2 = xL
    elseif xU1 >= 2.0*pi
        arg1 = pi - 2.0*pi*kL
        arg2 = 2.0*pi - 2.0*pi*kL
    else
        arg1 = pi - 2.0*pi*kL
        arg2 = (cos(xL1) >= cos(xU1)) ? xL : xU
    end
    return arg1, arg2
end
@inline function cos_kernel(x::MC{N, Diff}, y::Interval{Float64}, cv_tp1::Float64,
                              cv_tp2::Float64, cc_tp1::Float64, cc_tp2::Float64) where N
    xL = x.Intv.lo
    xU = x.Intv.hi
    xLc = y.lo
    xUc = y.hi
    eps_min, eps_max = cos_arg(x.Intv.lo, x.Intv.hi)
    midcc = mid3v(x.cv, x.cc, eps_max)
    midcv = mid3v(x.cv, x.cc, eps_min)
    cc, dcc, cc_tp1, cc_tp2 = cc_cos(midcc, x.Intv.lo, x.Intv.hi, cc_tp1, cc_tp2)
    cv, dcv, cv_tp1, cv_tp2 = cv_cos(midcv, x.Intv.lo, x.Intv.hi, cv_tp1, cv_tp2)
    gcc1, gdcc1, cc_tp1, cc_tp2 = cc_cos(x.cv, x.Intv.lo, x.Intv.hi, cc_tp1, cc_tp2)
    gcv1, gdcv1, cv_tp1, cv_tp2 = cv_cos(x.cv, x.Intv.lo, x.Intv.hi, cv_tp1, cv_tp2)
    gcc2, gdcc2, cc_tp1, cc_tp2 = cc_cos(x.cc, x.Intv.lo, x.Intv.hi, cc_tp1, cc_tp2)
    gcv2, gdcv2, cv_tp1, cv_tp2 = cv_cos(x.cc, x.Intv.lo, x.Intv.hi, cv_tp1, cv_tp2)
    cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
    cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
    return MC{N, Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst), cv_tp1, cv_tp2, cc_tp1, cc_tp2
end
@inline function cos_kernel(x::MC{N, T}, y::Interval{Float64}, cv_tp1::Float64,
                              cv_tp2::Float64, cc_tp1::Float64, cc_tp2::Float64) where {N,T<:Union{NS,MV}}
    xL = x.Intv.lo
    xU = x.Intv.hi
    xLc = y.lo
    xUc = y.hi
    eps_min, eps_max = cos_arg(x.Intv.lo, x.Intv.hi)
    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
    cc, dcc, cc_tp1, cc_tp2 = cc_cos(midcc, x.Intv.lo, x.Intv.hi, cc_tp1, cc_tp2)
    cv, dcv, cv_tp1, cv_tp2 = cv_cos(midcv, x.Intv.lo, x.Intv.hi, cv_tp1, cv_tp2)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
    cv,cc,cv_grad,cc_grad = cut(xLc,xUc,cv,cc,cv_grad,cc_grad)
    return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst), cv_tp1, cv_tp2, cc_tp1, cc_tp2
end
@inline function cos(x::MC)
    y, tp1, tp2, tp3, tp4 = cos_kernel(x, cos(x.Intv), Inf, Inf, Inf, Inf)
    return y
end
@inline function sin_kernel(x::MC, y::Interval{Float64}, cv_tp1::Float64,
                            cv_tp2::Float64, cc_tp1::Float64, cc_tp2::Float64)
    cos_kernel(x-pi/2.0, y, cv_tp1, cv_tp2, cc_tp1, cc_tp2)
end
@inline function sin(x::MC)
    y, tp1, tp2, tp3, tp4 = sin_kernel(x, sin(x.Intv), Inf, Inf, Inf, Inf)
    return y
end

@inline sinh_deriv(x::Float64) = cosh(x)
@inline sinh_env(x::Float64, y::Float64, z::Float64) = (x-y)*cosh(x)-(sinh(x)-sinh(y))
@inline function cv_sinh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return sinh(x), sinh_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(sinh, sinh_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, sinh_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, sinh_env, xU, 0.0))
    end
    (x <= p) && (return dline_seg(sinh, sinh_deriv, x, xL, p)..., p)
    return sinh(x), sinh_deriv(x), p
end
@inline function cc_sinh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(sinh, sinh_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return sinh(x), sinh_deriv(x), p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, sinh_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, sinh_env, xL, 0.0))
    end
    (x <= p) && (return sinh(x), sinh_deriv(x), p)
    return dline_seg(sinh, sinh_deriv, x, p, xU)..., p
end

@inline atanh_deriv(x::Float64) = 1.0/(1.0 - x^2)
@inline atanh_env(x::Float64, y::Float64, z::Float64) = (NaNMath.atanh(x) - NaNMath.atanh(y))*(1.0 - x^2) - x + y
@inline function cv_atanh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL <= -1.0) && (return NaN, NaN, NaN)
    (xU >= 1.0) && (return NaN, NaN, NaN)
    (xL >= 0.0) && (return NaNMath.atanh(x), atanh_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(NaNMath.atanh, atanh_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, atanh_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, atanh_env, xU, 0.0))
    end
    (x <= p) && (return dline_seg(NaNMath.atanh, atanh_deriv, x, xL, p)..., p)
    return NaNMath.atanh(x), atanh_deriv(x), p
end
@inline function cc_atanh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL <= -1.0) && (return NaN, NaN, NaN)
    (xU >= 1.0) && (return NaN, NaN, NaN)
    (xL >= 0.0) && (return dline_seg(NaNMath.atanh, atanh_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return NaNMath.atanh(x), atanh_deriv(x), p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, atanh_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, atanh_env, xL, 0.0))
    end
    (x <= p) && (return NaNMath.atanh(x), atanh_deriv(x), p)
    return dline_seg(NaNMath.atanh, atanh_deriv, x, p, xU)..., p
end

@inline tanh_deriv(x::Float64) = sech(x)^2
@inline tanh_env(x::Float64, y::Float64, z::Float64) = (x - y) - (tanh(x) - tanh(y))/(1.0 - tanh(x)^2)
@inline tanh_envd(x::Float64, y::Float64, z::Float64)= 2.0*tanh(x)/(1.0 - tanh(x)^2)*(tanh(x) - tanh(y))
@inline function cv_tanh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(tanh, tanh_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return tanh(x), sech(x)^2, p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, tanh_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, tanh_env, xU, 0.0))
    end
    (x <= p) && (return tanh(x), sech(x)^2, p)
    return dline_seg(tanh, tanh_deriv, x, p, xU)..., p
end
@inline function cc_tanh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return tanh(x), sech(x)^2, p)
    (xU <= 0.0) && (return dline_seg(tanh, tanh_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, tanh_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, tanh_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(tanh, tanh_deriv, x, xL, p)..., p)
    return tanh(x), sech(x)^2, p
end

@inline atan_deriv(x::Float64) = 1.0/(1.0 + x^2)
@inline atan_env(x::Float64, y::Float64, z::Float64) = (x - y) - (1.0 + sqr(x))*(atan(x) - atan(y))
@inline function cv_atan(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(atan, atan_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return atan(x), 1.0/(1.0+x^2), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, atan_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, atan_env, xU, 0.0))
    end
    (x <= p) && (return atan(x), 1.0/(1.0+x^2), p)
    return dline_seg(atan, atan_deriv, x, p, xU)..., p
end
@inline function cc_atan(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return atan(x), 1.0/(1.0+x^2), p)
    (xU <= 0.0) && (return dline_seg(atan, atan_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, atan_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, atan_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(atan, atan_deriv, x, xL, p)..., p)
    return atan(x), 1.0/(1.0+x^2), p
end

@inline asin_deriv(x::Float64) = 1.0/NaNMath.sqrt(1.0 - x^2)
@inline function asin_env(x::Float64, y::Float64, z::Float64)
    return (NaNMath.asin(x) - NaNMath.asin(y))*NaNMath.sqrt(1.0-x^2) - x + y
end
@inline function cv_asin(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return NaNMath.asin(x), 1.0/NaNMath.sqrt(1.0-x^2), p)
    (xU <= 0.0) && (return dline_seg(NaNMath.asin, asin_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, asin_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, asin_env, xL, 0.0))
  end
  (x <= p) && (return dline_seg(NaNMath.asin, asin_deriv, x, xL, p)..., p)
  return NaNMath.asin(x), 1.0/NaNMath.sqrt(1.0-x^2, p)
end
@inline function cc_asin(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(NaNMath.asin, asin_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return NaNMath.asin(x), 1.0/NaNMath.sqrt(1.0-x^2), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, asin_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, asin_env, xU, 0.0))
    end
    (x <= p) && (return NaNMath.asin(x), 1.0/NaNMath.sqrt(1.0-x^2), p)
    return dline_seg(NaNMath.asin, asin_deriv, x, p, xU)..., p
end

@inline tan_deriv(x::Float64) = sec(x)^2
@inline function tan_env(x::Float64, y::Float64, z::Float64)
    return (x - y) - (NaNMath.tan(x) - NaNMath.tan(y))/(1.0 + NaNMath.tan(x)^2)
end
@inline function tan_envd(x::Float64, y::Float64, z::Float64)
    return 2.0*NaNMath.tan(x)/(1.0 + NaNMath.tan(x)^2)*(NaNMath.tan(x) - NaNMath.tan(y))
end
@inline function cv_tan(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return NaNMath.tan(x), sec(x)^2, p)
    (xU <= 0.0) && (return dline_seg(NaNMath.tan, tan_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, tan_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, tan_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(NaNMath.tan, tan_deriv, x, xL, p)..., p)
    return NaNMath.tan(x), sec(x)^2, p
end
@inline function cc_tan(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(NaNMath.tan, tan_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return NaNMath.tan(x), sec(x)^2, p)
    if p === Inf
        p, flag = secant(0.0, xL, xL, 0.0, tan_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, tan_env, xU, 0.0))
    end
    (x <= p) && (return NaNMath.tan(x), sec(x)^2, p)
    return dline_seg(NaNMath.tan, tan_deriv, x, p, xU)..., p
end

@inline acos_deriv(x::Float64) = -1.0/NaNMath.sqrt(1.0-x^2)
@inline acos_env(x::Float64, y::Float64, z::Float64) = -(NaNMath.acos(x) - NaNMath.acos(y))*NaNMath.sqrt(1-x^2) - x + y
@inline function cc_acos(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(NaNMath.acos, acos_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return NaNMath.acos(x), -1.0/NaNMath.sqrt(1.0-x^2), p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, acos_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, acos_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(NaNMath.acos, acos_deriv, x, xL, p)..., p)
    return NaNMath.acos(x), -1.0/NaNMath.sqrt(1.0-x^2), p
end
@inline function cv_acos(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return NaNMath.acos(x), -1.0/NaNMath.sqrt(1.0-x^2), p)
    (xU <= 0.0) && (return dline_seg(NaNMath.acos, acos_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, acos_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, acos_env, xU, 0.0))
    end
    (x <= p) && (return NaNMath.acos(x), -1.0/NaNMath.sqrt(1.0-x^2), p)
    return dline_seg(acos, acos_deriv, x, p, xU)..., p
end

@inline asinh_deriv(x::Float64) = 1.0/sqrt(1.0 + x^2)
@inline asinh_env(x::Float64, y::Float64, z::Float64) = (x - y) - sqrt(1.0 + sqr(x))*(asinh(x) - asinh(y))
@inline function cv_asinh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(asinh, asinh_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return asinh(x), 1.0/sqrt(1.0 + x^2), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, asinh_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, asinh_env, xU, 0.0))
    end
    (x <= p) && (return asinh(x), 1.0/sqrt(1.0+x^2), p)
    return dline_seg(asinh, asinh_deriv, x, p, xU)..., p
end
@inline function cc_asinh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return asinh(x), 1.0/sqrt(1.0 + x^2), p)
    (xU <= 0.0) && (return dline_seg(asinh, asinh_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, asinh_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, asinh_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(asinh, asinh_deriv, x, xL, p)..., p)
    return asinh(x), 1.0/sqrt(1.0+x^2), p
end

@inline erf_deriv(x::Float64) = (2.0/sqrt(pi))*exp(-x^2)
@inline erf_deriv2(x::Float64) = 4.0*x*exp(-x^2)/sqrt(pi)
@inline function erf_env(x::Float64, y::Float64, z::Float64)
    (x - y) - (erf(x) - erf(y))/erf_deriv(x)
end
@inline function cv_erf(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(erf, erf_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return erf(x), erf_deriv(x), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, erf_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, erf_env, xU, 0.0))
    end
    (x <= p) && (return erf(x), erf_deriv(x), p)
    return dline_seg(erf, erf_deriv, x, p, xU)..., p
end
@inline function cc_erf(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return erf(x), erf_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(erf, erf_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, erf_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, erf_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(erf, erf_deriv, x, xL, p)..., p)
    return erf(x), erf_deriv(x), p
end

@inline erfinv_deriv(x::Float64) = (sqrt(pi)/2.0)*exp(erfinv(x)^2)
@inline erfinv_deriv2(x::Float64) = (sqrt(pi)/2.0)*exp(2.0*erfinv(x)^2)*erfinv(x)
@inline function erfinv_env(x::Float64, y::Float64, z::Float64)
    (x - y) - (erfinv(x) - erfinv(y))/erfinv_deriv(x)
end
@inline function cv_erfinv(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return erfinv(x), erfinv_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(erfinv, erfinv_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, erfinv_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, erfinv_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(erfinv, erfinv_deriv, x, xL, p)..., p)
    return erfinv(x), erfinv_deriv(x), p
end
@inline function cc_erfinv(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(erfinv, erfinv_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return erfinv(x), erf_deriv(x), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, erfinv_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, erfinv_env, xU, 0.0))
    end
    (x >= p) && (return dline_seg(erfinv, erfinv_deriv, x, p, xU)..., p)
    return erfinv(x), erfinv_deriv(x), p
end

@inline erfc_deriv(x::Float64) = (-2.0/sqrt(pi))*exp(-x^2)
@inline erfc_deriv2(x::Float64) = (4.0/sqrt(pi))*exp(-x^2)*x
@inline function erfc_env(x::Float64, y::Float64, z::Float64)
    (x - y) - (erfc(x) - erfc(y))/erfc_deriv(x)
end
@inline function cv_erfc(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return erfc(x), erfc_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(erfc, erfc_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, erfc_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, erfc_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(erfc, erfc_deriv, x, xL, p)..., p)
    return erfc(x), erfc_deriv(x), p
end
@inline function cc_erfc(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(erfc, erfc_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return erfc(x), erfc_deriv(x), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, erfc_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, erfc_env, xU, 0.0))
    end
    (x >= p) && (return dline_seg(erfc, erfc_deriv, x, p, xU)..., p)
    return erfc(x), erfc_deriv(x), p
end

@inline cbrt_deriv(x::Float64) = 1.0/(3.0*cbrt(x)^2)
@inline cbrt_deriv2(x::Float64) = -2.0/(9.0*cbrt(x)^5)
@inline function cbrt_env(x::Float64, y::Float64, z::Float64)
    (x - y) - (cbrt(x) - cbrt(y))/cbrt_deriv(x)
end
@inline function cv_cbrt(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(cbrt, cbrt_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return cbrt(x), cbrt_deriv(x), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, cbrt_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, cbrt_env, xU, 0.0))
    end
    (x <= p) && (return cbrt(x), cbrt_deriv(x), p)
    return dline_seg(cbrt, cbrt_deriv, x, p, xU)..., p
end
@inline function cc_cbrt(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return cbrt(x), cbrt_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(cbrt, cbrt_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, cbrt_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, cbrt_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(cbrt, cbrt_deriv, x, xL, p)..., p)
    return cbrt(x), cbrt_deriv(x), p
end

# Defines interval version of cbrt if necessary
# Copy of recent version from IntervalArithmetic.jl
if VERSION < v"1.3-"
    function cbrt(x::BigFloat, r::RoundingMode)
            setrounding(BigFloat, r) do
                cbrt(x)
            end
        end
    cbrt(a::Interval{Float64}) = atomic(Interval{Float64}, cbrt(big_val(a)))
    function cbrt(a::Interval{BigFloat})
        isempty(a) && return a
        @round(cbrt(a.lo), cbrt(a.hi))
    end
end

# basic method overloading operator (sinh, tanh, atanh, asinh), convexoconcave or concavoconvex
eps_min_dict = Dict{Symbol,Symbol}(:sinh => :xL, :tanh => :xL, :asinh => :xL,
                                 :atanh => :xL, :tan => :xL, :acos => :xU,
                                 :asin => :xL, :atan => :xL, :erf => :xL,
                                 :cbrt => :xL, :erfinv => :xL, :erfc => :xU)
eps_max_dict = Dict{Symbol,Symbol}(:sinh => :xU, :tanh => :xU, :asinh => :xU,
                                 :atanh => :xU, :tan => :xU, :acos => :xL,
                                 :asin => :xU, :atan => :xU, :erf => :xU,
                                 :cbrt => :xU, :erfinv => :xU, :erfc => :xL)

for expri in (:sinh, :tanh, :asinh, :atanh, :tan, :acos, :asin, :atan,
              (:(SpecialFunctions.erf), :erf), :cbrt,
              (:(SpecialFunctions.erfinv), :erfinv),
              (:(SpecialFunctions.erfc), :erfc))
    if expri isa Symbol
        expri_name = expri
        expri_sym = expri
    else
        expri_name = expri[1]
        expri_sym = expri[2]
    end
    expri_kernel = Symbol(String(expri_sym)*"_kernel")
    expri_cv = Symbol("cv_"*String(expri_sym))
    expri_cc = Symbol("cc_"*String(expri_sym))
    eps_min = eps_min_dict[expri_sym]
    eps_max = eps_max_dict[expri_sym]
    @eval @inline function ($expri_kernel)(x::MC{N, T}, y::Interval{Float64},
                            cv_p::Float64, cc_p::Float64) where {N,T<:Union{NS,MV}}
        xL = x.Intv.lo
        xU = x.Intv.hi
        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
        cv_grad = mid_grad(x.cv_grad, x.cc_grad, cv_id)*dcv
        cc_grad = mid_grad(x.cv_grad, x.cc_grad, cc_id)*dcc
        cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
        return MC{N, T}(cv, cc, y, cv_grad, cc_grad, x.cnst), cv_p, cc_p
    end
    @eval @inline function ($expri_kernel)(x::MC{N, Diff}, y::Interval{Float64},
                            cv_p::Float64, cc_p::Float64) where N
        xL = x.Intv.lo
        xU = x.Intv.hi
        midcv, cv_id = mid3(x.cv, x.cc, $eps_min)
        midcc, cc_id = mid3(x.cv, x.cc, $eps_max)
        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
        gcv1, gdcv1, cv_p = $(expri_cv)(x.cv, xL, xU, cv_p)
        gcc1, gdcc1, cc_p = $(expri_cc)(x.cv, xL, xU, cc_p)
        gcv2, gdcv2, cv_p = $(expri_cv)(x.cc, xL, xU, cv_p)
        gcc2, gdcc2, cc_p = $(expri_cc)(x.cc, xL, xU, cc_p)
        cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
        cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
        return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst), cv_p, cc_p
    end
    @eval @inline function ($expri_name)(x::MC{N,T}) where {N, T<:RelaxTag}
        z, tp1, tp2 = ($expri_kernel)(x, ($expri_name)(x.Intv), Inf, Inf)
        return z
    end
end

# cosh convex
@inline cv_cosh(x::Float64, xL::Float64, xU::Float64) = cosh(x), sinh(x)
@inline cc_cosh(x::Float64, xL::Float64, xU::Float64) = dline_seg(cosh, sinh, x, xL, xU)
@inline function cosh_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}}
    xL = x.Intv.lo
    xU = x.Intv.hi
    eps_max = abs(xU) > abs(xL) ?  xU : xL
    eps_min = in(0.0, x) ? 0.0 : (abs(xU) > abs(xL) ?  xL : xU)
    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
    cc, dcc = cc_cosh(midcc, xL, xU)
    cv, dcv = cv_cosh(midcv, xL, xU)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
    cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
    return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function cosh_kernel(x::MC{N,Diff}, y::Interval{Float64}) where N
    xL = x.Intv.lo
    xU = x.Intv.hi
    eps_max = abs(xU) > abs(xL) ?  xU : xL
    eps_min = in(0.0, x) ? 0.0 : (abs(xU) > abs(xL) ?  xL : xU)
    midcc = mid3v(x.cv, x.cc, eps_max)
    midcv = mid3v(x.cv, x.cc, eps_min)
    cc, dcc = cc_cosh(midcc, xL, xU)
    cv, dcv = cv_cosh(midcv, xL, xU)
    gcc1, gdcc1 = cc_cosh(x.cv, xL, xU)
    gcv1, gdcv1 = cv_cosh(x.cv, xL, xU)
    gcc2, gdcc2 = cc_cosh(x.cc, xL, xU)
    gcv2, gdcv2 = cv_cosh(x.cc, xL, xU)
    cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
    cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
    return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline cosh(x::MC) = cosh_kernel(x, cosh(x.Intv))

@inline deg2rad(x::MC, y::Interval{Float64}) = mult_kernel(x, pi/180.0, y)
@inline rad2deg(x::MC, y::Interval{Float64}) = mult_kernel(x, 180.0/pi, y)
@inline deg2rad(x::MC) = deg2rad(x, (pi/180.0)*x.Intv)
@inline rad2deg(x::MC) = rad2deg(x, (180.0/pi)*x.Intv)

# TODO: ADD efficient kernels for below (if applicable)
@inline sec(x::MC)= inv(cos(x))
@inline csc(x::MC)= inv(sin(x))
@inline cot(x::MC)= inv(tan(x))

@inline asec(x::MC) = acos(inv(x))
@inline acsc(x::MC) = asin(inv(x))
@inline acot(x::MC) = atan(inv(x))

@inline sech(x::MC) = inv(cosh(x))
@inline csch(x::MC) = inv(sinh(x))
@inline coth(x::MC) = inv(tanh(x))

@inline acsch(x::MC) = log(sqrt(1.0 + inv(sqr(x))) + inv(x))
@inline acoth(x::MC) = 0.5*(log(1.0 + inv(x)) - log(1.0 - inv(x)))

@inline sind(x::MC) = sin(deg2rad(x))
@inline cosd(x::MC) = cos(deg2rad(x))
@inline tand(x::MC) = tan(deg2rad(x))
@inline secd(x::MC) = inv(cosd(x))
@inline cscd(x::MC) = inv(sind(x))
@inline cotd(x::MC) = inv(tand(x))

@inline asind(x::MC) = rad2deg(asin(x))
@inline acosd(x::MC) = rad2deg(acos(x))
@inline atand(x::MC) = rad2deg(atan(x))
@inline asecd(x::MC) = rad2deg(asec(x))
@inline acscd(x::MC) = rad2deg(acsc(x))
@inline acotd(x::MC) = rad2deg(acot(x))

@inline deg2rad_kernel(x::MC, y::Interval{Float64}) = deg2rad(x,y)
@inline rad2deg_kernel(x::MC, y::Interval{Float64}) = rad2deg(x,y)

@inline sinpi(x::MC) = sin(pi*x)
@inline cospi(x::MC) = cos(pi*x)

@inline erfcinv(x::MC) = erfinv(1.0 - x)

for expri in (:sec, :csc, :cot, :asec, :acsc, :acot, :sech, :csch, :coth,
              :acsch, :acoth, :sind, :cosd, :tand, :secd, :cscd, :cotd,
              :asind, :acosd, :atand, :asecd, :acscd, :acotd, :sinpi, :cospi)

     expri_kernel = Symbol(String(expri)*"_kernel")
     @eval @inline ($expri)(x::MC, y::Interval{Float64}) = ($expri)(x)
 end
 @eval @inline erfcinv_kernel(x::MC, y::Interval{Float64}) = erfcinv(x)

erf(x) = SpecialFunctions.erf(x)
erfc(x) = SpecialFunctions.erfc(x)
erfinv(x) = SpecialFunctions.erfinv(x)
erfcinv(x) = SpecialFunctions.erfcinv(x)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function DiffusionDefinition.diffusion_message(::Val{:FitzHughNagumo})
    println("* * * * * * * * * * * * *")
    println("FitzHugh–Nagumo diffusion")
    println("* * * * * * * * * * * * *")
    println("Description")
    println("------------")
    println("A model developed to mimic the evolution of a neuron's membrane ")
    println("potential. Represented by two-dimensional process (Yₜ, Xₜ)")
    println("dYₜ = (Yₜ - (Yₜ)³ - Xₜ + s)/ϵ dt,")
    println("dXₜ = (γYₜ - Xₜ - β)dt + σdWₜ.")
    println("The model has five parameters:")
    println("\tϵ, s, γ, β, σ")
    println("\nConstructors")
    println("------------")
    println("You may initialize the process either by passing the parameters")
    println("in the order above, as in:")
    println("ϵ, s, γ, β, σ = ...")
    println("P = FitzHughNagumo(ϵ, s, γ, β, σ)")
    println("or by keyword arguments")
    println("p = (ϵ=..., s=..., γ=..., β=..., σ=...)")
    println("P = FitzHughNagumo(;p...)")
    println()
end

@diffusion_process FitzHughNagumo{T} begin
    :dimensions
    process --> 2
    wiener --> 1

    :parameters
    (ϵ, s, γ, β, σ) --> T

    :additional
    constdiff --> true
end

DiffusionDefinition.b(t, x, P::FitzHughNagumo) = @SVector [
    (x[1] - x[2] - x[1]^3 + P.s)/P.ϵ,
    P.γ*x[1] - x[2] + P.β
]

DiffusionDefinition.σ(t, x, P::FitzHughNagumo) = @SMatrix [0.0; P.σ]


text:
"""
$(SIGNATURES)

Clips the solution by the given threshold `λ` and ceils the entries to the corresponding decimal.
"""
function clip_by_threshold!(x::AbstractArray, λ::T, rounding::Bool = true) where T <: Real
    #dplace = ceil(Int, -log10(λ))
    for i in eachindex(x)
        x[i] = abs(x[i]) < λ ? zero(eltype(x)) : x[i]
        #x[i] = rounding ? round(x[i], digits = dplace) : x[i]
    end
    return
end

# Evaluate the results for pareto
G(opt::AbstractOptimizer{T} where T) = f->f[1] < 1 ? Inf : norm(f, 2) # 2*f[1]-2*log(f[2])
G(opt::AbstractSubspaceOptimizer{T} where T) = f->f[1] < 2 ? Inf : norm(f, 2) # 2*f[1]-2*log(f[2])
# Evaluate F
function F(opt::AbstractOptimizer{T} where T)
    f(x, A, y::AbstractArray) = [norm(x, 0); norm(y .- A*x, 2)] # explicit
    f(x, A, y::AbstractArray, λ) = [norm(x, 0, λ); norm(y .- A*x, 2, λ)]
    f(x, A) = [norm(x,0); norm(A*x, 2)] # implicit
    f(x, A, λ::Number) = [norm(x,0,λ); norm(A*x, 2,λ)] # implicit
    return f
end

# Derive the best n linear independent columns of a matrix
function linear_independent_columns(A::AbstractMatrix{T}, rtol::T = convert(T, 0.1)) where T
    iszero(rtol) && return A
    rA = rank(A)
    @static if VERSION < v"1.7.0"
        qr_ = qr(A, Val(true))
    else
        qr_ = qr(A, ColumnNorm())
    end
    r_ = abs.(diag(qr_.R))
    r = findlast(r_ .>= rtol*first(r_))
    r = max(min(r, rA), 1)
    A_ = zeros(T, size(A,1), r)
    idxs = sort(qr_.p[1:r])
    for i in 1:r
        A_[:, i] .= A[:, idxs[i]]
    end
    return A_
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni.
#

get_color(p::Pigment, uv::Vec2d) = throw(MethodError(get_color, p))
#"struct Pigment is abstract and cannot be used in get_color()"

get_color(p::UniformPigment, uv::Vec2d) = p.color

function get_color(p::CheckeredPigment, uv::Vec2d)
    u = floor(uv.u * p.num_steps)
    v = floor(uv.v * p.num_steps)
    ( (u%2) == (v%2) ) ? p.color1 : p.color2
end

function get_color(p::ImagePigment, uv::Vec2d)
    col = floor(uv.u * p.image.width)
    row = floor(uv.v * p.image.height)
    (col < p.image.width) || (col = p.image.width - 1)
    (row < p.image.height) || (row = p.image.height - 1)

    return get_pixel(p.image, convert(Int64, col), convert(Int64, row))
end

"""
    get_color(p::UniformPigment, uv::Vec2d) :: RGB{Float32}
    get_color(p::CheckeredPigment, uv::Vec2d) :: RGB{Float32}
    get_color(p::ImagePigment, uv::Vec2d) :: RGB{Float32}

Return the RGB color of the pigment `p` at the specified (`u`,`v`) coordinates.

See also: [`Pigment`](@ref), [`UniformPigment`](@ref),
[`CheckeredPigment`](@ref), [`ImagePigment`](@ref), [`Vec2d`](@ref)
"""
get_color

##########################################################################################92

evaluate(b::BRDF, n::Normal, in::Vec, out::Vec, uv::Vec2d) = BLACK

evaluate(b::DiffuseBRDF, n::Normal, in::Vec, out::Vec, uv::Vec2d) =
    get_color(b.pigment, uv) * (b.reflectance / pi)

function evaluate(b::SpecularBRDF, n::Normal, in::Vec, out::Vec, uv::Vec2d)
    θ_in = acos(n ⋅ normalize(in))
    θ_out = acos(n ⋅ normalize(out))

    if abs(θ_in - θ_out) < b.theresold_angle_rad
        return get_color(b.pigment, uv)
    else
        return BLACK
    end
end

"""
    evaluate(b::BRDF, n::Normal, in::Vec, out::Vec, uv::Vec2d) :: RGB{Float32}
    evaluate(b::DiffuseBRDF, n::Normal, in::Vec, out::Vec, uv::Vec2d) :: RGB{Float32}
    evaluate(b::SpecularBRDF, n::Normal, in::Vec, out::Vec, uv::Vec2d) :: RGB{Float32}

Return the RGB color with the specified BRDF `b` and spatial
configuation of durface normal `n`, incident ray direction `in`,
leaving ray direction `out`, (`u`,`v`) coordinates on surface.

See also: [`BRDF`](@ref), [`DiffuseBRDF`](@ref), [`SpecularBRDF`](@ref)
[`Normal`](@ref), [`Vec`](@ref), [`Vec2d`](@ref)
"""
evaluate


text:
module Proc
export cubeAnomalies, removeMSC, gapFillMSC, normalizeTS, simpleAnomalies,
  sampleLandPoints, getMSC, filterTSFFT, getNpY, getMedSC, DATfitOnline,
  spatialinterp, extractLonLats, cubePCA, rotation_matrix, transformPCA, explained_variance,exportcube
using ..DAT, ..Cubes

import Dates.year
"""
    getNpY(cube::AbstractCubeData)
    getNpY(cube::InputCube)

Get the number of time steps per year
"""
function getNpY(cube::AbstractCubeData)
    axlist = caxes(cube)
    isTime = [isa(a,TimeAxis) for a in axlist]
    timax = axlist[isTime][1]
    years = year.(timax.values)
    years[end] > years[1] + 1 || error("Must have at least 3 years to calculate MSC")
    return count(i -> i == years[1] + 1, years)
end
getNpY(cube::InputCube)=getNpY(cube.cube)

include("MSC.jl")
include("Stats.jl")
include("CubeIO.jl")
include("TSDecomposition.jl")
include("remap.jl")
include("Shapes.jl")
using .ReSample, .MSC, .Stats, .CubeIO,
  .TSDecomposition

end


text:
#= Print all the subsequences of the string.
=#

## Function

function subsequence(s, ans)
    if (length(s) == 0)
        if (length(ans) >= 1)
            println(ans)
        end
        return
    end
    ch = s[1]
    rest_of_string = SubString(s, 2)
    subsequence(rest_of_string, ans)
    subsequence(rest_of_string, ans * ch)
end

## Input

println("Enter the string: ")
s = readline()

## Calling the Function

println("The subsequences are :")
subsequence(s, "")

#=
Sample Test Case:
Input:
    Enter the string: abc
Output:
    The subsequences are :
    c
    b
    bc
    a
    ac
    ab
    abc
Time complexity: O ( 2^N * l) (where l is the length of string)
=#

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
module MultiQuad

using QuadGK, Cuba, HCubature, Unitful, FastGaussQuadrature, LinearAlgebra, LRUCache

export quad, dblquad, tplquad

function _quadgk_wrapper_1D(arg::Function, x1, x2; kwargs...)
    return quadgk(arg, x1, x2; kwargs...)
end

function _cuba_wrapper_1D(arg::Function, x1, x2; method::Symbol, kwargs...)
    if method == :suave
        integrate = suave
    elseif method == :vegas
        integrate = vegas
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    units = unit(arg(x1)) * unit(x1)

    arg2(a) = ustrip(units, (x2 - x1) * arg((x2 - x1) * a + x1))::Float64

    function integrand(x, f)
        f[1] = arg2(x[1])
    end

    result, err = integrate(integrand, 1, 1; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

#order, method, optional

function get_weights_fastgaussquadrature(order::Int, method::Symbol, optional::Real=NaN)
    if method == :gausslegendre
        xw = gausslegendre
        return xw(order)

    elseif method == :gausshermite
        xw = gausshermite
        return xw(order)

    elseif method == :gausslaguerre
        xw = gausslaguerre
        if isnan(optional)
            return xw(order)
        else
            α = optional
            return xw(order, α)
        end

    elseif method == :gausschebyshev
        xw = gausschebyshev
        if isnan(optional)
            return xw(order)
        else
            kind = optional
            return xw(order, kind)
        end

    # elseif method == :gaussjacobi
    #     xw = gaussjacobi
    #     return xw(order)

    elseif method == :gaussradau
        xw = gaussradau
        return xw(order)

    elseif method == :gausslobatto
        xw = gausslobatto
        return xw(order)

    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

end

const lru_xw_fastgaussquadrature = LRU{
    Tuple{Int,Symbol,Real},
    Tuple{Vector{Float64},Vector{Float64}},
}(maxsize = Int(1e5))

function get_weights_fastgaussquadrature_cached(order::Int, method::Symbol, optional::Real=NaN)
    get!(lru_xw_fastgaussquadrature, (order, method, optional)) do
        get_weights_fastgaussquadrature(order, method, optional)
    end
end

function _fastgaussquadrature_wrapper_1D(arg::Function, x1, x2; method::Symbol, order::Int, kwargs...)
    if ! (method in [:gausslegendre, :gausschebyshev, :gaussradau, :gausslobatto])
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    if haskey(kwargs, :α)
        optional = kwargs[:α]
    elseif haskey(kwargs, :kind)
        optional = kwargs[:kind]
    else
        optional = NaN
    end

    x, w = get_weights_fastgaussquadrature_cached(order, method, optional)

    b=x2
    a=x1

    b_a_2 = (b-a)/2
    a_plus_b_2 = (a+b)/2

    arg_gauss(x) = arg(b_a_2*x + a_plus_b_2)

    integral = b_a_2 * dot(w, arg_gauss.(x))
    return integral, NaN

end

function _fastgaussquadrature_wrapper_1D(arg::Function; method::Symbol, order::Int, kwargs...)
    if ! (method in [:gausshermite , :gausslaguerre])
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    if haskey(kwargs, :α)
        optional = kwargs[:α]
    elseif haskey(kwargs, :kind)
        optional = kwargs[:kind]
    else
        optional = NaN
    end

    x, w = get_weights_fastgaussquadrature_cached(order, method, optional)

    integral = dot(w, arg.(x))
    return integral, NaN

end

@doc raw"""
    quad(arg::Function, x1, x2[; method = :quadgk, oder=1000, kwargs...])

Performs the integral ``\int_{x1}^{x2}f(x)dx``

Available integration methods:
- `:suave`
- `:vegas`
- `:quadgk`

There are several fixed order quadrature methods available based on [`FastGaussQuadrature`](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)

- `:gausslegendre`
- `:gausshermite`
- `:gausslaguerre`
- `:gausschebyshev`
- `:gaussradau`
- `:gausslobatto`

If a specific integration routine is not needed, it is suggested to use `:quadgk` (the default option).
For highly oscillating integrals, it is advised to use `:gausslegendre`with a high order (~10000).

Please note that `:gausshermite` and `:gausslaguerre` are used to specific kind of integrals with infinite bounds.

For more informations on these integration techniques, please follow the official documentation [here](https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/gaussquadrature/).

See [QuadGK](https://github.com/JuliaMath/QuadGK.jl) and [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.
# Examples
```jldoctest
using MultiQuad

func(x) = x^2*exp(-x)
quad(func, 0, 4) # equivalent to quad(func, 0, 4, method=:quadgk)
quad(func, 0, 4, method=:gausslegendre, order=1000)

# for certain kinds of integrals with infinite bounds, it may be possible to use a specific integration routine
func(x) = x^2*exp(-x)
quad(func, 0, Inf, method=:quadgk)[1]
# gausslaguerre computes the integral of f(x)*exp(-x) from 0 to infinity
quad(x -> x^4, method=:gausslaguerre, order=10000)[1]
```

`quad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:

```jldoctest
using Unitful

func(x) = x^2
integral, error = quad(func, 1u"m", 5u"m")
```
"""
function quad(arg::Function, x1, x2; method::Symbol = :quadgk, order::Int=1000, kwargs...)

    if method in [:suave, :vegas]
        return _cuba_wrapper_1D(arg, x1, x2; method=method)
    elseif method == :quadgk
        return _quadgk_wrapper_1D(arg, x1, x2; kwargs...)
    elseif method in [:gausslegendre, :gausschebyshev, :gaussradau, :gausslobatto]
        return _fastgaussquadrature_wrapper_1D(arg, x1, x2; method=method, order=order, kwargs...)
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end
end

@doc raw"""
    quad(arg::Function; method[, oder=1000, kwargs...])
"""
function quad(arg::Function; method::Symbol, order::Int=1000, kwargs...)

    if method in [:gausshermite , :gausslaguerre]
        return _fastgaussquadrature_wrapper_1D(arg; method=method, order=order, kwargs...)
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end
end

function _cuba_wrapper_2D(arg::Function, x1, x2, y1::Function, y2::Function; method::Symbol, kwargs...)
    if method == :cuhre
        integrate = cuhre
    elseif method == :divonne
        integrate = divonne
    elseif method == :suave
        integrate = suave
    elseif method == :vegas
        integrate = vegas
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))

    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)

    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64

    function integrand2(x, f)
        f[1] = arg2(x[1], x[2])
    end

    result, err = integrate(integrand2, 2, 1; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

function _cuba_wrapper_2D(arg::Function, x1, x2, y1, y2; method::Symbol, kwargs...)
    if method == :cuhre
        integrate = cuhre
    elseif method == :divonne
        integrate = divonne
    elseif method == :suave
        integrate = suave
    elseif method == :vegas
        integrate = vegas
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    units = unit(arg(y1, x1)) * unit(x1) * unit(y1)

    arg1(a, x) = (y2 - y1) * arg((y2 - y1) * a + y1, x)

    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64

    function integrand2(x, f)
        f[1] = arg2(x[1], x[2])
    end

    result, err = integrate(integrand2, 2, 1; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

function _hcubature_wrapper_2D(arg::Function, x1, x2, y1::Function, y2::Function; kwargs...)
    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))

    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)

    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64

    function integrand(arr)
        return arg2(arr[1], arr[2])
    end

    min_arr = [0, 0]
    max_arr = [1, 1]
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

function _hcubature_wrapper_2D(arg::Function, x1, x2, y1, y2; kwargs...)
    units = unit(arg(y1, x1)) * unit(x1) * unit(y1)

    arg1(a, x) = (y2 - y1) * arg((y2 - y1) * a + y1, x)

    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64

    function integrand(arr)
        return arg2(arr[1], arr[2])
    end

    min_arr = [0, 0]
    max_arr = [1, 1]
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

@doc raw"""
    dblquad(arg::Function, x1, x2, y1::Function, y2::Function; method = :cuhre, kwargs...)

Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}f(y,x)dydx``

Available integration methods:
- `:cuhre`
- `:divonne`
- `:suave`
- `:vegas`
- `:hcubature`

See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.

# Examples
```jldoctest
func(y,x) = sin(x)*y^2
integral, error = dblquad(func, 1, 2, x->0, x->x^2, rtol=1e-9)
```

`dblquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:

```jldoctest
using Unitful

func(y,x) = x*y^2
integral, error = dblquad(func, 1u"m", 2u"m", x->0u"m^2", x->x^2, rtol=1e-9)
```
"""
function dblquad(
    arg::Function,
    x1,
    x2,
    y1::Function,
    y2::Function;
    method = :hcubature,
    kwargs...,
)

    if method in [:cuhre, :divonne, :suave, :vegas]
        return _cuba_wrapper_2D(arg, x1, x2, y1, y2; method=method)
    elseif method == :hcubature
        return _hcubature_wrapper_2D(arg, x1, x2, y1, y2; kwargs...)
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end
end

function dblquad(
    arg::Function,
    x1,
    x2,
    y1,
    y2;
    method = :hcubature,
    kwargs...,
)

    if method in [:cuhre, :divonne, :suave, :vegas]
        return _cuba_wrapper_2D(arg, x1, x2, y1, y2; method=method)
    elseif method == :hcubature
        return _hcubature_wrapper_2D(arg, x1, x2, y1, y2; kwargs...)
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end
end

function _cuba_wrapper_3D(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function;  method::Symbol, kwargs...)
    if method == :cuhre
        integrate = cuhre
    elseif method == :divonne
        integrate = divonne
    elseif method == :suave
        integrate = suave
    elseif method == :vegas
        integrate = vegas
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
            unit(z1(y1(x1), x1))

    arg0(a, y, x) =
        (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)

    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)

    arg2(a, b, c) =
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64

    function integrand2(x, f)
        f[1] = arg2(x[1], x[2], x[3])
    end

    result, err = integrate(integrand2, 3, 1; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

function _cuba_wrapper_3D(arg::Function, x1, x2, y1, y2, z1, z2;  method::Symbol, kwargs...)
    if method == :cuhre
        integrate = cuhre
    elseif method == :divonne
        integrate = divonne
    elseif method == :suave
        integrate = suave
    elseif method == :vegas
        integrate = vegas
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    units = unit(arg(z1, y1, x1)) * unit(x1) * unit(y1) * unit(z1)

    arg0(a, y, x) = (z2 - z1) * arg((z2 - z1) * a + z1, y, x)

    arg1(a, b, x) = (y2 - y1) * arg0(a, (y2 - y1) * b + y1, x)

    arg2(a, b, c) =
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64

    function integrand2(x, f)
        f[1] = arg2(x[1], x[2], x[3])
    end

    result, err = integrate(integrand2, 3, 1; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

function _hcubature_wrapper_3D(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; kwargs...)

    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
            unit(z1(y1(x1), x1))

    arg0(a, y, x) =
        (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)

    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)

    arg2(a, b, c) =
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64

    function integrand(arr)
        return arg2(arr[1], arr[2], arr[3])
    end

    min_arr = [0, 0, 0]
    max_arr = [1, 1, 1]
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

function _hcubature_wrapper_3D(arg::Function, x1, x2, y1, y2, z1, z2; kwargs...)

    units = unit(arg(z1, y1, x1)) * unit(x1) * unit(y1) * unit(z1)

    arg0(a, y, x) = (z2 - z1) * arg((z2 - z1) * a + z1, y, x)

    arg1(a, b, x) = (y2 - y1) * arg0(a, (y2 - y1) * b + y1, x)

    arg2(a, b, c) =
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64

    function integrand(arr)
        return arg2(arr[1], arr[2], arr[3])
    end

    min_arr = [0, 0, 0]
    max_arr = [1, 1, 1]
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

@doc raw"""
    tplquad(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; method = :cuhre, kwargs...)

Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}\int_{z1(x,y)}^{z2(x,y)}f(z,y,x)dzdydx``

Available integration methods:
- `:cuhre`
- `:divonne`
- `:suave`
- `:vegas`
- `:hcubature`

See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.

# Examples
```jldoctest
func(z,y,x) = sin(z)*y*x
integral, error = tplquad(func, 0, 4, x->x, x->x^2, (x,y)->2, (x,y)->3*x)
```

`tplquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:

```jldoctest
using Unitful

func(z,y,x) = sin(z)*y*x
integral, error = tplquad(func, 0u"m", 4u"m", x->0u"m^2", x->x^2, (x,y)->0, (x,y)->3)
```
"""
function tplquad(
    arg::Function,
    x1,
    x2,
    y1::Function,
    y2::Function,
    z1::Function,
    z2::Function;
    method = :hcubature,
    kwargs...,
)

    if method in [:cuhre, :divonne, :suave, :vegas]
        return _cuba_wrapper_3D(arg, x1, x2, y1, y2, z1, z2;  method=method, kwargs...)

    elseif method == :hcubature
        return _hcubature_wrapper_3D(arg, x1, x2, y1, y2, z1, z2; kwargs...)

    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end
end

function tplquad(
    arg::Function,
    x1,
    x2,
    y1,
    y2,
    z1,
    z2;
    method = :hcubature,
    kwargs...,
)

    if method in [:cuhre, :divonne, :suave, :vegas]
        return _cuba_wrapper_3D(arg, x1, x2, y1, y2, z1, z2;  method=method, kwargs...)

    elseif method == :hcubature
        return _hcubature_wrapper_3D(arg, x1, x2, y1, y2, z1, z2; kwargs...)

    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end
end

end # module


text:
module Ex5_11c
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.11c - Solution of an AFDIP
println("Example 5.11c")

# define s as an improper transfer function
s = rtf('s');
# define Gu(s), Gf(s) and Gw(s)
Gu = [(s+1)/(s+2); (s+2)/(s+3)];     # enter Gu(s)
Gw = [1/(s+2); 0];                   # enter Gw(s)
Gf = [(s+1)/(s+2) 0; 0 1];           # enter Gf(s)
p = 2; mu = 1; mw = 1; mf = 2;       # enter dimensions
SFDI = eye(mf) .> 0

# build the synthesis model with additive faults
sysf = fdimodset(dss([Gu Gf Gw]), c = 1:mu,f = mu .+ (1:mf), n = (mu+mf) .+ (1:mw));

Q, R, info = afdisyn(sysf, SFDI; smarg =-3, sdeg = -3)

# scale with -1 to match example
Q.sys[1] = -Q.sys[1]; R.sys[1] = -R.sys[1];

for i = 1:size(SFDI,1)
    println("Q[$i] = $(dss2rm(Q.sys[i], atol = 1.e-7))")
    println("Rf[$i] = $(dss2rm(R.sys[i][:,R.faults], atol = 1.e-7))")
    println("Rw[$i] = $(dss2rm(R.sys[i][:,R.noise], atol = 1.e-7))")
end

# check synthesis conditions: Qt*[Gu Gd;I 0] = 0 and Qt*[Gf; 0] = Rf
Rt = fdIFeval(Q,sysf) # form Qt*[Gu Gd Gf;I 0 0];
@test iszero(vcat(Rt.sys...)[:,[Rt.controls;Rt.disturbances]],atol=1.e-7) &&
      iszero(vcat(Rt.sys...)[:,[Rt.faults;Rt.noise]]-vcat(R.sys...),atol=1.e-7) &&
      fdif2ngap(R,SFDI)[1] ≈ info.gap

println("Q = ")
display(Q)
println("R = ")
display(R)
println("gap = $(fdif2ngap(R,SFDI)[1])")

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
"""
    This package contains a number of kernels frequently used in SPH, as well as their functions and derivatives an 2D and 3D.
"""

module SPHKernels

    export  kernel_value_1D,     𝒲₁,
            kernel_value_2D,     𝒲₂,
            kernel_value_3D,     𝒲₃,
            kernel_deriv_1D,    ∇𝒲₁,
            kernel_deriv_2D,    ∇𝒲₂,
            kernel_deriv_3D,    ∇𝒲₃,
            bias_correction_1D, δρ₁,
            bias_correction_2D, δρ₂,
            bias_correction_3D, δρ₃,
            SPHKernel,
            Cubic,
            Quintic,
            WendlandC2,
            WendlandC4,
            WendlandC6,
            WendlandC8

    """
        SPHKernel

    Supertype for all SPH kernels.
    """
    abstract type SPHKernel end

    include("bsplines/Cubic.jl")
    include("bsplines/Quintic.jl")
    include("wendland/C2.jl")
    include("wendland/C4.jl")
    include("wendland/C6.jl")
    include("wendland/C8.jl")

    # multiple dispatch for nicer look

    """
        𝒲₁( kernel::SPHKernel, u::Real, h_inv::Real)

    Evaluate 1D spline at position ``u = \\frac{x}{h}``.
    """
    𝒲₁( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_1D(kernel, u, h_inv)

    """
        ∇𝒲₁( kernel::SPHKernel, u::Real, h_inv::Real)

    Evaluate 1D derivative at position ``u = \\frac{x}{h}``.
    """
    ∇𝒲₁(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_1D(kernel, u, h_inv)

    """
        𝒲₂( kernel::SPHKernel, u::Real, h_inv::Real)

    Evaluate 2D spline at position ``u = \\frac{x}{h}``.
    """
    𝒲₂( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_2D(kernel, u, h_inv)

    """
        ∇𝒲₂( kernel::SPHKernel, u::Real, h_inv::Real)

    Evaluate 1D derivative at position ``u = \\frac{x}{h}``.
    """
    ∇𝒲₂(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_2D(kernel, u, h_inv)

    """
        𝒲₃( kernel::SPHKernel, u::Real, h_inv::Real)

    Evaluate 3D spline at position ``u = \\frac{x}{h}``.
    """
    𝒲₃( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_3D(kernel, u, h_inv)

    """
        ∇𝒲₃( kernel::SPHKernel, u::Real, h_inv::Real)

    Evaluate 1D derivative at position ``u = \\frac{x}{h}``.
    """
    ∇𝒲₃(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_3D(kernel, u, h_inv)

    """
        δρ₁(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)

    Corrects the 1D density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
    """
    δρ₁(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_1D(kernel, density, m, h_inv)

    """
        δρ₂(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)

    Corrects the 2D density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
    """
    δρ₂(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_2D(kernel, density, m, h_inv)

    """
        δρ₃(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)

    Corrects the 3D density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
    """
    δρ₃(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_3D(kernel, density, m, h_inv)

end # module