text:
# Context management

export
    CuPrimaryContext, CuContext, CuCurrentContext, activate,
    unsafe_reset!, isactive, flags, setflags!,
    device_synchronize, CuCurrentDevice,

## construction and destruction

@enum_without_prefix CUctx_flags CU_

"""
    CuPrimaryContext(dev::CuDevice)

Create a primary CUDA context for a given device.

Each primary context is unique per device and is shared with CUDA runtime API. It is meant
for interoperability with (applications using) the runtime API.
"""
struct CuPrimaryContext
    dev::CuDevice
end

"""
    CuContext(dev::CuDevice, flags=CTX_SCHED_AUTO)
    CuContext(f::Function, ...)

Create a CUDA context for device. A context on the GPU is analogous to a process on the CPU,
with its own distinct address space and allocated resources. When a context is destroyed,
the system cleans up the resources allocated to it.

When you are done using the context, call [`CUDA.unsafe_destroy!`](@ref) to mark it for
deletion, or use do-block syntax with this constructor.

"""
mutable struct CuContext
    handle::CUcontext
    valid::Bool

    function new_unique(handle)
        Base.@lock context_lock get!(valid_contexts, handle) do
            new(handle, true)
        end
    end

    function CuContext(dev::CuDevice, flags=0)
        handle_ref = Ref{CUcontext}()
        cuCtxCreate_v2(handle_ref, flags, dev)
        new_unique(handle_ref[])
    end

    """
        CuContext(pctx::CuPrimaryContext)

    Retain the primary context on the GPU, returning a context compatible with the driver API.
    The primary context will be released when the returned driver context is finalized.

    As these contexts are refcounted by CUDA, you should not call [`CUDA.unsafe_destroy!`](@ref)
    on them but use [`CUDA.unsafe_release!`](@ref) instead (available with do-block syntax as
    well).
    """
    function CuContext(pctx::CuPrimaryContext)
        handle_ref = Ref{CUcontext}()
        cuDevicePrimaryCtxRetain(handle_ref, pctx.dev)
        return new_unique(handle_ref[])
    end

    """
        CuCurrentContext()

    Return the current context, or `nothing` if there is no active context.
    """
    global function CuCurrentContext()
        handle_ref = Ref{CUcontext}()
        cuCtxGetCurrent(handle_ref)
        if handle_ref[] == C_NULL
            return nothing
        else
            new_unique(handle_ref[])
        end
    end

    """
        CuContext(ptr)

    Identify the context a CUDA memory buffer was allocated in.
    """
    function CuContext(x::Union{Ptr,CuPtr})
        new_unique(attribute(CUcontext, x, POINTER_ATTRIBUTE_CONTEXT))
    end
end

# the `valid` bit serves two purposes: make sure we don't double-free a context (in case we
# early-freed it ourselves before the GC kicked in), and to make sure we don't free derived
# resources after the owning context has been destroyed (which can happen due to
# out-of-order finalizer execution)
const valid_contexts = Dict{CUcontext,CuContext}()
const context_lock = ReentrantLock()
isvalid(ctx::CuContext) = ctx.valid
# NOTE: we can't just look up by the handle, because contexts derived from a primary one
#       have the same handle even though they might have been destroyed in the meantime.
function invalidate!(ctx::CuContext)
    Base.@lock context_lock delete!(valid_contexts, ctx.handle)
    ctx.valid = false
    return
end

"""
    unsafe_destroy!(ctx::CuContext)

Immediately destroy a context, freeing up all resources associated with it. This does not
respect any users of the context, and might make other objects unusable.
"""
function unsafe_destroy!(ctx::CuContext)
    if isvalid(ctx)
        cuCtxDestroy_v2(ctx)
        invalidate!(ctx)
    end
end

Base.unsafe_convert(::Type{CUcontext}, ctx::CuContext) = ctx.handle

# NOTE: we don't implement `isequal` or `hash` in order to fall back to `===` and `objectid`
#       as contexts are unique, and with primary device contexts identical handles might be
#       returned after resetting the context (device) and all associated resources.

function Base.show(io::IO, ctx::CuContext)
    if ctx.handle != C_NULL
        fields = [@sprintf("%p", ctx.handle), @sprintf("instance %x", objectid(ctx))]
        if !isvalid(ctx)
            push!(fields, "invalidated")
        end

        print(io, "CuContext(", join(fields, ", "), ")")
    else
        print(io, "CuContext(NULL)")
    end
end

## core context API

"""
    push!(CuContext, ctx::CuContext)

Pushes a context on the current CPU thread.
"""
Base.push!(::Type{CuContext}, ctx::CuContext) = cuCtxPushCurrent_v2(ctx)

"""
    pop!(CuContext)

Pops the current CUDA context from the current CPU thread.
"""
function Base.pop!(::Type{CuContext})
    handle_ref = Ref{CUcontext}()
    cuCtxPopCurrent_v2(handle_ref)
    # we don't return the context here, because it may be unused
    # (and constructing the unique object is expensive)
end

# perform some finalizer actions in a context
macro finalize_in_ctx(ctx, body)
    # XXX: should this not integrate with the high-level context management from state.jl?
    #      it might be good that the driver API wrappers don't need that runtime-esque
    #      state management, but it might be confusing that `context()` doesn't work here.
    quote
        ctx = $(esc(ctx))
        if isvalid(ctx)
            push!(CuContext, ctx)
            try
                $(esc(body))
            finally
                pop!(CuContext)
            end
        end
    end
end

"""
    activate(ctx::CuContext)

Binds the specified CUDA context to the calling CPU thread.
"""
activate(ctx::CuContext) = cuCtxSetCurrent(ctx)

function CuContext(f::Function, dev::CuDevice, args...)
    ctx = CuContext(dev, args...)    # implicitly pushes
    try
        f(ctx)
    finally
        pop!(CuContext)
        unsafe_destroy!(ctx)
    end
end

## primary context management

"""
    CUDA.unsafe_release!(ctx::CuContext)

Lower the refcount of a context, possibly freeing up all resources associated with it. This
does not respect any users of the context, and might make other objects unusable.
"""
function unsafe_release!(ctx::CuContext)
    if isvalid(ctx)
        dev = CuDevice(ctx)
        pctx = CuPrimaryContext(dev)
        if version() >= v"11"
            cuDevicePrimaryCtxRelease_v2(dev)
        else
            cuDevicePrimaryCtxRelease(dev)
        end
        isactive(pctx) || invalidate!(ctx)
    end
    return
end

function CuContext(f::Function, pctx::CuPrimaryContext)
    ctx = CuContext(pctx)
    try
        f(ctx)
    finally
        unsafe_release!(ctx)
    end
end

"""
    unsafe_reset!(pctx::CuPrimaryContext)

Explicitly destroys and cleans up all resources associated with a device's primary context
in the current process. Note that this forcibly invalidates all contexts derived from this
primary context, and as a result outstanding resources might become invalid.
"""
function unsafe_reset!(pctx::CuPrimaryContext)
    ctx = CuContext(pctx)
    invalidate!(ctx)
    if version() >= v"11"
        cuDevicePrimaryCtxReset_v2(pctx.dev)
    else
        cuDevicePrimaryCtxReset(pctx.dev)
    end
    return
end

function state(pctx::CuPrimaryContext)
    flags = Ref{Cuint}()
    active = Ref{Cint}()
    cuDevicePrimaryCtxGetState(pctx.dev, flags, active)
    return (flags[], active[] == one(Cint))
end

"""
    isactive(pctx::CuPrimaryContext)

Query whether a primary context is active.
"""
isactive(pctx::CuPrimaryContext) = state(pctx)[2]

"""
    flags(pctx::CuPrimaryContext)

Query the flags of a primary context.
"""
flags(pctx::CuPrimaryContext) = state(pctx)[1]

"""
    setflags!(pctx::CuPrimaryContext)

Set the flags of a primary context.
"""
function setflags!(pctx::CuPrimaryContext, flags)
    if version() >= v"11"
        cuDevicePrimaryCtxSetFlags_v2(pctx.dev, flags)
    else
        cuDevicePrimaryCtxSetFlags(pctx.dev, flags)
    end
end

## context properties

"""
    CuCurrentDevice()

Returns the current device, or `nothing` if there is no active device.
"""
function CuCurrentDevice()
    device_ref = Ref{CUdevice}()
    res = unsafe_cuCtxGetDevice(device_ref)
    if res == ERROR_INVALID_CONTEXT
        return nothing
    elseif res != SUCCESS
        throw_api_error(res)
    end
    return CuDevice(Bool, device_ref[])
end

"""
    CuDevice(::CuContext)

Returns the device for a context.
"""
function CuDevice(ctx::CuContext)
    push!(CuContext, ctx)
    dev = CuCurrentDevice()
    pop!(CuContext)
    return dev
end

"""
    synchronize(ctx::Context)

Synchronize the current context, waiting for all outstanding operations to complete.

!!! warning

    This is an operation that blocks in the driver, and should be avoided if possible.
    Instead, use [`device_synchronize()`](@ref) to perform synchronization in Julia.
"""
function synchronize(ctx::CuContext)
    push!(CuContext, ctx)
    try
        cuCtxSynchronize()
        check_exceptions()
    finally
        pop!(CuContext)
    end
end

## cache config

export cache_config, cache_config!

@enum_without_prefix CUfunc_cache CU_

function cache_config()
    config = Ref{CUfunc_cache}()
    cuCtxGetCacheConfig(config)
    return config[]
end

function cache_config!(config::CUfunc_cache)
    cuCtxSetCacheConfig(config)
end

## shared memory config

export shmem_config, shmem_config!

@enum_without_prefix CUsharedconfig CU_

function shmem_config()
    config = Ref{CUsharedconfig}()
    cuCtxGetSharedMemConfig(config)
    return config[]
end

function shmem_config!(config::CUsharedconfig)
    cuCtxSetSharedMemConfig(config)
end

## limits

export limit, limit!

@enum_without_prefix CUlimit CU_

function limit(lim::CUlimit)
    val = Ref{Csize_t}()
    cuCtxGetLimit(val, lim)
    return Int(val[])
end

limit!(lim::CUlimit, val) = cuCtxSetLimit(lim, val)


text:
global const COLSET = ["values", "movename", "moveparam", "accepted", "timestamp"]

"""
    readst(fname)

Read a search string from a csv-file
"""
function readst(fname)
    @assert endswith(fname, ".csv")
    st = CSV.read(fname, DataFrame)
    c = copy(colset)
    for v in names(st)
        if v ∈ colset
            @info "Found column $v"
            pop!(c)
        end
    end
    if !isempty(s)
        s = "Following columns were not found:\n"
        for v in c
            s * "- $(v)\n"
        end
        s = "Currently, PolyaUtils relies on the existence of columns with the aforementioned names. Not all functionality will be available to you with current naming. Either provide a new csv-file, or rename the appropriate columns using `rename(st, namedict)`. See the documentation of the `rename`-function in DataFrames for more information."
        @info s

    end
    return st
end

"""
    compute_deltas(values, acceptance)

"""
function compute_deltas(values, acceptance)
    # maintain a current value, which is updated every time a move is accepted
    v_curr = first(values)
    Δ = Vector{eltype(values)}(undef, length(values))
    for (i,v) in enumerate(values)
        Δ[i] = v - v_curr
        if acceptance[i]
            v_curr = v
        end
    end
    return Δ
end

#=
function compute_tsli(times, acceptance)
    t_curr = first(times)
end=#
"""
    vcurr(df)

Compute the current state objective value for each step in the search trace
"""
function vcurr(df)
    I = Iteratros.zip(df.accepted, df.values)
    accumulate( (c,n) -> n[1] ? n[2] : c, I, init  = first(df.values))
end

"""
    vbest(df)

Compute the best state objective value for each step in the search trace
"""
vbest(df) = accumulate(min, df.vcurr, init = first(df.vcurr))

"""
    add_columns!(df)

Adds the columns required to generate plots. Following columns are added:
- iteration: Iteration number in the trace, where each iteration is associated with exactly one move
- delta: Difference in objective value between the current solution and the proposed alternative
-
"""
function add_columns!(df)
    V = eltype(df.value)
    df.iteration = 1:size(df,1)
    df.vcurr = [first(df.value)] * vcurr(df) # Current objective value
    df.vbest = [first(df.value)] * vbest(df) # Best objective value
    df.Δv = df.vcurr - df.value              # Objective delta
    rankdict = rankingdict(df.values)
    ranker = v -> rankdict[v]
    df["ranking"] = ranker.(df.values)       # Rank the
    return df
end

finalcols = ["iteration", "movetype", "moveparam", "value", "accepted", "timestamp", "improvementstep"]

improving(df) =  df[df.Δv .<= zero(eltype(df.Δv)) .& df.accepted, :]
strictlyimproving(df) = improving(df)[df.Δv .< zero(eltype(df.Δv)), :]

ext: jl

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz


text:
@testset "Slow Manifold System" begin
  function slow_manifold(du, u, p, t)
    du[1] = p[1]*u[1]
    du[2] = p[2]*(u[2]-u[1]^2) + p[3]*sin(t^2)
  end

  u0 = [3.0; -2.0]
  tspan = (0.0, 5.0)
  p = [-2.0; -0.5; 1.0]

  problem = ODEProblem(slow_manifold, u0, tspan, p)
  solution = solve(problem, Tsit5(), saveat = 0.01)

  ufun(u,p,t) = sin(t^2)

  prob = ContinuousDataDrivenProblem(solution, U = ufun)

  @variables u[1:2] y[1:1] t
  Ψ = Basis([u; u[1]^2; y], u, controls = y, iv = t)

  for alg in [DMDPINV(); DMDSVD(); TOTALDMD()]
    res = solve(prob, Ψ, alg, digits = 1)
    b = result(res)
    m = metrics(res)
    @test isapprox(eigvals(b), [2*p[1]; p[1]; p[2]], atol = 1e-1)
    @test all(m[:L₂] .< 3e-1)

    # TODO This does not work right now, but it should
    #sdict = Dict([y[1] => sin(t^2)])

    #for (i,eq) in enumerate(equations(b))
    #  lhs, rhs = eq.lhs, eq.rhs
    #  b[i] = Num(lhs) ~ substitute(Num(rhs), sdict)
    #end

    #_prob = ODEProblem((args...)->b(args...), u0, tspan, parameters(res))
    #_sol = solve(_prob, Tsit5(), saveat = solution.t)
    #@test norm(solution - _sol)/size(X, 2) < 5e-1
  end
end

@testset "Slow Manifold Discrete System" begin
  function slow_manifold(du, u, p, t)
    du[1] = p[1]*u[1]
    du[2] = p[2]*u[2]+(p[1]^2 - p[2])*u[1]^2 + p[3]*sin(t^2)
  end

  u0 = [3.0; -2.0]
  tspan = (0.0, 10.0)
  p = [-0.8; -0.7; 1.0]

  problem = DiscreteProblem(slow_manifold, u0, tspan, p)
  solution = solve(problem, FunctionMap())

  ufun(u,p,t) = sin(t^2)

  ddprob = DiscreteDataDrivenProblem(solution, U = ufun)

  @variables u[1:2] y[1:1] t

  Ψ = Basis([u[1]; u[1]^2; u[2]-u[1]^2; y], u, controls = y, iv = t)

  for alg in [DMDPINV(); DMDSVD(); TOTALDMD()]
    res = solve(ddprob, Ψ, alg, digits = 1)
    b = result(res)
    m = metrics(res)
    @test isapprox(eigvals(b), [p[1]; p[2]; p[1]^2], atol = 1e-1)
    @test all(m[:L₂] .< 3e-1)
  end
end


text:
"""
function ptcsoli( F!, x0, FS, FPS, Jvec = dirder; rtol = 1.e-6, atol = 1.e-12,
    maxit = 20, lmaxit = -1, lsolver = "gmres", eta = 0.1, fixedeta = true,
    Pvec = nothing, PvecKnowsdelta = false, pside = "right", delta0 = 1.e-6,
    dx = 1.e-7, pdata = nothing, printerr = true, keepsolhist = false,
)

C. T. Kelley, 2021

Julia versions of the nonlinear solvers from my SIAM books.
Herewith: some new stuff ==> ptcsoli

PTC finds the steady-state solution of u' = -F(u), u(0) = u_0.
The - sign is a convention.

You must allocate storage for the function and Krylov basis in advance
--> in the calling program <-- ie. in FS and FPS

Inputs:\n
- F!: function evaluation, the ! indicates that F! overwrites FS, your
    preallocated storage for the function.\n
    So, FV=F!(FV,x) or FV=F!(FV,x,pdata) returns FV=F(x)

- x0: initial iterate\n

- FS: Preallocated storage for function. It is an N x 1 column vector.\n
You may dimension it as (n,) or (n,1). (n,) is best, but the
solvers can deal with it either way.

- FPS: preallocated storage for the Krylov basis. It is an N x m matrix where
       you plan to take at most m-1 GMRES iterations before a restart. \n

- Jvec: Jacobian vector product, If you leave this out the
    default is a finite difference directional derivative.\n
    So, FP=Jvec(v,FS,x) or FP=Jvec(v,FS,x,pdata) returns FP=F'(x) v. \n
    (v, FS, x) or (v, FS, x, pdata) must be the argument list,
    even if FP does not need FS.
    One reason for this is that the finite-difference derivative
    does and that is the default in the solver.

- Precision: Lemme tell ya 'bout precision. I designed this code for
    full precision functions and linear algebra in any precision you want.
    You can declare FPS as Float64 or Float32 and ptcsoli
    will do the right thing. Float16 support is there, but not working well.

    If the Jacobian is reasonably well conditioned, you can cut the cost
    of orthogonalization and storage (for GMRES) in half with no loss.
    There is no benefit if your linear solver is not GMRES or if
    othogonalization and storage of the Krylov vectors is only a
    small part of the cost of the computation. So if your preconditioner
    is good and you only need a few Krylovs/Newton, reduced precision won't
    help you much.

    BiCGSTAB does not benefit from reduced precsion.

----------------------

Keyword Arguments (kwargs):\n

rtol and atol: relative and absolute error tolerances\n

delta0: initial pseudo time step. The default value of 1.e-3 is a bit conservative
and is one option you really should play with. Look at the example
where I set it to 1.0!\n

maxit: limit on nonlinear iterations, default=100. \n
This is coupled to delta0. If your choice of delta0 is too small (conservative)
then you'll need many iterations to converge and will need a larger
value of maxit

For PTC you'll need more iterations than for a straight-up
nonlinear solve. This is part of the price for finding the
stable solution. \n

lmaxit: limit on linear iterations. If lmaxit > m-1, where FPS has
m columns, and you need more than m-1 linear iterations, then GMRES
will restart.

The default is -1. For GMRES this means that you'll take m-1 iterations, where
size(V) = (n,m), and get no restarts. For BiCGSTAB you'll then get the default
of 10 iterations.

lsolver: the linear solver, default = "gmres"\n
Your choices will be "gmres" or "bicgstab". However,
gmres is the only option for now. \n

eta and fixed eta: eta > 0 or there's an error.

The linear solver terminates when ||F'(x)s + F(x) || <= etag || F(x) ||

where

etag = eta if fixedeta=true

etag = Eisenstat-Walker as implemented in book if fixedeta=false

The default, which may change, is eta=.1, fixedeta=true \n

Pvec: Preconditioner-vector product. The rules are similar to Jvec
    So, Pv=Pvec(v,x) or Pv=Pvec(v,x,pdata) returns P(x) v where
    P(x) is the preconditioner. You must use x as an input even
    if your preconditioner does not depend on x.\n

PvecKnowsdelta: If you want your preconditioner-vector product to depend on
    the pseudo-timestep delta, put an array deltaval in your precomputed
    data. Initialize it as
    deltaval = zeros(1,)
    and let ptcsoli know about it by setting the kwarg
    PvecKnowsdelta = true
    ptcsoli will update the value in deltaval with every change
    to delta with pdata.deltaval[1]=delta
    so your preconditioner-vector product can get to it.\n

pside: apply preconditioner on pside, default = "right". I do not
      recommend "left". The problem with "left" for ptcsoli is
      that it can fail to satisfy the inexact Newton condition for
      the unpreconditioned equation, especially early in the iteration
      and lead to an incorrect result (unstable solution or wrong
      branch of steady state).
      See Chapter 3 for the story on this. \n

dx: default = 1.e-7\n
difference increment in finite-difference derivatives
      h=dx*norm(x)+1.e-8 \n

pdata:\n
precomputed data for the function/Jacobian-vector/Preconditioner-vector
products.  Things will go better if you use this rather than hide the data
in global variables within the module for your function/Jacobian

If you use pdata in any of F!, Jvec, or Pvec, you must use in in all of them.
precomputed data for the function/Jacobian.
Things will go better if you use this rather than hide the data
in global variables within the module for your function/Jacobian. \n

printerr: default = true\n
I print a helpful message when the solver fails. To suppress that
message set printerr to false. \n

keepsolhist: default = false\n
Set this to true to get the history of the iteration in the output
tuple. This is on by default for scalar equations and off for systems.
Only turn it on if you have use for the data, which can get REALLY LARGE.\n

Output:\n
A named tuple (solution, functionval, history, stats, idid,
               errcode, solhist)
where

solution = converged result
functionval = F(solution)
history = the vector of residual norms (||F(x)||) for the iteration
stats = named tuple of the history of (ifun, ijac, ikfail), the number
of functions/jacobian-vector prodcuts/linear solver filures at each iteration.

I do not count the function values for a finite-difference derivative
because they count toward a Jacobian-vector product.

Linear solver failures need not cause the nonlinear iteration to fail.
You get a warning and that is all. \n

idid=true if the iteration succeeded and false if not. \n

errcode = 0 if if the iteration succeeded \n
        = -1 if the initial iterate satisfies the termination criteria
        = 10 if no convergence after maxit iterations \n

solhist:\n
This is the entire history of the iteration if you've set
keepsolhist=true\n
solhist is an N x K array where N is the length of x and K is the number
of iteration + 1. So, for scalar equations, it's a row vector.

### Example for ptcsol

#### The buckling beam problem.
You'll need to use TestProblems for
this to work. The preconditioner is a solver for the high order term.

```jldoctest
julia> using SIAMFANLEquations.TestProblems

julia> function PreCondBeam(v, x, bdata)
          J = bdata.D2
          ptv = J\v
       end
PreCondBeam (generic function with 1 method)

julia> n=63; maxit=1000; delta0 = 0.01; lambda = 20.0;

julia> bdata = beaminit(n, 0.0, lambda);

julia> x = bdata.x; u0 = x .* (1.0 .- x) .* (2.0 .- x); u0 .*= exp.(-10.0 * u0);

julia> FS = copy(u0); FPJV=zeros(n,20);

julia> pout = ptcsoli( FBeam!, u0, FS, FPJV; delta0 = delta0, pdata = bdata,
       eta = 1.e-2, rtol = 1.e-10, maxit = maxit, Pvec = PreCondBeam);

julia> # It takes a few iterations to get there.
       length(pout.history)
25

julia> [pout.history[1:5] pout.history[21:25]]
5×2 Array{Float64,2}:
 6.31230e+01  1.79578e+00
 7.45926e+00  2.65964e-01
 8.73598e+00  6.58278e-03
 2.91936e+01  8.35069e-06
 3.47969e+01  5.11594e-09

julia> # We get the nonnegative stedy state.
       norm(pout.solution,Inf)
2.19086e+00

n=63; maxit=1000; delta0 = 0.01; lambda = 20.0;

julia> # Use BiCGSTAB for the linear solver

julia> FS = copy(u0); FPJV=zeros(n,);

julia> pout = ptcsoli( FBeam!, u0, FS, FPJV; delta0 = delta0, pdata = bdata,
       eta = 1.e-2, rtol = 1.e-10, maxit = maxit,
       Pvec = PreCondBeam, lsolver="bicgstab");

julia> # Same number of iterations as GMRES, but each one costs double

julia> # the Jacobian-vector products and much less storage

julia> length(pout.history)
25

julia> [pout.history[1:5] pout.history[21:25]]
5×2 Matrix{Float64}:
 6.31230e+01  1.68032e+00
 7.47081e+00  2.35073e-01
 8.62095e+00  5.18262e-03
 2.96495e+01  3.23715e-06
 3.51504e+01  3.33107e-10

```

"""
function ptcsoli(
    F!,
    x0,
    FS,
    FPS,
    Jvec = dirder;
    rtol = 1.e-6,
    atol = 1.e-12,
    maxit = 20,
    lmaxit = -1,
    lsolver = "gmres",
    eta = 0.1,
    fixedeta = true,
    Pvec = nothing,
    PvecKnowsdelta = false,
    pside = "right",
    delta0 = 1.e-6,
    dx = 1.e-7,
    pdata = nothing,
    printerr = true,
    keepsolhist = false,
)
    itc = 0
    idid = true
    #
    #   Initialize the iteration
    #   As with the other codes, ItRules packages all the details of
    #   the problem so it's easy to pass them around.
    #
    (ItRules, x, n) = PTCKrylovinit(
        x0,
        dx,
        F!,
        Jvec,
        delta0,
        Pvec,
        PvecKnowsdelta,
        pside,
        lsolver,
        eta,
        fixedeta,
        lmaxit,
        maxit,
        printerr,
        pdata,
    )
    keepsolhist ? (solhist = solhistinit(n, maxit, x)) : (solhist = [])
    #
    # First Evaluation of the function. Initialize the iteration history.
    # Fix the tolerances for convergence and define the derivative FPF
    # outside of the main loop for scoping.
    #
    FS = EvalF!(F!, FS, x, pdata)
    resnorm = norm(FS)
    ItData = ItStatsPTCK(resnorm)
    tol = rtol * resnorm + atol
    etag = eta
    ke_report = false
    residratio = 1.0
    #
    # Preallocate a vector for the step
    #
    step = copy(x)
    #
    # If the initial iterate satisfies the termination criteria, tell me.
    #
    toosoon = (resnorm <= tol)
    #
    # The main loop stops on convergence or too many iterations.
    #
    delta = delta0
    while resnorm > tol && itc < maxit
        residm = resnorm
        newjac = 0
        newikfail = 0
        #
        # Comppute the Jacobian-vector product; update x, F(x), and delta.
        #
        etag = forcing(itc, residratio, etag, ItRules, tol, resnorm)
        (x, delta, FS, resnorm, Lstats) =
            PTCUpdatei(FPS, FS, x, ItRules, step, resnorm, delta, etag)
        resdiratio = resnorm / residm
        newjac = Lstats.lits
        linok = Lstats.idid
        linok || (ke_report = Krylov_Error(lmaxit, ke_report); newikfail = 1)
        #
        # Keep the books
        #
        updateStats!(ItData, resnorm, newjac, newikfail)
        itc += 1
        ~keepsolhist || (@views solhist[:, itc+1] .= x)
    end
    (idid, errcode) = PTCOK(resnorm, tol, toosoon, ItRules, printerr)
    itout = CloseIteration(x, FS, ItData, idid, errcode, keepsolhist, solhist)
    return itout
end


text:
function test_ln2_all()
    print(rpad("Math.Series.ln(2)...", 50, ' '))

    ln2_N2a = Brainstorm._Math._Series._Ln2._Machin.ln2_N2a
    ln2_N2b = Brainstorm._Math._Series._Ln2._Machin.ln2_N2b
    ln2_N2c = Brainstorm._Math._Series._Ln2._Machin.ln2_N2c

    ln2_N3a = Brainstorm._Math._Series._Ln2._Machin.ln2_N3a
    ln2_N3b = Brainstorm._Math._Series._Ln2._Machin.ln2_N3b
    ln2_N3c = Brainstorm._Math._Series._Ln2._Machin.ln2_N3c
    ln2_N3d = Brainstorm._Math._Series._Ln2._Machin.ln2_N3d
    ln2_N3e = Brainstorm._Math._Series._Ln2._Machin.ln2_N3e

    ln2_N4a = Brainstorm._Math._Series._Ln2._Machin.ln2_N4a
    ln2_N4b = Brainstorm._Math._Series._Ln2._Machin.ln2_N4b
    ln2_N4c = Brainstorm._Math._Series._Ln2._Machin.ln2_N4c

    ln2_N5a = Brainstorm._Math._Series._Ln2._Machin.ln2_N5a

    setprecision(400_000) do
        expected = [trunc(BigInt, log(big(2))*(big(10)^d)) for d in (1, 10, 100, 1_000, 10_000, 100_000)]

        @test all([ln2_N2a(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([ln2_N2b(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([ln2_N2c(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)

        @test all([ln2_N3a(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([ln2_N3b(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([ln2_N3c(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([ln2_N3d(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([ln2_N3e(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)

        @test all([ln2_N4a(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([ln2_N4b(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
        @test all([ln2_N4c(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)

        @test all([ln2_N5a(d) for d in (1, 10, 100, 1_000, 10_000, 100_000)] - expected .== 0)
    end

    println("PASS")
end


text:
function f = fVelW(x,Param)
switch lower(Param.ProfVelW)
  case 'linear'
    f=x(1);
  otherwise
    f=0;
end
end


text:
"""
    initialize_sbm_gwf_model(config::Config)

Initial part of the sbm_gwf model concept. The model contains:
    - the vertical SBM concept
    - the following lateral components:
        - 1-D kinematic wave for river flow
        - 1-D kinematic wave for overland flow
        - unconfined aquifer with groundwater flow in four directions (adjacent cells)

The unconfined aquifer contains a recharge, river and a drain (optional) boundary.

The initial part reads the input settings and data as defined in the Config object.
Will return a Model that is ready to run.
"""
function initialize_sbm_gwf_model(config::Config)

    # unpack the paths to the NetCDF files
    static_path = input_path(config, config.input.path_static)

    reader = prepare_reader(config)
    clock = Clock(config, reader)
    Δt = clock.Δt

    reinit = get(config.model, "reinit", true)::Bool
    do_reservoirs = get(config.model, "reservoirs", false)::Bool
    do_lakes = get(config.model, "lakes", false)::Bool
    do_drains = get(config.model, "drains", false)::Bool
    do_constanthead = get(config.model, "constanthead", false)::Bool

    kw_river_tstep = get(config.model, "kw_river_tstep", 0)
    kw_land_tstep = get(config.model, "kw_land_tstep", 0)
    kinwave_it = get(config.model, "kin_wave_iteration", false)::Bool

    nc = NCDataset(static_path)

    subcatch_2d =
        ncread(nc, config, "subcatchment"; optional = false, allow_missing = true)
    # indices based on catchment
    inds, rev_inds = active_indices(subcatch_2d, missing)
    n = length(inds)
    modelsize_2d = size(subcatch_2d)

    river_2d = ncread(
        nc,
        config,
        "river_location";
        optional = false,
        type = Bool,
        fill = false,
    )
    river = river_2d[inds]
    riverwidth_2d = ncread(
        nc,
        config,
        "lateral.river.width";
        optional = false,
        type = Float,
        fill = 0,
    )
    riverwidth = riverwidth_2d[inds]
    riverlength_2d = ncread(
        nc,
        config,
        "lateral.river.length";
        optional = false,
        type = Float,
        fill = 0,
    )
    riverlength = riverlength_2d[inds]

    altitude =
        ncread(nc, config, "altitude"; optional = false, sel = inds, type = Float)
    # read x, y coordinates and calculate cell length [m]
    y_nc = read_y_axis(nc)
    x_nc = read_x_axis(nc)
    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
    cellength = abs(mean(diff(x_nc)))

    sizeinmetres = get(config.model, "sizeinmetres", false)::Bool
    xl, yl = cell_lengths(y, cellength, sizeinmetres)
    riverfrac = river_fraction(river, riverlength, riverwidth, xl, yl)

    # initialize vertical SBM concept
    sbm = initialize_sbm(nc, config, riverfrac, inds)

    inds_riv, rev_inds_riv = active_indices(river_2d, 0)
    nriv = length(inds_riv)

    # reservoirs
    pits = zeros(Bool, modelsize_2d)
    if do_reservoirs
        reservoirs, resindex, reservoir, pits =
            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(Δt))
    else
        reservoir = ()
        reservoirs = nothing
        resindex = fill(0, nriv)
    end

    # lakes
    if do_lakes
        lakes, lakeindex, lake, pits =
            initialize_natural_lake(config, nc, inds_riv, nriv, pits, tosecond(Δt))
    else
        lake = ()
        lakes = nothing
        lakeindex = fill(0, nriv)
    end

    # overland flow (kinematic wave)
    βₗ = ncread(
        nc,
        config,
        "lateral.land.slope";
        optional = false,
        sel = inds,
        type = Float,
    )
    clamp!(βₗ, 0.00001, Inf)
    ldd_2d = ncread(nc, config, "ldd"; optional = false, allow_missing = true)

    ldd = ldd_2d[inds]

    dl = map(detdrainlength, ldd, xl, yl)
    dw = (xl .* yl) ./ dl
    sw = map(det_surfacewidth, dw, riverwidth, river)

    olf = initialize_surfaceflow_land(
        nc,
        config,
        inds;
        sl = βₗ,
        dl = dl,
        width = sw,
        wb_pit = pits[inds],
        iterate = kinwave_it,
        tstep = kw_land_tstep,
        Δt = Δt,
    )

    graph = flowgraph(ldd, inds, pcr_dir)

    # river flow (kinematic wave)
    riverlength = riverlength_2d[inds_riv]
    riverwidth = riverwidth_2d[inds_riv]
    minimum(riverlength) > 0 || error("river length must be positive on river cells")
    minimum(riverwidth) > 0 || error("river width must be positive on river cells")

    ldd_riv = ldd_2d[inds_riv]
    graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)

    # the indices of the river cells in the land(+river) cell vector
    index_river = filter(i -> !isequal(river[i], 0), 1:n)
    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)

    rf = initialize_surfaceflow_river(
        nc,
        config,
        inds_riv;
        dl = riverlength,
        width = riverwidth,
        wb_pit = pits[inds_riv],
        reservoir_index = resindex,
        reservoir = reservoirs,
        lake_index = lakeindex,
        lake = lakes,
        river = river,
        iterate = kinwave_it,
        tstep = kw_river_tstep,
        Δt = Δt,
    )

    # unconfined aquifer
    if do_constanthead
        constanthead = ncread(
            nc,
            config,
            "lateral.subsurface.constant_head";
            sel = inds,
            type = Float,
            fill = mv,
        )
        index_constanthead = filter(i -> !isequal(constanthead[i], mv), 1:n)
        constant_head = ConstantHead(constanthead[index_constanthead], index_constanthead)
    else
        constant_head = ConstantHead{Float}(Float[], Int64[])
    end

    conductivity = ncread(
        nc,
        config,
        "lateral.subsurface.conductivity";
        sel = inds,
        type = Float,
    )
    specific_yield = ncread(
        nc,
        config,
        "lateral.subsurface.specific_yield";
        sel = inds,
        type = Float,
    )

    connectivity = Connectivity(inds, rev_inds, xl, yl)
    initial_head = altitude .- Float(0.10) # cold state for groundwater head
    initial_head[index_river] = altitude[index_river]

    if do_constanthead
        initial_head[constant_head.index] = constant_head.head
    end

    aquifer = UnconfinedAquifer(
        initial_head,
        conductivity,
        altitude,
        altitude .- sbm.soilthickness ./ Float(1000.0),
        xl .* yl,
        specific_yield,
        zeros(Float, connectivity.nconnection),  # conductance
    )

    # reset zi and satwaterdepth with groundwater head from unconfined aquifer
    sbm.zi .= (altitude .- initial_head) .* 1000.0
    sbm.satwaterdepth .= (sbm.soilthickness .- sbm.zi) .* (sbm.θₛ .- sbm.θᵣ)

    # river boundary of unconfined aquifer
    infiltration_conductance = ncread(
        nc,
        config,
        "lateral.subsurface.infiltration_conductance";
        sel = inds_riv,
        type = Float,
    )
    exfiltration_conductance = ncread(
        nc,
        config,
        "lateral.subsurface.exfiltration_conductance";
        sel = inds_riv,
        type = Float,
    )
    river_bottom = ncread(
        nc,
        config,
        "lateral.subsurface.river_bottom";
        sel = inds_riv,
        type = Float,
    )

    river_flux = fill(mv, nriv)
    river_stage = fill(mv, nriv)
    river = River(
        river_stage,
        infiltration_conductance,
        exfiltration_conductance,
        river_bottom,
        river_flux,
        index_river,
    )

    # recharge boundary of unconfined aquifer
    r = fill(mv, n)
    recharge = Recharge(r, zeros(Float, n), collect(1:n))

    # drain boundary of unconfined aquifer (optional)
    if do_drains
        drain_2d =
            ncread(nc, config, "lateral.subsurface.drain"; type = Bool, fill = false)

        drain = drain_2d[inds]
        # check if drain occurs where overland flow is not possible (sw = 0.0)
        # and correct if this is the case
        false_drain = filter(i -> !isequal(drain[i], 0) && sw[i] == Float(0), 1:n)
        n_false_drain = length(false_drain)
        if n_false_drain > 0
            drain_2d[inds[false_drain]] .= 0
            drain[false_drain] .= 0
            @info "$n_false_drain drain locations are removed that occur where overland flow
             is not possible (overland flow width is zero)"
        end
        inds_drain, rev_inds_drain = active_indices(drain_2d, 0)

        drain_elevation = ncread(
            nc,
            config,
            "lateral.subsurface.drain_elevation";
            sel = inds,
            type = Float,
            fill = mv,
        )
        drain_conductance = ncread(
            nc,
            config,
            "lateral.subsurface.drain_conductance";
            sel = inds,
            type = Float,
            fill = mv,
        )
        index_drain = filter(i -> !isequal(drain[i], 0), 1:n)
        drain_flux = fill(mv, length(index_drain))
        drains = Drainage(
            drain_elevation[index_drain],
            drain_conductance[index_drain],
            drain_flux,
            index_drain,
        )
        drain = (indices = inds_drain, reverse_indices = rev_inds_drain)
        aquifer_boundaries = AquiferBoundaryCondition[recharge, river, drains]
    else
        aquifer_boundaries = AquiferBoundaryCondition[recharge, river]
        drain = ()
    end

    gwf = GroundwaterFlow(aquifer, connectivity, constant_head, aquifer_boundaries)

    # map GroundwaterFlow and its boundaries
    if do_drains
        subsurface_map = (
            flow = gwf,
            recharge = gwf.boundaries[1],
            river = gwf.boundaries[2],
            drain = gwf.boundaries[3],
        )
    else
        subsurface_map =
            (flow = gwf, recharge = gwf.boundaries[1], river = gwf.boundaries[2])
    end

    # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
    # subdomain is equal to the complete domain
    toposort = topological_sort_by_dfs(graph)
    toposort_riv = topological_sort_by_dfs(graph_riv)
    index_pit_land = findall(x -> x == 5, ldd)
    index_pit_river = findall(x -> x == 5, ldd_riv)
    subbas_order, indices_subbas, topo_subbas =
        kinwave_set_subdomains(config, graph, toposort, index_pit_land)
    subriv_order, indices_subriv, topo_subriv =
        kinwave_set_subdomains(config, graph_riv, toposort_riv, index_pit_river)

    modelmap =
        (vertical = sbm, lateral = (subsurface = subsurface_map, land = olf, river = rf))
    indices_reverse = (
        land = rev_inds,
        river = rev_inds_riv,
        reservoir = isempty(reservoir) ? nothing : reservoir.reverse_indices,
        lake = isempty(lake) ? nothing : lake.reverse_indices,
        drain = isempty(drain) ? nothing : rev_inds_drain,
    )
    writer = prepare_writer(
        config,
        modelmap,
        indices_reverse,
        x_nc,
        y_nc,
        nc,
        extra_dim = (name = "layer", value = Float64.(1:sbm.maxlayers)),
    )
    close(nc)

    # for each domain save:
    # - the directed acyclic graph (graph),
    # - the traversion order (order),
    # - upstream_nodes,
    # - subdomains for the kinematic wave domains for parallel execution (execution order of
    #   subbasins (subdomain_order), traversion order per subbasin (topo_subdomain) and
    #   Vector indices per subbasin matching the traversion order of the complete domain
    #   (indices_subdomain))
    # - the indices that map it back to the two dimensional grid (indices)

    # for the land domain the x and y length [m] of the grid cells are stored
    # for reservoirs and lakes indices information is available from the initialization
    # functions
    land = (
        graph = graph,
        upstream_nodes = filter_upsteam_nodes(graph, olf.wb_pit),
        subdomain_order = subbas_order,
        topo_subdomain = topo_subbas,
        indices_subdomain = indices_subbas,
        order = toposort,
        indices = inds,
        reverse_indices = rev_inds,
        xl = xl,
        yl = yl,
        altitude = altitude,
    )
    river = (
        graph = graph_riv,
        upstream_nodes = filter_upsteam_nodes(graph_riv, rf.wb_pit),
        subdomain_order = subriv_order,
        topo_subdomain = topo_subriv,
        indices_subdomain = indices_subriv,
        order = toposort_riv,
        indices = inds_riv,
        reverse_indices = rev_inds_riv,
    )

    model = Model(
        config,
        (; land, river, reservoir, lake, drain, index_river, frac_toriver),
        (subsurface = subsurface_map, land = olf, river = rf),
        sbm,
        clock,
        reader,
        writer,
        SbmGwfModel(),
    )

    # read and set states in model object if reinit=false
    if reinit == false
        instate_path = input_path(config, config.state.path_input)
        @info "Set initial conditions from state file `$instate_path`."
        state_ncnames = ncnames(config.state)
        set_states(instate_path, model, state_ncnames, type = Float, dimname = :layer)
        # update kinematic wave volume for river and land domain
        @unpack lateral = model
        # makes sure land cells with zero flow width are set to zero q and h
        for i in eachindex(lateral.land.width)
            if lateral.land.width[i] <= 0.0
                lateral.land.q[i] = 0.0
                lateral.land.h[i] = 0.0
            end
        end
        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl

        if do_lakes
            # storage must be re-initialized after loading the state with the current
            # waterlevel otherwise the storage will be based on the initial water level
            lakes.storage .=
                initialize_storage(lakes.storfunc, lakes.area, lakes.waterlevel, lakes.sh)
        end
    else
        @info "Set initial conditions from default values."
    end

    return model
end

"update the sbm_gwf model for a single timestep"
function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
    @unpack lateral, vertical, network, clock, config = model

    inds_riv = network.index_river

    # extract water levels h_av [m] from the land and river domains
    # this is used to limit open water evaporation
    vertical.waterlevel_land .= lateral.land.h_av .* 1000.0
    vertical.waterlevel_river[inds_riv] .= lateral.river.h_av .* 1000.0

    # vertical sbm concept is updated until snow state, after that (optional)
    # snow transport is possible
    update_until_snow(vertical, config)

    # lateral snow transport
    if get(config.model, "masswasting", false)::Bool
        lateral_snow_transport!(
            vertical.snow,
            vertical.snowwater,
            lateral.land.sl,
            network.land,
        )
    end

    # update vertical sbm concept until recharge [mm]
    update_until_recharge(vertical, config)

    # set river stage (groundwater) to average h from kinematic wave
    lateral.subsurface.river.stage .= lateral.river.h_av .+ lateral.subsurface.river.bottom

    # determine stable time step for groundwater flow
    Δt_gw = stable_timestep(lateral.subsurface.flow.aquifer) # time step in day (Float64)
    Δt_sbm = (vertical.Δt / tosecond(basetimestep)) # vertical.Δt is in seconds (Float64)
    if Δt_gw < Δt_sbm
        @warn(
            "stable time step Δt $Δt_gw for groundwater flow is smaller than sbm Δt $Δt_sbm"
        )
    end

    Q = zeros(vertical.n)
    # exchange of recharge between vertical sbm concept and groundwater flow domain
    # recharge rate groundwater is required in units [m d⁻¹]
    lateral.subsurface.recharge.rate .= vertical.recharge ./ 1000.0 .* (1.0 / Δt_sbm)
    # update groundwater domain
    update(lateral.subsurface.flow, Q, Δt_sbm)

    # determine excess water depth [m] (exfiltwater) in groundwater domain (head > surface)
    # and reset head
    exfiltwater =
        (
            lateral.subsurface.flow.aquifer.head .-
            min.(lateral.subsurface.flow.aquifer.head, lateral.subsurface.flow.aquifer.top)
        ) .* storativity(lateral.subsurface.flow.aquifer)
    lateral.subsurface.flow.aquifer.head .=
        min.(lateral.subsurface.flow.aquifer.head, lateral.subsurface.flow.aquifer.top)

    # update vertical sbm concept (runoff, ustorelayerdepth and satwaterdepth)
    update_after_subsurfaceflow(
        vertical,
        (network.land.altitude .- lateral.subsurface.flow.aquifer.head) .* 1000.0, # zi [mm] in vertical concept SBM
        exfiltwater .* 1000.0,
    )

    ssf_toriver = zeros(vertical.n)
    ssf_toriver[inds_riv] = -lateral.subsurface.river.flux ./ lateral.river.Δt
    surface_routing(model, ssf_toriver = ssf_toriver)

    write_output(model)

    # update the clock
    advance!(clock)

    return model
end


text:
function complete_lower_bound_grad(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    #
    # ∫ q(α) q(β) q(f) log𝓝(y| Af + b, S) + log𝓝(f|0,K) + log𝓝(a|0,v²ₐI) + log𝓝(b|0,v²ᵦ) dα dβ df
    #                           (1)               (2)             (3)             (4)

    gradα, gradΣα, gradβ, gradΣβ = exact_expectation_inside_log_of_lower_bound_grad(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    #
    # H[q(α)]  +  H[q(β)]  +  H[q(f)]
    #

    # remember we parametrise the n-th entry of the diagonal covariance
    # Σa with elements exp(pₙ), i.e. gradient calculated below is dependent
    # on this parametrisation (see "unpack" function)
    gradΣα += 0.5 * ones(size(Σα, 1))

    # same comment as above applies here
    gradΣβ += 0.5 * ones(size(Σβ, 1))

    return gradα, gradΣα, gradβ, gradΣβ

end


text:
# ************************************************************************
# FAUST Architecture File
# Copyright (C) 2021 GRAME, Centre National de Creation Musicale
# ---------------------------------------------------------------------

# This is sample code. This file is provided as an example of minimal
# FAUST architecture file. Redistribution and use in source and binary
# forms, with or without modification, in part or in full are permitted.
# In particular you can create a derived work of this FAUST architecture
# and distribute that work under terms of your choice.

# This sample code is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# ************************************************************************

const FAUSTFLOAT = Float32

# Generated code
<<includeIntrinsic>>
<<includeclass>>

# Retrieve the application name
mutable struct NameMeta <: FMeta
    name::String
end

import Faust.declare!

function declare!(m::NameMeta, key::String, value::String)
    if (key == "name")
        m.name = value;
    end
end

# Testing
using Plots

samplerate = Int32(44100)
block_size = Int32(512 * 10)

test!() = begin
    # Init DSP
    my_dsp = mydsp{REAL}()
    init!(my_dsp, samplerate)

    m = NameMeta("")
    metadata!(my_dsp, m)
    println("Application name: ", m.name, "\n")

    println("getNumInputs: ", getNumInputs(my_dsp))
    println("getNumOutputs: ", getNumOutputs(my_dsp), "\n")

    # Create a MapUI controller
    map_ui = MapUI(my_dsp)
    buildUserInterface!(my_dsp, map_ui)

    # Print all zones
    println("Path/UIZone dictionary: ", getZoneMap(map_ui), "\n")

    inputs = zeros(FAUSTFLOAT, block_size, getNumInputs(my_dsp))
    outputs = zeros(FAUSTFLOAT, block_size, getNumOutputs(my_dsp))
    compute!(my_dsp, block_size, inputs, outputs)

    # display the outputs
    display(plot(outputs, layout = (getNumOutputs(my_dsp), 1)))

end

test!()


text:
module Err
export PatternCompilationError, InternalException, SyntaxError, UnknownExtension, PatternUnsolvedException

PatternUnsolvedException = ErrorException

struct PatternCompilationError <: Exception
    line::Union{LineNumberNode,Nothing}
    msg::AbstractString
end

struct InternalException <: Exception
    msg :: String
end

struct SyntaxError <: Exception
    msg :: String
end

struct UnknownExtension <: Exception
    ext :: Union{String, Symbol}
end

end


text:
struct Occupancy{V} <: Measurement
    val::V
end

Base.rand(::Type{Occupancy}, x::Int, y::Int) = Base.rand(Occupancy, 0.5, x, y)
Base.rand(::Type{Occupancy}, p::Float64, x::Int, y::Int) = begin
    z = zeros(Occupancy, x, y)
    for i in eachindex(z)
        z[i] = rand() < p
    end
    return z
end

Float64(t::Occupancy{T}) where {T} = Float64(t.val)

Base.zero(::Occupancy{T}) where {T} = Occupancy{T}(zero(T))
Base.zero(::Type{Occupancy}) = Occupancy{Bool}(false)
Base.zero(::Type{<:Occupancy{T}}) where {T} = Occupancy(zero(T))

Base.show(io::IO, t::Occupancy{T}) where {T} = Base.show(io, T(t.val))

Base.convert(::Type{Occupancy}, v::V) where {V<:Bool} = Occupancy{V}(v)

Base.isless(a::Occupancy, b::Occupancy) = isless(a.val, b.val)
Base.oneunit(::Type{<:Occupancy{T}}) where {T} = Occupancy(one(T))
Base.one(::Type{<:Occupancy{T}}) where {T} = Occupancy(one(T))
Base.one(::Occupancy{T}) where {T} = Occupancy{T}(one(T))

Base.:*(x::Occupancy, v::Number) = Occupancy(x.val * v)
Base.:*(v::Number, x::Occupancy) = Occupancy(x.val * v)
Base.:/(x::Occupancy, v::Number) = Occupancy(x.val / v)
Base.:+(x1::Occupancy, x2::Occupancy) = Occupancy(x1.val & x2.val)
Base.:-(x1::Occupancy, x2::Occupancy) = Occupancy(x1.val & x2.val)
Base.:-(o::Occupancy{T}, n::Number) where {T} = Occupancy{T}(zero(T))


text:
@discord_object struct Invite
    code::String
    guild::Guild
    channel::DiscordChannel
    inviter::User
    target_user::User
    target_user_type::TargetUserType.TargetUserTypeEnum
    approximate_presence_count::Int
    approximate_member_count::Int
    # https://discord.com/developers/docs/resources/invite#invite-metadata-object
    uses::Int
    max_uses::Int
    max_age::Int
    temporary::Bool
    created_at::Union{String, DateTime}
end


text:
# Autogenerated wrapper script for ReadStat_jll for x86_64-apple-darwin
export libreadstat, readstat

using Libiconv_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("ReadStat")
JLLWrappers.@declare_library_product(libreadstat, "@rpath/libreadstat.1.dylib")
JLLWrappers.@declare_executable_product(readstat)
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libreadstat,
        "lib/libreadstat.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        readstat,
        "bin/readstat",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export
    Invertible

"""
    struct Invertible{I,O} <: SFunc{Tuple{I},O,Nothing}

An invertible sfunc, with both a `forward` and a `inverse` function provided.

# Additional supported operators
- `support`
- `support_quality`
- `sample`
- `logcpdf`
- `bounded_probs`
- `make_factors`
- `compute_pi`
- `send_lambda`
"""
struct Invertible{I,O} <: SFunc{Tuple{I},O}
    forward :: Function
    inverse :: Function
end

@impl begin
    struct InvertibleSupport end
    function support(sf::Invertible{I,O}, parranges, size, curr) where {I,O}
        invsup = map(sf.forward, parranges[1])
        return invsup
    end
end

@impl begin
    struct InvertibleSupportQuality end

    function support_quality(::Invertible{I,O}, parranges) where {I,O}
        return :CompleteSupport
    end
end

@impl begin
    struct InvertibleSample end
    function sample(sf::Invertible{I,O}, i::Tuple{I}) where {I,O}
        return sf.forward(i[1])
    end
end

@impl begin
    struct InvertibleLogcpdf end
    function logcpdf(sf::Invertible{I,O}, i::Tuple{I}, o::O) where {I,O}
        return sf.forward(i[1]) == o ? 0.0 : -Inf
    end
end

@impl begin
    struct InvertibleBoundedProbs end

    function bounded_probs(sf::Invertible{I,O}, range::Vector{<:O},
            parranges) where {I,O}
        result = Float64[]
        for i in parranges[1]
            for o in range
                p = sf.forward(i) == o ? 1.0 : 0.0
                push!(result, p)
            end
        end
        return (result, result)
    end
end

@impl begin
    struct InvertibleMakeFactors end

    function make_factors(sf::Invertible{I,O}, range::Vector{<:O}, parranges,
            id, parids) where {I,O}
        dims = (length(parranges[1]), length(range))
        keys = (parids[1], id)
        entries = bounded_probs(sf, range, parranges)[1]
        facts = [Factor(dims, keys, entries)]
        return (facts, facts)
    end
end

@impl begin
    struct InvertibleComputePi end

    function compute_pi(sf::Invertible{I,O}, range::Vector{<:O}, parranges,
            incoming_pis) where {I,O}
        ps = [cpdf(incoming_pis[1], (), sf.inverse(o)) for o in range]
        return Cat(range, ps)
    end
end

@impl begin
    struct InvertibleSendLambda end

    function send_lambda(sf::Invertible{I,O}, lambda, range, parranges, incoming_pis,
            parent_idx) where {I,O}
        @assert parent_idx == 1
        ls = [get_log_score(lambda, sf.forward(x)) for x in parranges[1]]
        return LogScore(parranges[1], ls)
    end
end
#=
@impl begin
    struct InvertibleInitialStats end
    initial_stats(::Invertible) = nothing
end

@impl begin
    struct InvertibleAccumulateStats end
    accumulate_stats(::Invertible, existing_stats, new_stats) = nothing
end

@impl begin
    struct InvertibleExpectedStats end
    expected_stats(::Invertible, range, parranges, incoming_pis, child_lambda) =
        nothing
end

@impl begin
    struct InvertibleMaximizeStats end
    maximize_stats(::Invertible, stats) = nothing
end
=#


text:
"""
"""
function lazy_append(a::AbstractArray,b::AbstractArray)
  AppendedArray(a,b)
end

"""
"""
function lazy_split(f::AbstractArray,n::Integer)
  _lazy_split(f,n)
end

function lazy_split(f::Fill,n::Integer)
  l = length(f)
  @assert n <= l
  a = Fill(f.value,n)
  b = Fill(f.value,l-n)
  a,b
end

function lazy_split(f::CompressedArray,n::Integer)
  _a , _b =_lazy_split(f,n)
  a = _compact_values_ptrs(_a)
  b = _compact_values_ptrs(_b)
  a,b
end

function _lazy_split(f,n)
  l = length(f)
  @assert n <= l
  ids_a = collect(1:n)
  ids_b = collect((n+1):l)
  a = lazy_map(Reindex(f),ids_a)
  b = lazy_map(Reindex(f),ids_b)
  a,b
end

function _compact_values_ptrs(a)
  v = a.values
  p = a.ptrs
  touched = fill(false,length(v))
  for i in p
    touched[i] = true
  end
  if all(touched) || length(p) == 0
    return a
  else
    j_to_i = findall(touched)
    i_to_j = zeros(Int,length(v))
    i_to_j[j_to_i] = 1:length(j_to_i)
    _v = v[j_to_i]
    _p = i_to_j[p]
    return CompressedArray(_v, _p)
  end
end

struct AppendedArray{T,A,B} <: AbstractVector{T}
  a::A
  b::B
  function AppendedArray(a::AbstractArray,b::AbstractArray)
    A = typeof(a)
    B = typeof(b)
    ai = testitem(a)
    bi = testitem(b)
    T = eltype(collect((ai,bi)))
    new{T,A,B}(a,b)
  end
end

Base.IndexStyle(::Type{<:AppendedArray}) = IndexLinear()

function Base.getindex(v::AppendedArray,i::Integer)
  l = length(v.a)
  if i > l
    v.b[i-l]
  else
    v.a[i]
  end
end

Base.size(v::AppendedArray) = (length(v.a)+length(v.b),)

function array_cache(v::AppendedArray)
  ca = array_cache(v.a)
  cb = array_cache(v.b)
  (ca,cb)
end

function getindex!(cache,v::AppendedArray,i::Integer)
  ca, cb = cache
  l = length(v.a)
  if i > l
    getindex!(cb,v.b,i-l)
  else
    getindex!(ca,v.a,i)
  end
end

function lazy_map(::typeof(evaluate),::Type{T},b::Fill,a::AppendedArray...) where T
  f = b.value
  la = map(ai->length(ai.a),a)
  lb = map(ai->length(ai.b),a)
  if all(la .== first(la)) && all(lb .== first(lb))
    c_a = lazy_map(f,map(ai->ai.a,a)...)
    c_b = lazy_map(f,map(ai->ai.b,a)...)
    lazy_append(c_a,c_b)
  else
    s = _common_size(a...)
    LazyArray(T,Fill(f,s...),a...)
  end
end

# In that case we don't implement it for ::Type{T} variant since we want to
# avoid to run type inference since the first entry in a can have non-concrete
# eltype
function lazy_map(k::typeof(evaluate),a::AppendedArray...)
  @check all(map(i->length(i)==length(a[1]),a))
  la = map(ai->length(ai.a),a)
  if all(la .== first(la))
    c_a = lazy_map(k,map(ai->ai.a,a)...)
    c_b = lazy_map(k,map(ai->ai.b,a)...)
    lazy_append(c_a,c_b)
  else
    ai = map(testitem,a)
    T = return_type(k, ai...)
    lazy_map(k,T,a...)
  end
end

# In that case we don't implement it for ::Type{T} variant since we want to
# avoid to run type inference since a can have non-concrete
# eltype
function lazy_map(k::typeof(evaluate),a::AppendedArray,b::AbstractArray)
  @assert length(a) == length(b)
  n = length(a.a)
  c = lazy_append(lazy_split(b,n)...)
  lazy_map(evaluate,a,c)
end

# Very important optimization to compute error norms efficiently e.g. in EmbeddedFEM
function Base.sum(a::AppendedArray)
  sum(a.a) + sum(a.b)
end


text:
# Julia internal methods overload.
import Base: ==
==(pb1::ParamBox, pb2::ParamBox) = (pb1[] == pb2[] &&
                                    pb1.canDiff[] == pb2.canDiff[] &&
                                    pb1.index[] == pb2.index[] &&
                                    typeof(pb1.map) === typeof(pb2.map))

import Base: show
const nSigShown = 10
function show(io::IO, pb::ParamBox)
    c = pb.canDiff[] ? :green : :light_black
    print(io, typeof(pb))
    print(io, "(", round(pb.data[], sigdigits=nSigShown), ")")
    print(io, "[")
    printstyled(io, "∂", color=c)
    print(io, "][")
    printstyled(io, "$(pb|>getVar)", color=:cyan)
    print(io, "]")
end

function show(io::IO, gf::GaussFunc)
    print(io, typeof(gf), "(xpn=")
    show(io, gf.xpn)
    print(io, ", con=")
    show(io, gf.con)
    print(io, ")")
end

function show(io::IO, bf::BasisFunc)
    xyz1 = bf.ijk[1]
    xyz2 = ""
    print(io, typeof(bf))
    print(io, "(gauss, subshell, center)[")
    printstyled(io, xyz1, color=:cyan)
    print(io, xyz2)
    print(io, "][", round(bf.center[1](), sigdigits=nSigShown), ", ",
                    round(bf.center[2](), sigdigits=nSigShown), ", ",
                    round(bf.center[3](), sigdigits=nSigShown), "]")
end

function show(io::IO, bf::BasisFuncs)
    OON = typeof(bf).parameters[3]
    SON = SubshellDimList[bf.subshell]
    if OON == 1
        xyz1 = bf.ijk[1]
        xyz2 = ""
    else
        xyz1 = "$(bf.ijk |> length)"
        xyz2 = "/$(SON)"
    end
    print(io, typeof(bf))
    print(io, "(gauss, subshell, center)[")
    printstyled(io, xyz1, color=:cyan)
    print(io, xyz2)
    print(io, "][", round(bf.center[1](), sigdigits=nSigShown), ", ",
                    round(bf.center[2](), sigdigits=nSigShown), ", ",
                    round(bf.center[3](), sigdigits=nSigShown), "]")
end

function show(io::IO, bfm::BasisFuncMix)
    print(io, typeof(bfm))
    print(io, "(BasisFunc, param)")
end

function show(io::IO, gtb::GTBasis)
    print(io, typeof(gtb))
    print(io, "(basis, S, Te, eeI, getVne, getHcore)")
end

function show(io::IO, box::GridBox)
    print(io, typeof(box))
    print(io, "(num, len, coord)")
end

function show(io::IO, vars::HFtempVars)
    print(io, typeof(vars))
    print(io, "(shared.Etots=[", round(vars.shared.Etots[1], sigdigits=nSigShown),", … , ",
                                 round(vars.shared.Etots[end], sigdigits=nSigShown), "], "*
              "shared.Dtots, Cs, Es, Ds, Fs)")
end

function show(io::IO, vars::HFfinalVars)
    print(io, typeof(vars))
    print(io, "(E0HF=", round(vars.E0HF, sigdigits=nSigShown), ", C, F, D, Emo, occu, temp"*
              ", isConverged)")
end

import Base: +
+(bfm1::CompositeGTBasisFuncs{<:Any, 1}, bfm2::CompositeGTBasisFuncs{<:Any, 1}) =
add(bfm1, bfm2)

import Base: *
*(bfm1::GaussFunc, bfm2::GaussFunc) = mul(bfm1, bfm2)

*(bfm::GaussFunc, coeff::Real) = mul(bfm, coeff)

*(coeff::Real, bfm::GaussFunc) = mul(coeff, bfm)

*(bfm1::CompositeGTBasisFuncs{<:Any, 1}, bfm2::CompositeGTBasisFuncs{<:Any, 1}) =
mul(bfm1, bfm2)

*(bfm::CompositeGTBasisFuncs{<:Any, 1}, coeff::Real) =
mul(bfm, coeff)

*(coeff::Real, bfm::CompositeGTBasisFuncs{<:Any, 1}) =
mul(coeff, bfm)

# Iteration Interface
import Base: iterate, size, length, ndims
iterate(pb::ParamBox) = (pb.data[], nothing)
iterate(::ParamBox, _) = nothing
size(::ParamBox) = ()
length(::ParamBox) = 1
ndims(::ParamBox) = 0
size(::ParamBox, d::Integer) = d == 1 ? 1 : throw(BoundsError())

iterate(gf::GaussFunc) = (gf, nothing)
iterate(::GaussFunc, _) = nothing
size(::GaussFunc) = ()
length(::GaussFunc) = 1

iterate(bf::BasisFunc) = (bf, nothing)
iterate(::BasisFunc, _) = nothing
size(::BasisFunc) = ()
length(::BasisFunc) = 1

iterate(bfm::BasisFuncMix) = (bfm, nothing)
iterate(::BasisFuncMix, _) = nothing
size(::BasisFuncMix) = ()
length(::BasisFuncMix) = 1

function iterate(bfs::CompositeGTBasisFuncs{<:Any, N}) where {N}
    item, state = iterate(bfs.ijk)
    (BasisFunc(bfs.center, bfs.gauss, ijkOrbitalList[item], bfs.normalizeGTO), state)
end
function iterate(bfs::CompositeGTBasisFuncs{<:Any, N}, state) where {N}
    iter = iterate(bfs.ijk, state)
    iter !== nothing ? (BasisFunc(bfs.center, bfs.gauss, ijkOrbitalList[iter[1]],
                                  bfs.normalizeGTO), iter[2]) : nothing
end
size(::CompositeGTBasisFuncs{<:Any, N}) where {N} = (N,)
length(::CompositeGTBasisFuncs{<:Any, N}) where {N} = N

size(x::SpatialOrbital, d::Integer) = d == 1 ? length(x) : throw(BoundsError())

# Indexing Interface
import Base: getindex, setindex!, firstindex, lastindex, eachindex, axes
getindex(pb::ParamBox) = pb.data[]
getindex(pb::ParamBox, ::Val{:first}) = getindex(pb)
getindex(pb::ParamBox, ::Val{:last}) = getindex(pb)
setindex!(pb::ParamBox, d) = begin pb.data[] = d end
firstindex(::ParamBox) = Val(:first)
lastindex(::ParamBox) = Val(:last)
axes(::ParamBox) = ()

getindex(gf::GaussFunc) = gf.param |> collect
getindex(gf::GaussFunc, ::Val{:first}) = getindex(gf)
getindex(gf::GaussFunc, ::Val{:last}) = getindex(gf)
firstindex(::GaussFunc) = Val(:first)
lastindex(::GaussFunc) = Val(:last)

getindex(bf::BasisFunc) = bf.gauss |> collect
getindex(bf::BasisFunc, ::Val{:first}) = getindex(bf)
getindex(bf::BasisFunc, ::Val{:last}) = getindex(bf)
firstindex(::BasisFunc) = Val(:first)
lastindex(::BasisFunc) = Val(:last)

getindex(bfm::BasisFuncMix) = (collect ∘ flatten)( getfield.(bfm.BasisFunc, :gauss) )
getindex(bfm::BasisFuncMix, ::Val{:first}) = getindex(bfm)
getindex(bfm::BasisFuncMix, ::Val{:last}) = getindex(bfm)
firstindex(::BasisFuncMix) = Val(:first)
lastindex(::BasisFuncMix) = Val(:last)

getindex(bfs::BasisFuncs, i) =
BasisFunc(bfs.center, bfs.gauss, ijkOrbitalList[bfs.ijk[i]], bfs.normalizeGTO)
getindex(bfs::BasisFuncs{<:Any, <:Any, N}, ::Colon) where {N} = [getindex(bfs, i) for i=1:N]
firstindex(bfs::BasisFuncs) = 1
lastindex(::BasisFuncs{<:Any, <:Any, N}) where {N} = N
eachindex(bfs::BasisFuncs{<:Any, <:Any, N}) where {N} = Base.OneTo(lastindex(bfs))
getindex(bfs::BasisFuncs) = getfield.(bfs[:], :gauss) |> flatten

# Broadcasting Interface
import Base: broadcastable
broadcastable(pb::ParamBox) = Ref(pb)
broadcastable(gf::GaussFunc) = Ref(gf)
broadcastable(bf::CompositeGTBasisFuncs{<:Any, 1}) = Ref(bf)
broadcastable(bfs::BasisFuncs) = getindex(bfs, :)

# Quiqbox methods overload.
## Method overload of `hasBoolRelation` from Tools.jl.
function hasBoolRelation(boolFunc::F, pb1::ParamBox, pb2::ParamBox;
                         ignoreFunction::Bool=false,
                         ignoreContainer::Bool=false,
                         decomposeNumberCollection::Bool=false) where {F<:Function}
    if ignoreContainer
        # boolFunc(pb1.data, pb2.data) &&
        boolFunc(pb1(), pb2())
    elseif boolFunc(pb1.canDiff[], pb2.canDiff[])
        if ignoreFunction
            return boolFunc(pb1.data, pb2.data)
        else
            return (boolFunc(pb1.map, pb2.map) && boolFunc(pb1.data, pb2.data))
        end
    else
        false
    end
end


text:
const ROADWAY = gen_straight_roadway(1, 20.0)

function create_vehicle(x::Float64, y::Float64, θ::Float64 = 0.0; id::Int64=1)
    s = VehicleState(VecSE2(x, y, θ), ROADWAY, 0.0)
    return Entity(s, VehicleDef(), id)
end

const VEH_REF = create_vehicle(0.0, 0.0, 0.0, id=1)

function no_collisions(xspace, yspace, thetaspace)
    for (x,y,th) in zip(xspace, yspace, thetaspace)
        veh2 = create_vehicle(x,y,th,id=2)
        @assert !collision_checker(VEH_REF, veh2) "Error for veh ", veh2
    end
    return true
end

function collisions(xspace, yspace, thetaspace)
    for (x,y,th) in zip(xspace, yspace, thetaspace)
        veh2 = create_vehicle(x,y,th,id=2)
        @assert collision_checker(VEH_REF, veh2) "Error for veh ", veh2
    end
    return true
end

@testset "Minkowski" begin
    @test AutomotiveSimulator.cyclic_shift_left!([1,2,3,4], 1, 4) == [2,3,4,1]
    @test AutomotiveSimulator.cyclic_shift_left!([1,2,3,4], 2, 4) == [3,4,1,2]
    @test AutomotiveSimulator.cyclic_shift_left!([1,2,3,4,5,6,7], 1, 4) == [2,3,4,1,5,6,7]
    @test AutomotiveSimulator.cyclic_shift_left!([1,2,3,4,5,6,7], 2, 4) == [3,4,1,2,5,6,7]

    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)
    scene = Scene(Entity{VehicleState, VehicleDef, Int64})

    col = get_first_collision(trajdata[1])
    @test col.is_colliding
    @test col.A == 1
    @test col.B == 2

    @test get_first_collision(trajdata[1], CPAMemory()).is_colliding == true
    scene = trajdata[1]
    @test is_collision_free(scene) == false
    @test is_collision_free(scene, [1]) == false
    @test is_colliding(scene[1], scene[2])
    @test get_distance(scene[1], scene[2]) == 0
    @test is_collision_free(trajdata[2])
    get_distance(scene[1], scene[2])

    roadway = gen_straight_roadway(2, 100.0)
    veh1 = Entity(VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 10.0), VehicleDef(), 1)
    veh2 = Entity(VehicleState(VecSE2(10.0, 0.0, 0.0), roadway, 5.0), VehicleDef(), 2)
    scene = Scene([veh1, veh2])
    @test is_collision_free(scene)
    @test get_distance(veh1, veh2) ≈ 6.0
end

@testset "parallel axis" begin
    ## Series of test 1: Far field
    thetaspace = LinRange(0.0, 2*float(pi), 10)
    xspace = LinRange(15.0, 300.0, 10)
    yspace = LinRange(10.0, 300.0, 10)
    @test no_collisions(xspace, yspace, thetaspace)

    ## Series of test 2: Superposition
    thetaspace = LinRange(0.0, 2*float(pi), 10)
    xspace = LinRange(-2.0, 2.0, 10)
    yspace = LinRange(1.5, 1.5, 10)

    @test collisions(xspace,yspace,thetaspace)

    ## Close but no collisions
    veh2 = create_vehicle(0.0, 2.2, 0.0, id=2)
    @test !collision_checker(VEH_REF, veh2)
    veh2 = create_vehicle(3.1, 2.2, 1.0pi, id=2)
    @test !collision_checker(VEH_REF, veh2)

    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)
    scene = trajdata[1]
    @test collision_checker(scene, 1)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: acm

using AWS.Compat
using AWS.UUIDs
"""
    AddTagsToCertificate()

Adds one or more tags to an ACM certificate. Tags are labels that you can use to identify and organize your AWS resources. Each tag consists of a key and an optional value. You specify the certificate on input by its Amazon Resource Name (ARN). You specify the tag by using a key-value pair.  You can apply a tag to just one certificate if you want to identify a specific characteristic of that certificate, or you can apply the same tag to multiple certificates if you want to filter for a common relationship among those certificates. Similarly, you can apply the same tag to multiple resources if you want to specify a relationship among those resources. For example, you can add the same tag to an ACM certificate and an Elastic Load Balancing load balancer to indicate that they are both used by the same website. For more information, see Tagging ACM certificates.  To remove one or more tags, use the RemoveTagsFromCertificate action. To view all of the tags that have been applied to the certificate, use the ListTagsForCertificate action.

# Required Parameters
- `CertificateArn`: String that contains the ARN of the ACM certificate to which the tag is to be applied. This must be of the form:  arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012  For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.
- `Tags`: The key-value pair that defines the tag. The tag value is optional.

"""

add_tags_to_certificate(CertificateArn, Tags; aws_config::AWSConfig=global_aws_config()) = acm("AddTagsToCertificate", Dict{String, Any}("CertificateArn"=>CertificateArn, "Tags"=>Tags); aws_config=aws_config)
add_tags_to_certificate(CertificateArn, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("AddTagsToCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn, "Tags"=>Tags), args)); aws_config=aws_config)

"""
    DeleteCertificate()

Deletes a certificate and its associated private key. If this action succeeds, the certificate no longer appears in the list that can be displayed by calling the ListCertificates action or be retrieved by calling the GetCertificate action. The certificate will not be available for use by AWS services integrated with ACM.   You cannot delete an ACM certificate that is being used by another AWS service. To delete a certificate that is in use, the certificate association must first be removed.

# Required Parameters
- `CertificateArn`: String that contains the ARN of the ACM certificate to be deleted. This must be of the form:  arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012  For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.

"""

delete_certificate(CertificateArn; aws_config::AWSConfig=global_aws_config()) = acm("DeleteCertificate", Dict{String, Any}("CertificateArn"=>CertificateArn); aws_config=aws_config)
delete_certificate(CertificateArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("DeleteCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn), args)); aws_config=aws_config)

"""
    DescribeCertificate()

Returns detailed metadata about the specified ACM certificate.

# Required Parameters
- `CertificateArn`: The Amazon Resource Name (ARN) of the ACM certificate. The ARN must have the following form:  arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012  For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.

"""

describe_certificate(CertificateArn; aws_config::AWSConfig=global_aws_config()) = acm("DescribeCertificate", Dict{String, Any}("CertificateArn"=>CertificateArn); aws_config=aws_config)
describe_certificate(CertificateArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("DescribeCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn), args)); aws_config=aws_config)

"""
    ExportCertificate()

Exports a private certificate issued by a private certificate authority (CA) for use anywhere. The exported file contains the certificate, the certificate chain, and the encrypted private 2048-bit RSA key associated with the public key that is embedded in the certificate. For security, you must assign a passphrase for the private key when exporting it.  For information about exporting and formatting a certificate using the ACM console or CLI, see Export a Private Certificate.

# Required Parameters
- `CertificateArn`: An Amazon Resource Name (ARN) of the issued certificate. This must be of the form:  arn:aws:acm:region:account:certificate/12345678-1234-1234-1234-123456789012
- `Passphrase`: Passphrase to associate with the encrypted exported private key. If you want to later decrypt the private key, you must have the passphrase. You can use the following OpenSSL command to decrypt a private key:   openssl rsa -in encrypted_key.pem -out decrypted_key.pem

"""

export_certificate(CertificateArn, Passphrase; aws_config::AWSConfig=global_aws_config()) = acm("ExportCertificate", Dict{String, Any}("CertificateArn"=>CertificateArn, "Passphrase"=>Passphrase); aws_config=aws_config)
export_certificate(CertificateArn, Passphrase, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("ExportCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn, "Passphrase"=>Passphrase), args)); aws_config=aws_config)

"""
    GetCertificate()

Retrieves an Amazon-issued certificate and its certificate chain. The chain consists of the certificate of the issuing CA and the intermediate certificates of any other subordinate CAs. All of the certificates are base64 encoded. You can use OpenSSL to decode the certificates and inspect individual fields.

# Required Parameters
- `CertificateArn`: String that contains a certificate ARN in the following format:  arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012  For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.

"""

get_certificate(CertificateArn; aws_config::AWSConfig=global_aws_config()) = acm("GetCertificate", Dict{String, Any}("CertificateArn"=>CertificateArn); aws_config=aws_config)
get_certificate(CertificateArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("GetCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn), args)); aws_config=aws_config)

"""
    ImportCertificate()

Imports a certificate into AWS Certificate Manager (ACM) to use with services that are integrated with ACM. Note that integrated services allow only certificate types and keys they support to be associated with their resources. Further, their support differs depending on whether the certificate is imported into IAM or into ACM. For more information, see the documentation for each service. For more information about importing certificates into ACM, see Importing Certificates in the AWS Certificate Manager User Guide.   ACM does not provide managed renewal for certificates that you import.  Note the following guidelines when importing third party certificates:   You must enter the private key that matches the certificate you are importing.   The private key must be unencrypted. You cannot import a private key that is protected by a password or a passphrase.   If the certificate you are importing is not self-signed, you must enter its certificate chain.   If a certificate chain is included, the issuer must be the subject of one of the certificates in the chain.   The certificate, private key, and certificate chain must be PEM-encoded.   The current time must be between the Not Before and Not After certificate fields.   The Issuer field must not be empty.   The OCSP authority URL, if present, must not exceed 1000 characters.   To import a new certificate, omit the CertificateArn argument. Include this argument only when you want to replace a previously imported certifica   When you import a certificate by using the CLI, you must specify the certificate, the certificate chain, and the private key by their file names preceded by file://. For example, you can specify a certificate saved in the C:temp folder as file://C:tempcertificate_to_import.pem. If you are making an HTTP or HTTPS Query request, include these arguments as BLOBs.    When you import a certificate by using an SDK, you must specify the certificate, the certificate chain, and the private key files in the manner required by the programming language you're using.    The cryptographic algorithm of an imported certificate must match the algorithm of the signing CA. For example, if the signing CA key type is RSA, then the certificate key type must also be RSA.   This operation returns the Amazon Resource Name (ARN) of the imported certificate.

# Required Parameters
- `Certificate`: The certificate to import.
- `PrivateKey`: The private key that matches the public key in the certificate.

# Optional Parameters
- `CertificateArn`: The Amazon Resource Name (ARN) of an imported certificate to replace. To import a new certificate, omit this field.
- `CertificateChain`: The PEM encoded certificate chain.
- `Tags`: One or more resource tags to associate with the imported certificate.  Note: You cannot apply tags when reimporting a certificate.
"""

import_certificate(Certificate, PrivateKey; aws_config::AWSConfig=global_aws_config()) = acm("ImportCertificate", Dict{String, Any}("Certificate"=>Certificate, "PrivateKey"=>PrivateKey); aws_config=aws_config)
import_certificate(Certificate, PrivateKey, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("ImportCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Certificate"=>Certificate, "PrivateKey"=>PrivateKey), args)); aws_config=aws_config)

"""
    ListCertificates()

Retrieves a list of certificate ARNs and domain names. You can request that only certificates that match a specific status be listed. You can also filter by specific attributes of the certificate. Default filtering returns only RSA_2048 certificates. For more information, see Filters.

# Optional Parameters
- `CertificateStatuses`: Filter the certificate list by status value.
- `Includes`: Filter the certificate list. For more information, see the Filters structure.
- `MaxItems`: Use this parameter when paginating results to specify the maximum number of items to return in the response. If additional items exist beyond the number you specify, the NextToken element is sent in the response. Use this NextToken value in a subsequent request to retrieve additional items.
- `NextToken`: Use this parameter only when paginating results and only in a subsequent request after you receive a response with truncated results. Set it to the value of NextToken from the response you just received.
"""

list_certificates(; aws_config::AWSConfig=global_aws_config()) = acm("ListCertificates"; aws_config=aws_config)
list_certificates(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("ListCertificates", args; aws_config=aws_config)

"""
    ListTagsForCertificate()

Lists the tags that have been applied to the ACM certificate. Use the certificate's Amazon Resource Name (ARN) to specify the certificate. To add a tag to an ACM certificate, use the AddTagsToCertificate action. To delete a tag, use the RemoveTagsFromCertificate action.

# Required Parameters
- `CertificateArn`: String that contains the ARN of the ACM certificate for which you want to list the tags. This must have the following form:  arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012  For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.

"""

list_tags_for_certificate(CertificateArn; aws_config::AWSConfig=global_aws_config()) = acm("ListTagsForCertificate", Dict{String, Any}("CertificateArn"=>CertificateArn); aws_config=aws_config)
list_tags_for_certificate(CertificateArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("ListTagsForCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn), args)); aws_config=aws_config)

"""
    RemoveTagsFromCertificate()

Remove one or more tags from an ACM certificate. A tag consists of a key-value pair. If you do not specify the value portion of the tag when calling this function, the tag will be removed regardless of value. If you specify a value, the tag is removed only if it is associated with the specified value.  To add tags to a certificate, use the AddTagsToCertificate action. To view all of the tags that have been applied to a specific ACM certificate, use the ListTagsForCertificate action.

# Required Parameters
- `CertificateArn`: String that contains the ARN of the ACM Certificate with one or more tags that you want to remove. This must be of the form:  arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012  For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.
- `Tags`: The key-value pair that defines the tag to remove.

"""

remove_tags_from_certificate(CertificateArn, Tags; aws_config::AWSConfig=global_aws_config()) = acm("RemoveTagsFromCertificate", Dict{String, Any}("CertificateArn"=>CertificateArn, "Tags"=>Tags); aws_config=aws_config)
remove_tags_from_certificate(CertificateArn, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("RemoveTagsFromCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn, "Tags"=>Tags), args)); aws_config=aws_config)

"""
    RenewCertificate()

Renews an eligable ACM certificate. At this time, only exported private certificates can be renewed with this operation. In order to renew your ACM PCA certificates with ACM, you must first grant the ACM service principal permission to do so. For more information, see Testing Managed Renewal in the ACM User Guide.

# Required Parameters
- `CertificateArn`: String that contains the ARN of the ACM certificate to be renewed. This must be of the form:  arn:aws:acm:region:123456789012:certificate/12345678-1234-1234-1234-123456789012  For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.

"""

renew_certificate(CertificateArn; aws_config::AWSConfig=global_aws_config()) = acm("RenewCertificate", Dict{String, Any}("CertificateArn"=>CertificateArn); aws_config=aws_config)
renew_certificate(CertificateArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("RenewCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn), args)); aws_config=aws_config)

"""
    RequestCertificate()

Requests an ACM certificate for use with other AWS services. To request an ACM certificate, you must specify a fully qualified domain name (FQDN) in the DomainName parameter. You can also specify additional FQDNs in the SubjectAlternativeNames parameter.  If you are requesting a private certificate, domain validation is not required. If you are requesting a public certificate, each domain name that you specify must be validated to verify that you own or control the domain. You can use DNS validation or email validation. We recommend that you use DNS validation. ACM issues public certificates after receiving approval from the domain owner.

# Required Parameters
- `DomainName`:  Fully qualified domain name (FQDN), such as www.example.com, that you want to secure with an ACM certificate. Use an asterisk (*) to create a wildcard certificate that protects several sites in the same domain. For example, *.example.com protects www.example.com, site.example.com, and images.example.com.   The first domain name you enter cannot exceed 64 octets, including periods. Each subsequent Subject Alternative Name (SAN), however, can be up to 253 octets in length.

# Optional Parameters
- `CertificateAuthorityArn`: The Amazon Resource Name (ARN) of the private certificate authority (CA) that will be used to issue the certificate. If you do not provide an ARN and you are trying to request a private certificate, ACM will attempt to issue a public certificate. For more information about private CAs, see the AWS Certificate Manager Private Certificate Authority (PCA) user guide. The ARN must have the following form:   arn:aws:acm-pca:region:account:certificate-authority/12345678-1234-1234-1234-123456789012
- `DomainValidationOptions`: The domain name that you want ACM to use to send you emails so that you can validate domain ownership.
- `IdempotencyToken`: Customer chosen string that can be used to distinguish between calls to RequestCertificate. Idempotency tokens time out after one hour. Therefore, if you call RequestCertificate multiple times with the same idempotency token within one hour, ACM recognizes that you are requesting only one certificate and will issue only one. If you change the idempotency token for each call, ACM recognizes that you are requesting multiple certificates.
- `Options`: Currently, you can use this parameter to specify whether to add the certificate to a certificate transparency log. Certificate transparency makes it possible to detect SSL/TLS certificates that have been mistakenly or maliciously issued. Certificates that have not been logged typically produce an error message in a browser. For more information, see Opting Out of Certificate Transparency Logging.
- `SubjectAlternativeNames`: Additional FQDNs to be included in the Subject Alternative Name extension of the ACM certificate. For example, add the name www.example.net to a certificate for which the DomainName field is www.example.com if users can reach your site by using either name. The maximum number of domain names that you can add to an ACM certificate is 100. However, the initial quota is 10 domain names. If you need more than 10 names, you must request a quota increase. For more information, see Quotas.  The maximum length of a SAN DNS name is 253 octets. The name is made up of multiple labels separated by periods. No label can be longer than 63 octets. Consider the following examples:     (63 octets).(63 octets).(63 octets).(61 octets) is legal because the total length is 253 octets (63+1+63+1+63+1+61) and no label exceeds 63 octets.    (64 octets).(63 octets).(63 octets).(61 octets) is not legal because the total length exceeds 253 octets (64+1+63+1+63+1+61) and the first label exceeds 63 octets.    (63 octets).(63 octets).(63 octets).(62 octets) is not legal because the total length of the DNS name (63+1+63+1+63+1+62) exceeds 253 octets.
- `Tags`: One or more resource tags to associate with the certificate.
- `ValidationMethod`: The method you want to use if you are requesting a public certificate to validate that you own or control domain. You can validate with DNS or validate with email. We recommend that you use DNS validation.
"""

request_certificate(DomainName; aws_config::AWSConfig=global_aws_config()) = acm("RequestCertificate", Dict{String, Any}("DomainName"=>DomainName); aws_config=aws_config)
request_certificate(DomainName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("RequestCertificate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DomainName"=>DomainName), args)); aws_config=aws_config)

"""
    ResendValidationEmail()

Resends the email that requests domain ownership validation. The domain owner or an authorized representative must approve the ACM certificate before it can be issued. The certificate can be approved by clicking a link in the mail to navigate to the Amazon certificate approval website and then clicking I Approve. However, the validation email can be blocked by spam filters. Therefore, if you do not receive the original mail, you can request that the mail be resent within 72 hours of requesting the ACM certificate. If more than 72 hours have elapsed since your original request or since your last attempt to resend validation mail, you must request a new certificate. For more information about setting up your contact email addresses, see Configure Email for your Domain.

# Required Parameters
- `CertificateArn`: String that contains the ARN of the requested certificate. The certificate ARN is generated and returned by the RequestCertificate action as soon as the request is made. By default, using this parameter causes email to be sent to all top-level domains you specified in the certificate request. The ARN must be of the form:   arn:aws:acm:us-east-1:123456789012:certificate/12345678-1234-1234-1234-123456789012
- `Domain`: The fully qualified domain name (FQDN) of the certificate that needs to be validated.
- `ValidationDomain`: The base validation domain that will act as the suffix of the email addresses that are used to send the emails. This must be the same as the Domain value or a superdomain of the Domain value. For example, if you requested a certificate for site.subdomain.example.com and specify a ValidationDomain of subdomain.example.com, ACM sends email to the domain registrant, technical contact, and administrative contact in WHOIS and the following five addresses:   admin@subdomain.example.com   administrator@subdomain.example.com   hostmaster@subdomain.example.com   postmaster@subdomain.example.com   webmaster@subdomain.example.com

"""

resend_validation_email(CertificateArn, Domain, ValidationDomain; aws_config::AWSConfig=global_aws_config()) = acm("ResendValidationEmail", Dict{String, Any}("CertificateArn"=>CertificateArn, "Domain"=>Domain, "ValidationDomain"=>ValidationDomain); aws_config=aws_config)
resend_validation_email(CertificateArn, Domain, ValidationDomain, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("ResendValidationEmail", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn, "Domain"=>Domain, "ValidationDomain"=>ValidationDomain), args)); aws_config=aws_config)

"""
    UpdateCertificateOptions()

Updates a certificate. Currently, you can use this function to specify whether to opt in to or out of recording your certificate in a certificate transparency log. For more information, see  Opting Out of Certificate Transparency Logging.

# Required Parameters
- `CertificateArn`: ARN of the requested certificate to update. This must be of the form:  arn:aws:acm:us-east-1:account:certificate/12345678-1234-1234-1234-123456789012
- `Options`: Use to update the options for your certificate. Currently, you can specify whether to add your certificate to a transparency log. Certificate transparency makes it possible to detect SSL/TLS certificates that have been mistakenly or maliciously issued. Certificates that have not been logged typically produce an error message in a browser.

"""

update_certificate_options(CertificateArn, Options; aws_config::AWSConfig=global_aws_config()) = acm("UpdateCertificateOptions", Dict{String, Any}("CertificateArn"=>CertificateArn, "Options"=>Options); aws_config=aws_config)
update_certificate_options(CertificateArn, Options, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = acm("UpdateCertificateOptions", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("CertificateArn"=>CertificateArn, "Options"=>Options), args)); aws_config=aws_config)


text:
@doc raw"
`struct RK_H0{T <: AbstractFloat} <: ReproducingKernel_0`

Defines a type of reproducing kernel of Bessel Potential space ``H^{n/2 + 1/2}_ε (R^n)`` ('Basic Matérn kernel'):
```math
V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|) \, .
```
# Fields
- `ε::T`: 'scaling parameter' from the Bessel Potential space definition,
           it may be omitted in the struct constructor otherwise it must be greater than zero
"
struct RK_H0{T <: AbstractFloat} <: ReproducingKernel_0
     ε::T
     RK_H0() = new{Float64}(0.0)
     function RK_H0(ε::T) where T <: AbstractFloat
        if ε <= 0
          throw(DomainError(ε, "Parameter ε must be a positive number."))
        end
        new{T}(ε)
     end
     function RK_H0(ε::Integer)
        if ε <= 0
           throw(DomainError(ε, "Parameter ε must be a positive number."))
        end
        new{Float64}(convert(Float64, ε))
     end
end

@doc raw"
`struct RK_H1{T <: AbstractFloat} <: ReproducingKernel_1`

Defines a type of reproducing kernel of Bessel Potential space ``H^{n/2 + 3/2}_ε (R^n)`` ('Linear Matérn kernel'):
```math
V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|)
             (1 + \varepsilon |\xi  - \eta|) \, .
```
# Fields
- `ε::T`: 'scaling parameter' from the Bessel Potential space definition,
           it may be omitted in the struct constructor otherwise it must be greater than zero
"
struct RK_H1{T <: AbstractFloat} <: ReproducingKernel_1
     ε::T
     RK_H1() = new{Float64}(0.0)
     function RK_H1(ε::T) where T <: AbstractFloat
        if ε <= T(0.0)
           throw(DomainError(ε, "Parameter ε must be a positive number."))
        end
        new{T}(ε)
     end
     function RK_H1(ε::Integer)
        if ε <= 0
           throw(DomainError(ε, "Parameter ε must be a positive number."))
        end
        new{Float64}(convert(Float64, ε))
     end
end

@doc raw"
`struct RK_H2{T <: AbstractFloat} <: ReproducingKernel_2`

Defines a type of reproducing kernel of Bessel Potential space ``H^{n/2 + 5/2}_ε (R^n)`` ('Quadratic Matérn kernel'):
```math
V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|)
             (3 + 3\varepsilon |\xi  - \eta| + \varepsilon ^2 |\xi - \eta| ^2 ) \, .
```
# Fields
- `ε::T`: 'scaling parameter' from the Bessel Potential space definition,
           it may be omitted in the struct constructor otherwise it must be greater than zero
"
struct RK_H2{T <: AbstractFloat} <: ReproducingKernel_2
     ε::T
     RK_H2() = new{Float64}(0.0)
     function RK_H2(ε::T) where T <: AbstractFloat
        if ε <= T(0.0)
           throw(DomainError(ε, "Parameter ε must be a positive number."))
        end
        new{T}(ε)
     end
     function RK_H2(ε::Integer)
        if ε <= 0
           throw(DomainError(ε, "Parameter ε must be a positive number."))
        end
        new{Float64}(convert(Float64, ε))
     end
end

@inline function _rk(kernel::RK,
                     η::Vector{T},
                     ξ::Vector{T}
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
   defined::Bool  = false
   x::T = kernel.ε * norm(ξ - η)
   if isa(kernel, RK_H2)
      defined = true
      value = (T(3.0) + x * (T(3.0) + x)) * exp(-x)
   end
   if isa(kernel, RK_H1)
      defined = true
      value = (T(1.0) + x) * exp(-x)
   end
   if isa(kernel, RK_H0)
      defined = true
      value = exp(-x)
   end
   if !defined
      Throw(ArgumentError("`kernel`: Incorrect parameter type."))
   end
   return value
end

@inline function _∂rk_∂e(kernel::RK,
                         η::Vector{T},
                         ξ::Vector{T},
                         e::Vector{T},
                        ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
   value::T = T(0.0)
   defined::Bool  = false
   t = η .- ξ
   x = kernel.ε * norm(t)
   s = sum(t .* e)
   if isa(kernel, RK_H2)
      defined = true
      value = kernel.ε^2 * exp(-x) * (T(1.0) + x) * s
   end
   if isa(kernel, RK_H1)
      defined = true
      value = kernel.ε^2 * exp(-x) * s
   end
   if !defined
      Throw(ArgumentError("`kernel`: Incorrect parameter type."))
   end
   return value
end

@inline function _∂rk_∂η_k(kernel::RK,
                           η::Vector{T},
                           ξ::Vector{T},
                           k::Int
                          ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
#  Note: Derivative of spline built with reproducing kernel RK_H0 does not exist at the spline nodes.
   value::T = T(0.0)
   defined::Bool  = false
   normt = norm(η - ξ)
   x = kernel.ε * normt
   if isa(kernel, RK_H2)
      defined = true
      value = kernel.ε^2 * exp(-x) * (T(1.0) + x) * (ξ[k] - η[k])
   end
   if isa(kernel, RK_H1)
      defined = true
      value = kernel.ε^2 * exp(-x) * (ξ[k] - η[k])
   end
   if isa(kernel, RK_H0)
      if normt < sqrt(eps(T))
         Throw(ArgumentError("Derivative does not exist at this point."))
      end
      defined = true
      value = kernel.ε * exp(-x) * (ξ[k] - η[k]) / normt
   end
   if !defined
      Throw(ArgumentError("`kernel`: Incorrect parameter type."))
   end
   return value
end

@inline function _∂²rk_∂η_r_∂ξ_k(kernel::RK,
                                η::Vector{T},
                                ξ::Vector{T},
                                r::Int,
                                k::Int
                               ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
   value::T = T(0.0)
   defined::Bool  = false
   if isa(kernel, RK_H2)
      defined = true
      x = kernel.ε * norm(η - ξ)
      if r == k
         if x > T(0.0)
            value = kernel.ε^2 * exp(-x) * (T(1.0) + x - (kernel.ε * (ξ[r] - η[r]))^2)
         else
            value = kernel.ε^2
         end
      else
         if x > T(0.0)
            value = -kernel.ε^4 * exp(-x) * (ξ[r] - η[r]) * (ξ[k] - η[k])
         else
            value = T(0.0)
         end
      end
   end
   if isa(kernel, RK_H1)
      defined = true
      t = norm(η - ξ)
      x = kernel.ε * t
      if r == k
         if t > T(0.0)
            value = kernel.ε^2 * exp(-x) * (T(1.0) - kernel.ε * (ξ[r] - η[r])^2 / t)
         else
            value = kernel.ε^2
         end
      else
         if t > T(0.0)
            value = -kernel.ε^3 * exp(-x) * (ξ[r] - η[r]) * (ξ[k] - η[k]) / t
         else
            value = T(0.0)
         end
      end
   end
   if !defined
      Throw(ArgumentError("`kernel`: Incorrect parameter type."))
   end
   return value
end


text:
module XmlToDict

using LightXML
using DataStructures

export xmlparse_string, @x2d_str, XmlDictAccessor, getindex

function handle_node(xelem:: XMLElement; strip_text=true, wrapper_text=false)

    elem = OrderedDict{String,Any}()
    for a in attributes(xelem)
        elem["@"*name(a)] = value(a)
    end

    has_element_children = false
    for item in child_elements(xelem)
        cname = "$(name(item))"
        celem = handle_node(item)

        if cname in keys(elem)
            curr = elem[cname]
            if isa(curr, Array)
                push!(elem[cname],celem)
            else
                elem[cname] = [curr, celem]
            end
        else
            elem[cname] = celem
        end
        has_element_children = true
    end

    text = ""
    for item in child_nodes(xelem) |> xs->filter(is_textnode, xs)
        ctext = content(item)
        ctext_stripped = strip(ctext)
        text *= (strip_text || (!wrapper_text&&has_element_children)) ? ctext_stripped : ctext
    end
    # println("text: $text")

    if !has_element_children || (has_element_children && !wrapper_text && !isempty(text) )
        elem["#text"] = text
    end

    return elem
end

function xmlparse_string(xml_input)

    xdoc = parse_string(xml_input)
    xroot = root(xdoc)

    elem = OrderedDict{String,Any}()
    elem[name(xroot)] = handle_node(xroot)
    return elem
end

type XmlDictAccessor
    value::String
end

macro x2d_str(str)
    return :( XmlDictAccessor($str) )
end

function getindex(od::OrderedDict, name::XmlDictAccessor)
    if endswith(name.value, "+")
		rname = rstrip(name.value, '+')
        rvalue = od[rname]
		return isa(rvalue, Array) ? rvalue : [ rvalue ]
    elseif endswith(name.value, "*")
		rname = rstrip(name.value, '*')
		if rname in keys(od)
	        rvalue = od[rname]
			return isa(rvalue, Array) ? rvalue : [ rvalue ]
		else
			return []
		end
    else
        return od[name.value]
    end
end

end # module