text:
"""
    dynamicRM(trr::tree, RM::Array{Any,1} )

Implements the nested version of RM for the tree specified by trr.
The user has to specify a vector of risk measures of lenght N, where N is the number of stages of trr.
Furthermore the risk measures RM[i] should be of the form RM[i](states, prob).

"""
function dynamicRM(trr::Tree; RM::Array{Any,1})

    if !(length(RM) == maximum(stage(trr)) )
        println("Error, number of risk measures should be equal to the number of stages")
        return(nothing)
    end
# recursively go through tree and apply RM
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]
        T0.state[i] = RM[i](states,prob)
    end

    return(T0[1],T0)
end

"""
    VaR(trr::Tree,alpha::Float32)
implements the nested version of the Value-at-Risk at level ``\\alpha`` for the Tree ``trr``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version.
"""
function VaR(trr::Tree,alpha::Float32)
    # look for smallest x such that  P(-states <= x) >alpha , i.e.
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.VaR(states,prob,alpha)
    end
    return(T0.state[1],T0)
end

"""
    CTE(trr::Tree,alpha::Float32)
implements the nested version of the Conditional Tail Expectation at level ``\\alpha`` for the Tree ``trr``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version.
"""
function CTE(trr::Tree,alpha::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]
        println(states)
        println(prob)
        T0.state[i] = RiskMeasures.CTE(states,prob,alpha)
        println("risk value:")
        println(T0.state[i])
    end
    return(T0.state[1],T0)
end

"""
    EVaR2(trr::Tree,beta::Float32)
implements the nested version of the Entropic Value-at-Risk at level ``\\beta`` for the Tree ``trr``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version.
"""
function EVaR2(trr::Tree,beta::Float32)
	T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.EVaR2(states,prob,beta)[1]
    end
    return(T0.state[1],T0)
end

"""
    EVaR(trr::Tree,beta::Float32)
implements the nested version of the Entropic Value-at-Risk at level ``\\beta`` for the Tree ``trr``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version. This version uses JuMP and Ipopt.
"""
function EVaR(trr::Tree,beta::Float32)
	T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.EVaR(states,prob,beta)[1]
    end
    return(T0.state[1],T0)
end

"""
    AVaR(trr::Tree,beta::Float32)
implements the nested version of the Average Value-at-Risk at level ``\\alpha`` for the Tree ``trr``.
See [here](https://rubsc.github.io/riskMeasures.jl/dev/) for the univariate version.
"""
function AVaR(trr::Tree, alpha::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]
        println(states)
        println(prob)

        T0.state[i] = RiskMeasures.AVaR(states,prob,alpha)[1]
    end
    return(T0.state[1],T0)
end

#########################################
# From basicRM.jl
"""
    Expectation(states,prob)

Computes the expectation for the random variable ``Y`` with values given by `states` under
the probability measure ``Q`` given by `prob`:

```math
\\mathbb{E} Y = states \\cdot prob
```
"""
function Expectation(trr::Tree, alpha::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.Expectation(states,prob)
    end
    return(T0.state[1],T0)
end

"""
    entropic(states, prob,theta::Float32)

implements the entropic risk measure defined as
```math
\\rho_\\theta(Y) = \\theta \\cdot \\log \\mathbb{E} e^{\\frac{Y}{\\theta}},
```
where ``\\theta`` is greater 0 and ``Y`` is the random variable defined by `states` and `prob`.
"""
function entropic(trr::Tree, theta::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.entropic(states,prob,theta)
    end
    return(T0.state[1],T0)
end

#########################################################
#Mean Semi-deviation
#########################################################
"""
    mSD(states,prob,beta::Float32,p::Float32)

implements the mean semi-deviation of order `p` which is a coherent risk measure defined by
```math
mSD_\\beta^p (Y) = \\mathbb{E} Y  + \\beta \\lvert \\left( Y - \\mathbb{E}Y \\right)_+ \\rvert_p,
```
for the random variable ``Y`` defined by `states` and `prob`.
"""
function mSD(trr::Tree, beta::Float32,p::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures. mSD(states,prob,beta::Float32,p::Float32)
    end
    return(T0.state[1],T0)
end

"""
    meanVariance(states, prob,c::Float32)

implements the mean Variance risk measure defined by
```math
\\rho_c(Y) := \\mathbb{E}Y + c \\cdot \\mathbb{E} \\left( Y- \\mathbb{E}Y)^2 \\right),
```
where ``c >0`` and ``Y`` is the random variable defined by `states` and `prob`.
"""
function meanVariance(trr::Tree, c::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.meanVariance(states, prob,c::Float32)
    end
    return(T0.state[1],T0)
end

"""
    meanDeviation(states, prob,c::Float32,p::Float32)

implements the mean deviation risk measure of order ``p\\geq 1`` defined by
```math
\\rho_c^p(Y) := \\mathbb{E}Y + c \\cdot \\lVert \\left( Y- \\mathbb{E}Y \\right)^2 \\rVert_p,
```
where ``c >0`` and ``Y`` is the random variable defined by `states` and `prob`.
"""
function meanDeviation(trr::Tree, c::Float32,p::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.meanDeviation(states, prob,c::Float32,p::Float32)
    end
    return(T0.state[1],T0)
end

"""
    meanSemiVariance(states, prob,c::Float32,t::Float32)

implements the mean upper-semi variance risk measure from a target ``t`` defined by
```math
\\rho_{c,t}(Y) = \\mathbb{E}Y + c \\cdot \\mathbb{E} \\left( Y - t \\right)^2_+ ,
```
where ``c >0`` and ``Y`` is the random variable defined by `states` and `prob`.
"""
function meanSemiVariance(trr::Tree, c::Float32,t::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.meanSemiVariance(states, prob,c::Float32,t::Float32)
    end
    return(T0.state[1],T0)
end

"""
    meanSemiDevi(states, prob,c::Float32,target::Float32,p::Float32)

implements the mean upper-semi variance risk measure of order ``p \\geq 1`` from a target ``t`` defined by
```math
\\rho_{c,t}^p(Y) = \\mathbb{E}Y + c \\cdot \\lVert \\left( Y - t \\right)_+ \\rVert_p ,
```
where ``c >0`` and ``Y`` is the random variable defined by `states` and `prob`.
"""
function meanSemiDevi(trr::Tree, c::Float32,target::Float32,p::Float32)
    T0 = deepcopy(trr);

    # For every parent node going backwards
    tmp = reverse(unique(T0.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        #Risk measure w/ adjustment of values
        states = T0.state[T0.children[i+1]]
        prob = T0.probability[T0.children[i+1]]

        T0.state[i] = RiskMeasures.meanSemiDevi(states, prob,c::Float32,target::Float32,p::Float32)
    end
    return(T0.state[1],T0)
end


text:
export timebasisc0d1
export timebasiscxd0
export timebasisdelta
export timebasisshiftedlagrange

using Compat

type MonomialBasis{T,Degree,NF} <: RefSpace{T,NF} end

valuetype{T}(::MonomialBasis{T}) = T
degree{T,D}(::MonomialBasis{T,D}) = D
numfunctions(x::MonomialBasis) = degree(x)+1
@compat function (x::MonomialBasis)(p)
    t = cartesian(p)[1]
    valuetype(x)[t^d for d in 0 : degree(x)]
end

abstract type AbstractTimeBasisFunction end

"""
    TimeBasisFunction{N,D}

T: the value type of the time basis function
N: the number of intervals in the support (this included the semi infinite interval
    stretching to +∞)
D1: the degree of the TBF restricted to each of the intervals **plus one**
"""
type TimeBasisFunction{T,N,D1,D} <: AbstractTimeBasisFunction
    timestep::T
    numfunctions::Int
    polys::SVector{N,Polynomial{D1,T}}
end

function derive{T,N,D1,D}(tbf::TimeBasisFunction{T,N,D1,D})
    dpolys = Polynomial{D1-1,T}[ derive(p) for p in tbf.polys ]
    TimeBasisFunction{T,N,D1-1,D-1}(
        tbf.timestep,
        tbf.numfunctions,
        SVector{N,Polynomial{D1-1,T}}(dpolys...)
    )
end

type TimeBasisDelta{T} <: AbstractTimeBasisFunction
    timestep::T
    numfunctions::Int
end

type DiracBoundary{T} <: RefSpace{T,1} end
numfunctions(x::DiracBoundary) = 1

scalartype{T}(x::TimeBasisDelta{T}) = T
numfunctions(x::TimeBasisDelta) = x.numfunctions
refspace{T}(x::TimeBasisDelta{T}) = DiracBoundary{T}()
timestep(x::TimeBasisDelta) = x.timestep
numintervals(x::TimeBasisDelta) = 1

timebasisdelta(dt, Nt)= TimeBasisDelta(dt, Nt)

numintervals{T,N,D1,D}(tbf::TimeBasisFunction{T,N,D1,D}) = N
degree{T,N,D1,D}(tbf::TimeBasisFunction{T,N,D1,D}) = D
timestep(tbf::TimeBasisFunction) = tbf.timestep
scalartype{T,N,D1,D}(tbf::TimeBasisFunction{T,N,D1,D}) = T

numfunctions(t::TimeBasisFunction) = t.numfunctions
refspace{T,N,D1,D}(t::TimeBasisFunction{T,N,D1,D}) = MonomialBasis{T,D,D1}()

geometry(t::AbstractTimeBasisFunction) = [SegmentedAxis(timestep(t), numfunctions(t)+numintervals(t)-3)]

"""
    timebasisc0d1(type, timestep, numfunctions)

Build the space of continuous, piecewise linear time basis functions. The DoFs
are the time steps. `numfunctions` basis functions will be built in total.
"""
function timebasisc0d1(timestep, numfunctions, T::Type=Float64)
    i, z = one(T), zero(T)
    polys = SVector(
        Polynomial(SVector(i,  i/timestep)), # 1 + t
        Polynomial(SVector(i, -i/timestep)), # 1 - t
        Polynomial(SVector(z,  z))
    )
    TimeBasisFunction{T,3,2,1}(timestep, numfunctions, polys)
end

"""
    timebasiscxd0(timestep, numfunctions, T::Type=Float64)

Create a temporal basis based on shifted copies of the nodal continuous, piecewise
linear interpolant.
"""
function timebasiscxd0(timestep, numfunctions, T::Type=Float64)
    i, z = one(T), zero(T)
    polys = SVector(
        Polynomial(i), # 1
        Polynomial(z), # 0
    )
    TimeBasisFunction{T,2,1,0}(timestep, numfunctions, polys)
end

"""
    timebasisspline2(timestep, numfunctions, T::Type=Float64)

Create a temporal basis based on shifted copies of the quadratic spline. The
spline is the convolution of a cxd0 and a c0d1 basis function.
"""
function timebasisspline2(dt, numfunctions, T::Type=Float64)
    i, z = one(T), zero(T)
    polys = SVector{4,Polynomial{3,T}}(
        Polynomial(SVector(i/2, i/dt, i/2/dt/dt)),
        Polynomial(SVector(i/2, i/dt, -i/dt/dt)),
        Polynomial(SVector(2*i, -2*i/dt, i/2/dt/dt)),
        Polynomial(SVector(z, z, z))
    )
    TimeBasisFunction{T,4,3,2}(dt, numfunctions, polys)
end

function timebasisshiftedlagrange(dt, numfunctions, degree, T::Type=Float64)
    z, i = zero(T), one(T)
    c = Polynomial(SVector(i))
    t = Polynomial(SVector(z,i/dt))
    polys = Polynomial{degree+1,T}[]
    for k = 0:degree
        f = c
        for i in 1:k
            f = (1/i) * f * (i - t)
        end
        g = c
        for i in 1:(degree-k)
            g = (1/i) * g * (i + t)
        end
        push!(polys, f*g)
    end
    push!(polys, 0*(t^degree))
    @assert length(polys) == degree+2
    polys = SVector{degree+2,Polynomial{degree+1,T}}(polys...)
    TimeBasisFunction{T,degree+2,degree+1,degree}(dt, numfunctions, polys)
end

function (f::TimeBasisFunction)(t::Real)
    dt = timestep(f)
    ni = numintervals(f)
    t < -dt && return zero(t)
    i = floor(Int, (t/dt)) + 2
    p = i <= ni ? f.polys[i] : f.polys[end]
    return p(t)
end

function assemblydata(tbf::TimeBasisFunction)

    T = scalartype(tbf)
    Δt = timestep(tbf)
    z = (0, zero(T))

    t = Polynomial(zero(T), one(T))

    num_cells = numfunctions(tbf)
    num_refs  = degree(tbf)+1

    max_num_funcs = numintervals(tbf)
    numfuncs = zeros(Int, num_cells, num_refs)
    data = fill(z, max_num_funcs, num_refs, num_cells)

    els = [ simplex(point((i-1)*Δt),point(i*Δt)) for i in 1:num_cells ]

    for k in 1 : numfunctions(tbf)
        tk = (k-1) * Δt
        for i in 1 : numintervals(tbf)
            # Focus on interval [(i-2)Δt,(i-1)Δt]
            p = tbf.polys[i]
            q = substitute(p,t-tk)

            c = k + i - 2
            1 <= c <= num_cells || continue
            for d = 0 : degree(q)
                r = d + 1
                w = q[d]

                j = (numfuncs[c,r] += 1)
                data[j,r,c] = (k,w)
            end
        end
    end

    return els, AssemblyData(data)
end

function assemblydata(tbf::TimeBasisDelta)

    T = scalartype(tbf)
    Δt = timestep(tbf)

    z = zero(scalartype(tbf))
    w = one(scalartype(tbf))

    num_cells = numfunctions(tbf)
    num_refs  = 1

    max_num_funcs = 1
    num_funcs = zeros(Int, num_cells, num_refs)
    data = fill((0,z), max_num_funcs, num_refs, num_cells)

    els = [ simplex(point((i-0)*Δt),point((i+1)*Δt)) for i in 1:num_cells ]

    for k in 1 : numfunctions(tbf)
        data[1,1,k] = (k,w)
    end

    return els, AssemblyData(data)
end

function convolve(f::TimeBasisFunction, g::TimeBasisFunction)
    warn("BEAST.convolve only returns correct result for constant * continuous,linear")
    dt = timestep(f)
    fg = timebasisspline2(dt, numfunctions(f), scalartype(f))
    fg.polys = dt * fg.polys
    return fg
end

function convolve(δ::TimeBasisDelta, g::TimeBasisFunction)
    return g
end

ext: jl


text:
ts = time()

println("The system time is (in ISO 8601 format):")
println(strftime("    %F %T %Z", ts))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type AbstractLearner end

mutable struct Learner <: AbstractLearner
    model
    optimizer
    loss_func
    data::AbstractDatabunch
    callbacks
    function Learner(model,optimizer,loss_func,data;cbs=())
        cb_te = TrainEvalCallback()
        all_cbs = (cb_te,cbs...)
        ordered_cbs = Tuple(z[2] for z in sort([(cb.order, cb) for cb in all_cbs], by = x -> x[1]))
        new(model, optimizer, loss_func, data, ordered_cbs)
    end
end

model(l::Learner) = l.model
#model!(l::Learner, model) = l.model = model

optimizer(l::Learner) = l.optimizer
#optimizer!(l::Learner, optimizer) = l.optimizer = optimizer

loss_func(l::Learner) = l.loss_func
#loss_func!(l::Learner, loss_func) = l.loss_func = loss_func

data(l::Learner) = l.data
#data!(l::Learner, data) = l.data = data

callbacks(l::Learner) = l.callbacks
#callbacks!(l::Learner, callbacks) = l.callbacks = callbacks

trainevalcb(l::Learner) = l.callbacks[1]

# default function cb_blabla!(l::AbstractLearner;kwargs...)
for sym in [:before_fit, :before_epoch, :before_all_batches, :before_batch, :after_batch, :before_validate, :after_validate, :after_all_batches, :after_epoch, :after_fit]
    @eval function $(Symbol(string("cb_",sym,"!")))(l::AbstractLearner;kwargs...)
        te_cb = trainevalcb(l)
        for cb in l.callbacks
           $(Symbol(string(sym,"!")))(cb;te_cb=te_cb, kwargs...)
        end
    end
end


text:
using Gridap
using Gridap.Io
using GridapGmsh

model = GmshDiscreteModel("elasticFlag.msh")

writevtk(model,"elasticFlag")

fn = "elasticFlag.json"
to_json_file(model,fn)


text:
using SpecialFunctions: erf, erfc

"""
    ρzx(E0::Float64, Ec::Float64, elm::Element)

X-ray range in g/cm² from Eqn 6 Merlet 1994 or Eqn 7 in Merlet 1995
E0, Ec in keV, Z is the atomic number
"""
function ρzx(E0::Float64, Ec::Float64, elm::Element)::Float64
    @assert E0 > Ec "The beam energy must be larger than the edge energy!"
    @assert E0 > Ec && E0 <= 50.0
    a = 1.845e-6 * evalpoly(E0, (2.6, -0.216, 0.015, -0.000387, 0.00000501)) # Ok
    x =  2.2 - 0.0166 * E0 # Ok
    return (a * (E0^x - Ec^x) * (1.0 + 2.0 / (E0^2))) /
           ((1.078 - 0.015 * z(elm)^0.7)^(1.2 - 0.04 * E0)) # Ok
end

"""
    ρzm(ρzx::Float64, E0::Float64, Ec::Float64, em::Union{Element,Material})

Depth of the peak of the ϕ(ρz) curve in g/cm² (Merlet 1994 eqn 5 or Merlet 1995 eqn 7).
E0, Ec in keV, Z is the atomic number
"""
function ρzm(ρzx::Float64, E0::Float64, Ec::Float64, em::Union{Element,Material})::Float64
    U0, Z = E0 / Ec, z(em)
    # return ρzx * (0.1 + 0.35 * exp(-0.07 * Z)) / (1.0 + 10.0 / (U0^10.0)) # 1994
    return ρzx * (0.1 + 0.35 * exp(-0.07 * Z)) / (1.0 + 1000.0 / (Z * U0^10.0)) # 1995
end

"""
    τ(elm::Element, t::Float64)

Transmission coefficient of Zeller and Ruste 1976 (from Merlet 1994 & 1995)
"""
function τ(elm::Element, ρzm::Float64, ρzx::Float64)::Float64
    Z, t = z(elm), ρzm / ρzx
    τ1, τ2 = (1.0 - t)^(4.65 + 0.0356 * Z), (1.0 - t)^(1.112 + 0.00414 * Z^2)
    return τ1 + 4.65 * (τ2 - τ1) / evalpoly(Z, (3.54, 0.0356, -0.00414)) # Ok
end

mexp(sh::AtomicSubShell) = n(sh) == 1 ? 0.95 : 0.8

"""
    Φm(elm::Element, J::Float64, A::Float64, sh::AtomicSubShell, E0::Float64, ρzm::Float64, ρzx::Float64)

The maximum amplitude of the ϕ(ρz) curve (Located at ρzm).
Z, J, A are the atomic number, ionization potential and atomic weight (mass-fraction averaged except J which is log-averaged)
"""
function Φm(
    elm::Element,
    J::Float64, # Mean ionization potential (keV)
    A::Float64, # Atomic weight
    sh::AtomicSubShell,
    E0::Float64, # in keV
    ρzm::Float64,
    ρzx::Float64,
)::Float64
    Z, m, U0 = z(elm), mexp(sh), E0 / (0.001 * energy(sh)) # Ok
    @assert U0 > 1.0  && U0 < 1000.0
    @assert J > 0.021 && J < 1.0 # keV
    Ud = U0 * (1.0 - 1.03e5 * ρzm * Z / (E0^1.61 * J^0.3 * A)) # Ok
    pZ = -0.25 + 0.0787 * Z^0.3 # Ok
    a = 5.0 * (1.0 - 1.0 / (1.0 + Z)^0.8) * ((E0 / 30.0)^pZ) # Ok
    d1, d2 = 1.0 - m, 7.0 - 4.0 * exp(-0.1 * Z) - m # Ok
    b1 =
        (Ud / U0)^d1 *
        (log(Ud) / log(U0)) *
        (1.0 / d1) *
        (1.0 - (1.0 - 1.0 / (Ud^d1)) / (d1 * log(Ud))) # Ok
    b2 =
        (Ud / U0)^d2 *
        (log(Ud) / log(U0)) *
        (1.0 / d2) *
        (1.0 - (1.0 - 1.0 / (Ud^d2)) / (d2 * log(Ud))) # Ok
    QUd, QU0 = log(Ud) / (Ud^m), log(U0) / (U0^m) # Ok
    return τ(elm, ρzm, ρzx) * (QUd / QU0) +
           a * (0.28 * (1.0 - 0.5 * exp(-0.1 * Z)) * b1 + 0.165 * Z^0.6 * b2) # Ok
end

"""
    Φ0(mat::Material, E0::Float64, Ec::Float64)

The value of the ϕ(ρz) curve at ρz=0.
"""
function Φ0(mat::Material, E0::Float64, Ec::Float64)::Float64
    Z, U0 = z(mat), E0 / Ec
    @assert U0 > 1.0
    d1, d2, d3 = 0.02 * Z, 0.1 * Z, 0.4 * Z # Ok
    b1 = (1.0 / d1) * (1.0 - (1.0 - 1.0 / U0^d1) / (d1 * log(U0))) # Ok
    b2 = (1.0 / d2) * (1.0 - (1.0 - 1.0 / U0^d2) / (d2 * log(U0))) # Ok
    b3 = (1.0 / d3) * (1.0 - (1.0 - 1.0 / U0^d3) / (d3 * log(U0))) # Ok
    pZ = -0.25 + 0.0787 * Z^0.3 # Ok
    a =
        (1.87 * Z) *
        evalpoly(log(Z + 1.0), (0.0, -0.00391, 0.00721, -0.001067)) *
        (E0 / 30.0)^pZ # Ok
    return 1.0 + a * (0.27 * b1 + (1.1 + 5.0 / Z) * (b2 - 1.1 * b3)) # Ok
end

"""
An implementation of Merlet's XPhi matrix correction algorithm as described in Llovett 2010 and Merlet 1994/1995.
"""
struct XPhi <: MatrixCorrection
    subshell::AtomicSubShell
    material::Material
    E0::Float64 # Beam energy (eV)
    Φm::Float64  # Max amplitude
    ρzm::Float64 # ρz at max amplitude
    α::Float64 # Interior side shape parameter
    β::Float64 # Surface side shape parameter

    function XPhi(
        subshell::AtomicSubShell,
        material::Material,
        E0::Float64, # Beam energy (eV)
        Φm::Float64,  # Max amplitude
        ρzm::Float64, # ρz at max amplitude
        α::Float64, # Interior side shape parameter
        β::Float64 # Surface side shape parameter
    )
        new(subshell, material, E0, Φm, ρzm, α, β)
    end

    """
        XPhi(mat::Material, sh::AtomicSubShell, e0::Float64)

    Construct an object representing the XPhi molde for the specified material and atomic shell at the
    beam energy `e0` in eV.
    """
    function XPhi(mat::Material, sh::AtomicSubShell, e0::Float64)
        E0, Ec = 0.001 * e0, 0.001 * energy(sh)
        mk = keys(mat)
        # This implements the multi-element weighting on page 367 of 1994
        Φ0_v = Φ0(mat, E0, Ec)
        M = sum(nonneg(mat, elm) * z(elm) / a(elm, mat) for elm in mk)
        Φm_v =
            sum(
                nonneg(mat, elm) * z(elm) / a(elm, mat) *
                Φm(elm, 0.001 * J(Berger1982, elm), a(elm, mat), sh, E0,  ρzm(ρzx(E0, Ec, elm), E0, Ec, elm), ρzx(E0, Ec, elm))
                for elm in mk
            ) / M
        ρzx_v = sum(nonneg(mat, elm) * z(elm) / a(elm, mat) * ρzx(E0, Ec, elm) for elm in mk) / M
        ρzm_v = ρzm(ρzx_v, E0, Ec, mat)
        α_v, β_v = 0.46598 * (ρzx_v - ρzm_v), ρzm_v / sqrt(log(Φm_v / Φ0_v))
        @assert isapprox(Φm_v * exp(-(ρzm_v / β_v)^2), Φ0_v, rtol = 1.0e-6)
        @assert isapprox(exp(-((ρzx_v - ρzm_v) / α_v)^2), 0.01, atol = 1.0e-5)
        return new(sh, mat, e0, Φm_v, ρzm_v, α_v, β_v)
    end
end

NeXLCore.range(xp::XPhi) = xp.α * sqrt(-log(0.01)) + xp.ρzm
Base.max(xp::XPhi) = xp.Φm
ϕ0(xp::XPhi) = xp.Φm * exp(-(xp.ρzm / xp.β)^2)

"""
    Φ(xphi::XPhi, ρz::Float64) # W/o absorption
    Φ(xphi::XPhi, ρz::Float64, χ::Float64) # With absorption

The antsatz for ϕ(ρz) curve in the XPhi matrix correction algorithm.
"""
function Φ(xphi::XPhi, ρz::Float64)::Float64
    ρz >= 0.0 ? #
    (
        ρz < xphi.ρzm ? #
        xphi.Φm * exp(-((ρz - xphi.ρzm) / xphi.β)^2) : #
        xphi.Φm * exp(-((ρz - xphi.ρzm) / xphi.α)^2)
    ) : #
    0.0
end
Φ(xphi::XPhi, χ::Float64, ρz::Float64)::Float64 = Φ(xphi, ρz) * exp(-ρz * χ)

"""
    F(xphi::XPhi, 𝒜::Float64, ci::Float64, χ::Float64)

𝒜 - Instrumental, physical and other poorly known parameters
ci - Mass fraction of the i-th element

The generated intensity in the XPhi model.
"""
F(xphi::XPhi, 𝒜::Float64, ci::Float64)::Float64 = #
    0.5 * sqrt(π) * 𝒜 * ci * xphi.Φm * (xphi.α + xphi.β * erf(xphi.ρzm / xphi.β))

"""
    Fχ(xphi::XPhi, 𝒜::Float64, ci::Float64, χ::Float64)

𝒜 - Instrumental, physical and other poorly known parameters
ci - Mass fraction of the i-th element
χ - reduced mass absorption coefficient

The emitted intensity in the XPhi model.
"""
function Fχ(xphi::XPhi, 𝒜::Float64, ci::Float64, χ::Float64)::Float64
    return (0.5 * sqrt(π) * xphi.Φm) * (
        exp(χ * (0.25 * xphi.β^2 * χ - xphi.ρzm)) *
        xphi.β *
        (erf(0.5 * xphi.β * χ) + erf(xphi.ρzm / xphi.β - 0.5 * xphi.β * χ)) +
        exp(χ * (0.25 * xphi.α^2 * χ - xphi.ρzm)) * xphi.α * erfc(0.5 * xphi.α * χ)
    )
end

"""
    Fχ(xphi::XPhi, 𝒜::Float64, ci::Float64, χ::Float64, ρz0::Float64, ρz1::Float64)

𝒜 - Instrumental, physical and other poorly known parameters
ci - Mass fraction of the i-th element
χ - reduced mass absorption coefficient
ρz0, ρz1 - The mass depth range to integrate over

The emitted intensity in the XPhi model for the depth between ρz0 and ρz1
"""
function Fχp(
    xphi::XPhi,
    𝒜::Float64,
    ci::Float64,
    χ::Float64,
    ρz0::Float64,
    ρz1::Float64,
)::Float64
    @assert ρz0 >= 0.0
    @assert ρz0 < ρz1
    if ρz0 < xphi.ρzm && ρz1 > xphi.ρzm
        return Fχp(xphi, 𝒜, ci, χ, ρz0, xphi.ρzm) + Fχp(xphi, 𝒜, ci, χ, xphi.ρzm, ρz1)
    else
        ξ = (ρz0 < xphi.ρzm ? xphi.β : xphi.α)
        return (0.5 * sqrt(π) * xphi.Φm * ξ) *
               exp(χ * (0.25 * ξ^2 * χ - xphi.ρzm)) *
               (
                   erf((ρz1 - xphi.ρzm) / ξ + 0.5 * ξ * χ) -
                   erf((ρz0 - xphi.ρzm) / ξ + 0.5 * ξ * χ)
               )
    end
end

# The remainder of the code is to allow XPhi to fit within the framework with other
# matrix correction algorithms.

matrixcorrection(::Type{XPhi}, mat::Material, ashell::AtomicSubShell, e0::Float64) =
    XPhi(mat, ashell, e0)

function Fχ(xphi::XPhi, xray::CharXRay, θtoa::Float64)::Float64
    @assert inner(xray) == xphi.subshell
    return Fχ(xphi, 1.0, 1.0, χ(material(xphi), xray, θtoa))
end
F(xphi::XPhi) = F(xphi, 1.0, 1.0)
Fχp(xphi::XPhi, xray::CharXRay, θtoa::Real, t::Real) =
    Fχp(xphi, 1.0, 1.0, χ(material(xphi), xray, θtoa), 0.0, t)
ϕ(xphi::XPhi, ρz) = Φ(xphi, ρz)

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
module WeightedLS

using LinearAlgebra

greet() = print("Hello World!")

function normalLS(A::Array{Float64,2}, b::Array{Float64,1})
    return A\b
end

# On the modified Gram-Schmidt algorithm for weighted and constrained linear least squares problems
# https://github.com/borglab/gtsam/blob/develop/gtsam/linear/NoiseModel.cpp#L445
# solve min (b-Ax)^TW(b-Ax)
function weightedLS(A::Array{Float64,2}, b::Array{Float64,1}, W::Array{Float64,1})
    # TODO: input size check

    # get size, A is mxn, b is nx1, W is nx1
    m = size(A,1)
    n = size(A,2)
    max_rank = min(m,n)
    Ab = [A b]
    weights = W;
    outweights = zeros(n, 1)
    Rd = zeros(m, n+1)
    rd_row = 1
    rd = zeros(1, n+1)
    for j=1:n  # loop over columns
        a = Ab[1:m,j:j]  # get jth column of new system
        constraint_row = checkConstraint(dropdims(a,dims=2),W,m)
        if constraint_row > 0
            rd = Ab[constraint_row:constraint_row,:]      # get constraint row
            Rd[rd_row,:] = rd                             # row of R matrix
            outweights[rd_row] = weights[constraint_row]
            rd_row += 1
            if rd_row > max_rank
                break
            end

            m -= 1  # decrease number of rows?

            if constraint_row != m+1
                Ab[constraint_row,:] = Ab[m+1,:]
                weights[constraint_row,:] = weights[m+1,:]
            end
            a_reduced = Ab[1:m,j:j]

            a_reduced = a_reduced .* 1.0/rd[1,j]  # divide by the jth column
            Ab[1:m,j+1:n+1] -= a_reduced * rd[1:1,j+1:n+1]  # subtract off the rest of the row?
        else
            precision = 0.0
            pseudo = zeros(m, 1)
            rd = zeros(1, n+1)
            for i=1:m
                ai = a[i]
                if (abs(ai) > 1e-9)
                    pseudo[i,1] = weights[i]*ai
                    precision += pseudo[i]*ai
                else
                    pseudo[i,1] = 0.0
                end
            end
            if precision > 1e-8
                pseudo = pseudo ./ precision
                rd[1,j] = 1.0
                rd[1:1,j+1:n+1] = pseudo' * Ab[1:m,j+1:n+1]
                Rd[rd_row,:] = rd
                outweights[rd_row] = precision
                rd_row += 1
            else
                continue
            end
            if rd_row > max_rank
                break
            end
            Ab[1:m,j+1:n+1] -= a * rd[1:1,j+1:n+1]

        end
    end
    return Rd, outweights
end

function checkConstraint(a::Array{Float64,1}, W::Array{Float64,1}, m::Integer)
    max_element = 1e-9
    constraint_row = -1
    for i = 1:m
        if W[i] != Inf
            continue
        end
        abs_ai = abs(a[i])
        if abs_ai > max_element
            max_element = abs_ai
            constraint_row = i
        end
    end
    return constraint_row
end

end # module


text:
"""
`KernelFactors` is a module implementing separable filtering kernels,
each stored in terms of their factors. The following kernels are
supported:

  - `sobel`
  - `prewitt`
  - `ando3`, `ando4`, and `ando5` (the latter in 2d only)
  - `scharr`
  - `bickley`
  - `gaussian`
  - `IIRGaussian` (approximate gaussian filtering, fast even for large σ)

See also: [`Kernel`](@ref).
"""
module KernelFactors

using StaticArrays, OffsetArrays
using ..ImageFiltering: centered, dummyind
import ..ImageFiltering: _reshape, _vec, nextendeddims
using Base: tail, Indices, @pure, checkbounds_indices, throw_boundserror, @propagate_inbounds

# We would like to do `using ..ImageFiltering.Kernel` but we cannot because
# `kernelfactors.jl` is included before  `kernel.jl` in the ImageFiltering.jl
# file. Hence, ImageFiltering.Kernel does not yet exist when this module is
# loaded. The reason we want to have the Kernel module defined is so that that
# Documenter.jl (the documentation system) can parse a reference such as `See
# also: [`Kernel.ando3`](@ref)`. With the more general `using ImageFiltering`,
# we seem to sidestep this scope problem, although I don't actually understand
# the mechanism form why this works. - ZS
using ImageFiltering

abstract type IIRFilter{T} end

Base.eltype(kernel::IIRFilter{T}) where {T} = T

"""
    ReshapedOneD{N,Npre}(data)

Return an object of dimensionality `N`, where `data` must have
dimensionality 1. The axes are `0:0` for the first `Npre`
dimensions, have the axes of `data` for dimension `Npre+1`, and are
`0:0` for the remaining dimensions.

`data` must support `eltype` and `ndims`, but does not have to be an
AbstractArray.

ReshapedOneDs allow one to specify a "filtering dimension" for a
1-dimensional filter.
"""
struct ReshapedOneD{T,N,Npre,V}  # not <: AbstractArray{T,N} (more general, incl. IIR)
    data::V

    function ReshapedOneD{T,N,Npre,V}(data::V) where {T,N,Npre,V}
        ndims(V) == 1 || throw(DimensionMismatch("must be one dimensional, got $(ndims(V))"))
        new{T,N,Npre,V}(data)
    end
end

ReshapedOneD{N,Npre}(data::V) where {N,Npre,V} = ReshapedOneD{eltype(data),N,Npre,V}(data)
# Convenient loop constructor that uses dummy NTuple{N,Bool} to keep
# track of dimensions for type-stability
@inline function ReshapedOneD(pre::NTuple{Npre,Bool}, data, post::NTuple{Npost,Bool}) where {Npre,Npost}
    total = (pre..., true, post...)
    _reshapedvector(total, pre, data)
end
_reshapedvector(::NTuple{N,Bool}, ::NTuple{Npre,Bool}, data) where {N,Npre} = ReshapedOneD{eltype(data),N,Npre,typeof(data)}(data)

# Give ReshapedOneD many of the characteristics of AbstractArray
Base.eltype(A::ReshapedOneD{T}) where {T} = T
Base.ndims(A::ReshapedOneD{_,N}) where {_,N} = N
Base.BroadcastStyle(::Type{ReshapedOneD{T,N,Npre,V}}) where {T,N,Npre,V} = Broadcast.DefaultArrayStyle{N}()
Base.broadcastable(A::ReshapedOneD) = A
Base.isempty(A::ReshapedOneD) = isempty(A.data)

@inline Base.axes(A::ReshapedOneD{_,N,Npre}) where {_,N,Npre} = Base.fill_to_length((Base.ntuple(d->0:0, Val(Npre))..., UnitRange(Base.axes1(A.data))), 0:0, Val(N))

Base.iterate(A::ReshapedOneD) = iterate(A.data)
Base.iterate(A::ReshapedOneD, state) = iterate(A.data, state)

@inline @propagate_inbounds function Base.getindex(A::ReshapedOneD, i::Int)
    @boundscheck checkbounds(A.data, i)
    @inbounds ret = A.data[i]
    ret
end
@inline @propagate_inbounds function Base.getindex(A::ReshapedOneD{T,N,Npre}, I::Vararg{Int,N}) where {T,N,Npre}
    @boundscheck checkbounds_indices(Bool, axes(A), I) || throw_boundserror(A, I)
    @inbounds ret = A.data[I[Npre+1]]
    ret
end
@inline @propagate_inbounds function Base.getindex(A::ReshapedOneD, I::Union{CartesianIndex,Int}...)
    A[Base.IteratorsMD.flatten(I)...]
end

Base.convert(::Type{AA}, A::ReshapedOneD) where {AA<:AbstractArray} = convert(AA, reshape(A.data, axes(A)))

import Base: ==
==(A::ReshapedOneD, B::ReshapedOneD) = convert(AbstractArray, A) == convert(AbstractArray, B)
==(A::ReshapedOneD, B::AbstractArray) = convert(AbstractArray, A) == B
==(A::AbstractArray, B::ReshapedOneD) = A == convert(AbstractArray, B)

import Base: +, -, *, /
for op in (:+, :-, :*, :/)
    @eval begin
        ($op)(A::ReshapedOneD, B::ReshapedOneD) =
            broadcast($op, convert(AbstractArray, A), convert(AbstractArray, B))
        @inline ($op)(As::ReshapedOneD...) =
            broadcast($op, map(A->convert(AbstractArray, A), As)...)
    end
end
Base.BroadcastStyle(::Type{R}) where {R<:ReshapedOneD{T,N}} where {T,N} = Broadcast.DefaultArrayStyle{N}()

_reshape(A::ReshapedOneD{T,N}, ::Val{N}) where {T,N} = A
_vec(A::ReshapedOneD) = A.data
Base.vec(A::ReshapedOneD) = A.data  # is this OK? (note indices won't nec. start with 1)
nextendeddims(a::ReshapedOneD) = 1

"""
    iterdims(inds, v::ReshapedOneD{T,N,Npre}) -> Rpre, ind, Rpost

Return a pair `Rpre`, `Rpost` of CartesianIndicess that correspond to
pre- and post- dimensions for iterating over an array with axes
`inds`. `Rpre` corresponds to the `Npre` "pre" dimensions, and `Rpost`
to the trailing dimensions (not including the vector object wrapped in
`v`).  Concretely,

    Rpre  = CartesianIndices(inds[1:Npre])
    ind   = inds[Npre+1]
    Rpost = CartesianIndices(inds[Npre+2:end])

although the implementation differs for reason of type-stability.
"""
function iterdims(inds::Indices{N}, v::ReshapedOneD{T,N,Npre}) where {T,N,Npre}
    indspre, ind, indspost = _iterdims((), (), inds, v)
    CartesianIndices(indspre), ind, CartesianIndices(indspost)
end
@inline function _iterdims(indspre, ::Tuple{}, inds, v)
    _iterdims((indspre..., inds[1]), (), tail(inds), v)  # consume inds and push to indspre
end
@inline function _iterdims(indspre::NTuple{Npre}, ::Tuple{}, inds, v::ReshapedOneD{<:Any,N,Npre}) where {N,Npre}
    indspre, inds[1], tail(inds)   # return the "central" and trailing dimensions
end

function indexsplit(I::CartesianIndex{N}, v::ReshapedOneD{<:Any,N}) where N
    ipre, i, ipost = _iterdims((), (), Tuple(I), v)
    CartesianIndex(ipre), i, CartesianIndex(ipost)
end

#### FIR filters

## gradients

"""
```julia
    kern1, kern2 = sobel()
```
Return factored  Sobel filters for dimensions 1 and 2 of a two-dimensional
image. Each is a 2-tuple of one-dimensional filters.

# Citation
P.-E. Danielsson and O. Seger, "Generalized and separable sobel operators," in  *Machine Vision for Three-Dimensional Scenes*,  H. Freeman, Ed.  Academic Press, 1990,  pp. 347–379. [doi:10.1016/b978-0-12-266722-0.50016-6](https://doi.org/doi:10.1016/b978-0-12-266722-0.50016-6)

See also: [`Kernel.sobel`](@ref)  and [`ImageFiltering.imgradients`](@ref).
"""
function sobel()
    f1 = centered(SVector( 1.0, 2.0, 1.0)/4)
    f2 = centered(SVector(-1.0, 0.0, 1.0)/2)
    return kernelfactors((f2, f1)), kernelfactors((f1, f2))
end

"""
```julia
    kern = sobel(extended::NTuple{N,Bool}, d)
```
Return a factored Sobel filter for computing the gradient in `N` dimensions
along axis `d`. If `extended[dim]` is false, `kern` will have size 1 along that
dimension.

See also: [`Kernel.sobel`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function sobel(extended::NTuple{N,Bool}, d) where N
    gradfactors(extended, d, [-1, 0, 1]/2, [1, 2, 1]/4)
end

"""
```julia
    kern1, kern2 = prewitt()
```
Return factored Prewitt filters for dimensions 1 and 2 of your image.
Each is a 2-tuple of one-dimensional filters.

# Citation
J. M. Prewitt, "Object enhancement and extraction," *Picture processing and Psychopictorics*, vol. 10, no. 1, pp. 15–19, 1970.

See also: [`Kernel.prewitt`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function prewitt()
    f1 = centered(SVector( 1.0, 1.0, 1.0)/3)
    f2 = centered(SVector(-1.0, 0.0, 1.0)/2)
    return kernelfactors((f2, f1)), kernelfactors((f1, f2))
end

"""
```julia
    kern = prewitt(extended::NTuple{N,Bool}, d)
```
Return a factored Prewitt filter for computing the gradient in `N` dimensions
along axis `d`. If `extended[dim]` is false, `kern` will have size 1 along that
dimension.

See also: [`Kernel.prewitt`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function prewitt(extended::NTuple{N,Bool}, d) where N
    gradfactors(extended, d, [-1, 0, 1]/2, [1, 1, 1]/3)
end

"""
```julia
    kern1, kern2 = scharr()
```
Return factored Scharr filters for dimensions 1 and 2 of your image.  Each is a
2-tuple of one-dimensional filters.

# Citation
H. Scharr and  J. Weickert, "An anisotropic diffusion algorithm with optimized rotation invariance," *Mustererkennung 2000*, pp. 460–467, 2000. [doi:10.1007/978-3-642-59802-9_58](https://doi.org/doi:10.1007/978-3-642-59802-9_58)

See also: [`Kernel.scharr`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function scharr()
    f1 = centered(SVector( 3.0/32.0, 5.0/16.0, 3.0/32.0)*2)
    f2 = centered(SVector(-1.0, 0.0, 1.0)/2)
    return kernelfactors((f2, f1)), kernelfactors((f1, f2))
end

"""
```julia
    kern = scharr(extended::NTuple{N,Bool}, d)
```
Return a factored Scharr filter for computing the gradient in `N` dimensions
along axis `d`. If `extended[dim]` is false, `kern` will have size 1 along that
dimension.

See also: [`Kernel.scharr`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function scharr(extended::NTuple{N,Bool}, d) where N
    # The first factor is the central difference, and since we assume a pixel
    # spacing of one, we divide by 2.
    gradfactors(extended, d, [-1.0, 0.0, 1.0]/2,[3.0/32.0, 5.0/16.0, 3.0/32.0]*2,)
end

"""
```julia
    kern1, kern2 = bickley()
```
Return factored Bickley filters for dimensions 1 and 2 of your image.  Each is
a 2-tuple of one-dimensional filters.

# Citation
W. G. Bickley, "Finite difference formulae for the square lattice," *The Quarterly Journal of Mechanics and Applied Mathematics*, vol. 1, no. 1, pp. 35–42, 1948.  [doi:10.1093/qjmam/1.1.35](https://doi.org/doi:10.1137/12087092x)

See also: [`Kernel.bickley`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function bickley()
    f1 = centered(SVector( 1.0/12.0, 1.0/3.0, 1.0/12.0)*2)
    f2 = centered(SVector(-1.0, 0.0, 1.0)/2)
    return kernelfactors((f2, f1)), kernelfactors((f1, f2))
end

"""
```julia
    kern = bickley(extended::NTuple{N,Bool}, d)
```
Return a factored Bickley filter for computing the gradient in `N` dimensions
along axis `d`. If `extended[dim]` is false, `kern` will have size 1 along that
dimension.

See also: [`Kernel.bickley`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function bickley(extended::NTuple{N,Bool}, d) where N
    # The first factor is the central difference, and since we assume a pixel
    # spacing of one, we divide by 2.
    gradfactors(extended, d, [-1.0, 0.0, 1.0]/2, [1.0/12.0, 1.0/3.0, 1.0/12.0]*2)
end

# Consistent Gradient Operators
# Ando Shigeru
# IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March 2000
#
# TODO: These coefficients were taken from the paper. It would be nice
#       to resolve the optimization problem and use higher precision
#       versions, which might allow better separable approximations of
#       ando4 and ando5.

"""
```julia
    kern1, kern2 = ando3()
```
Return a factored form of Ando's "optimal" ``3 \\times 3`` gradient filters for dimensions 1 and 2 of your image.

# Citation
S. Ando, "Consistent gradient operators," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, no.3, pp. 252–265, 2000. [doi:10.1109/34.841757](https://doi.org/doi:10.1109/34.841757)

See also: [`Kernel.ando3`](@ref),[`KernelFactors.ando4`](@ref),
[`KernelFactors.ando5`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function ando3()
    f1 = centered(2*SVector(0.112737, 0.274526, 0.112737))
    f2 = centered(SVector(-1.0, 0.0, 1.0)/2)
    return kernelfactors((f2, f1)), kernelfactors((f1, f2))
end

"""
```julia
    kern = ando3(extended::NTuple{N,Bool}, d)
```
Return a factored Ando filter (size 3) for computing the gradient in
`N` dimensions along axis `d`.  If `extended[dim]` is false, `kern`
will have size 1 along that dimension.

See also: [`KernelFactors.ando4`](@ref), [`KernelFactors.ando5`](@ref) and
[`ImageFiltering.imgradients`](@ref).
"""
function ando3(extended::NTuple{N,Bool}, d) where N
    gradfactors(extended, d, [-1.0, 0.0, 1.0]/2, 2*[0.112737, 0.274526, 0.112737])
end

"""
```julia
    kern1, kern2 = ando4()
```
Return separable approximations of Ando's "optimal" 4x4 filters for dimensions 1
and 2 of your image.

# Citation
S. Ando, "Consistent gradient operators," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, no.3, pp. 252–265, 2000. [doi:10.1109/34.841757](https://doi.org/doi:10.1109/34.841757)

See also: [`Kernel.ando4`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function ando4()
    f1 = centered(SVector( 0.0919833,  0.408017, 0.408017, 0.0919833))
    f2 = centered(1.46205884*SVector(-0.0919833, -0.408017, 0.408017, 0.0919833))
    return kernelfactors((f2, f1)), kernelfactors((f1, f2))
end

"""
```julia
    kern = ando4(extended::NTuple{N,Bool}, d)
```
Return a factored Ando filter (size 4) for computing the gradient in
`N` dimensions along axis `d`.  If `extended[dim]` is false, `kern`
will have size 1 along that dimension.

# Citation
S. Ando, "Consistent gradient operators," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, no.3, pp. 252–265, 2000. [doi:10.1109/34.841757](https://doi.org/doi:10.1109/34.841757)

See also: [`Kernel.ando4`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function ando4(extended::NTuple{N,Bool}, d) where N
    if N == 2 && all(extended)
        return ando4()[d]
    else
        error("dimensions other than 2 are not yet supported")
    end
end

"""
```julia
    kern1, kern2 = ando5()
```
Return a separable approximations of Ando's "optimal" 5x5 gradient filters for
dimensions 1 and 2 of your image.

# Citation
S. Ando, "Consistent gradient operators," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, no.3, pp. 252–265, 2000. [doi:10.1109/34.841757](https://doi.org/doi:10.1109/34.841757)

See also: [`Kernel.ando5`](@ref) and [`ImageFiltering.imgradients`](@ref).
"""
function ando5()
    f1 = centered(SVector( 0.0357338, 0.248861, 0.43081, 0.248861, 0.0357338))
    f2 = centered(0.784406*SVector(-0.137424, -0.362576, 0.0, 0.362576, 0.137424))
    return kernelfactors((f2, f1)), kernelfactors((f1, f2))
end

"""
```julia
    kern = ando5(extended::NTuple{N,Bool}, d)
```
Return a factored Ando filter (size 5) for computing the gradient in
`N` dimensions along axis `d`.  If `extended[dim]` is false, `kern`
will have size 1 along that dimension.
"""
function ando5(extended::NTuple{N,Bool}, d) where N
    if N == 2 && all(extended)
        return ando5()[d]
    else
        error("dimensions other than 2 are not yet supported")
    end
end

## Gaussian

"""
    gaussian(σ::Real, [l]) -> g

Construct a 1d gaussian kernel `g` with standard deviation `σ`, optionally
providing the kernel length `l`. The default is to extend by two `σ`
in each direction from the center. `l` must be odd.
"""
function gaussian(σ::Real, l = 4*ceil(Int,σ)+1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l>>1
    g = σ == 0 ? [exp(0/(2*oftype(σ, 1)^2))] : [exp(-x^2/(2*σ^2)) for x=-w:w]
    centered(g/sum(g))
end

"""
    gaussian((σ1, σ2, ...), [l]) -> (g1, g2, ...)

Construct a multidimensional gaussian filter as a product of single-dimension
factors, with standard deviation `σd` along dimension `d`. Optionally
provide the kernel length `l`, which must be a tuple of the same
length.
"""
gaussian(σs::NTuple{N,Real}, ls::NTuple{N,Integer}) where {N} =
    kernelfactors( map((σ,l)->gaussian(σ,l), σs, ls) )
gaussian(σs::NTuple{N,Real}) where {N} = kernelfactors(map(σ->gaussian(σ), σs))

gaussian(σs::AbstractVector, ls::AbstractVector) = gaussian((σs...,), (ls...,))
gaussian(σs::AbstractVector) = gaussian((σs...,))

#### IIR

struct TriggsSdika{T,k,l,L} <: IIRFilter{T}
    a::SVector{k,T}
    b::SVector{l,T}
    scale::T
    M::SMatrix{l,k,T,L}
    asum::T
    bsum::T

    TriggsSdika{T,k,l,L}(a, b, scale, M) where {T,k,l,L} =
        new{T,k,l,L}(a, b, scale, M, sum(a), sum(b))
end
"""
    TriggsSdika(a, b, scale, M)

Defines a kernel for one-dimensional infinite impulse response (IIR)
filtering. `a` is a "forward" filter, `b` a "backward" filter, `M` is
a matrix for matching boundary conditions at the right edge, and
`scale` is a constant scaling applied to each element at the
conclusion of filtering.

# Citation

B. Triggs and M. Sdika, "Boundary conditions for Young-van Vliet
recursive filtering". IEEE Trans. on Sig. Proc. 54: 2365-2367
(2006).
"""
TriggsSdika(a::SVector{k,T}, b::SVector{l,T}, scale, M::SMatrix{l,k,T,L}) where {T,k,l,L} = TriggsSdika{T,k,l,L}(a, b, scale, M)

"""
    TriggsSdika(ab, scale)

Create a symmetric Triggs-Sdika filter (with `a = b = ab`). `M` is
calculated for you. Only length 3 filters are currently supported.
"""
function TriggsSdika(a::SVector{3,T}, scale) where T
    a1, a2, a3 = a[1], a[2], a[3]
    Mdenom = (1+a1-a2+a3)*(1-a1-a2-a3)*(1+a2+(a1-a3)*a3)
    M = @SMatrix([-a3*a1+1-a3^2-a2     (a3+a1)*(a2+a3*a1)  a3*(a1+a3*a2);
                  a1+a3*a2            -(a2-1)*(a2+a3*a1)  -(a3*a1+a3^2+a2-1)*a3;
                  a3*a1+a2+a1^2-a2^2   a1*a2+a3*a2^2-a1*a3^2-a3^3-a3*a2+a3  a3*(a1+a3*a2)]);
    TriggsSdika(a, a, scale, M/Mdenom)
end
Base.vec(kernel::TriggsSdika) = kernel
Base.ndims(kernel::TriggsSdika) = 1
Base.ndims(::Type{T}) where {T<:TriggsSdika} = 1
Base.axes1(kernel::TriggsSdika) = 0:0
Base.axes(kernel::TriggsSdika) = (Base.axes1(kernel),)
Base.isempty(kernel::TriggsSdika) = false

iterdims(inds::Indices{1}, kern::TriggsSdika) = (), inds[1], ()
_reshape(kern::TriggsSdika, ::Val{1}) = kern

# Note that there's a sign reversal between Young & Triggs.
"""
    IIRGaussian([T], σ; emit_warning::Bool=true)

Construct an infinite impulse response (IIR) approximation to a
Gaussian of standard deviation `σ`. `σ` may either be a single real
number or a tuple of numbers; in the latter case, a tuple of such filters
will be created, each for filtering a different dimension of an array.

Optionally specify the type `T` for the filter coefficients; if not
supplied, it will match `σ` (unless `σ` is not floating-point, in
which case `Float64` will be chosen).

# Citation

I. T. Young, L. J. van Vliet, and M. van Ginkel, "Recursive Gabor
Filtering". IEEE Trans. Sig. Proc., 50: 2798-2805 (2002).
"""
function IIRGaussian(::Type{T}, σ::Real; emit_warning::Bool = true) where T
    if emit_warning && σ < 1 && σ != 0
        @warn("σ is too small for accuracy")
    end
    m0 = convert(T,1.16680)
    m1 = convert(T,1.10783)
    m2 = convert(T,1.40586)
    q = convert(T,1.31564*(sqrt(1+0.490811*σ*σ) - 1))
    ascale = (m0+q)*(m1*m1 + m2*m2  + 2m1*q + q*q)
    B = (m0*(m1*m1 + m2*m2)/ascale)^2
    # This is what Young et al call -b, but in filt() notation would be called a
    a1 = q*(2*m0*m1 + m1*m1 + m2*m2 + (2*m0+4*m1)*q + 3*q*q)/ascale
    a2 = -q*q*(m0 + 2m1 + 3q)/ascale
    a3 = q*q*q/ascale
    a = SVector(a1,a2,a3)
    TriggsSdika(a, B)
end
IIRGaussian(σ::Real; emit_warning::Bool = true) = IIRGaussian(iirgt(σ), σ; emit_warning=emit_warning)

function IIRGaussian(::Type{T}, σs::NTuple{N,Real}; emit_warning::Bool = true) where {T,N}
    iirg(T, (), σs, tail(ntuple(d->true, Val(N))), emit_warning)
end
IIRGaussian(σs::Tuple; emit_warning::Bool = true) = IIRGaussian(iirgt(σs), σs; emit_warning=emit_warning)

IIRGaussian(σs::AbstractVector; kwargs...) = IIRGaussian((σs...,); kwargs...)
IIRGaussian(::Type{T}, σs::AbstractVector; kwargs...) where {T} = IIRGaussian(T, (σs...,); kwargs...)

iirgt(σ::AbstractFloat) = typeof(σ)
iirgt(σ::Real) = Float64
iirgt(σs::Tuple) = promote_type(map(iirgt, σs)...)

@inline function iirg(::Type{T}, pre, σs, post, emit_warning) where T
    kern = ReshapedOneD(pre, IIRGaussian(T, σs[1]; emit_warning=emit_warning), post)
    (kern, iirg(T, (pre..., post[1]), tail(σs), tail(post), emit_warning)...)
end
iirg(::Type{T}, pre, σs::Tuple{Real}, ::Tuple{}, emit_warning) where {T} =
    (ReshapedOneD(pre, IIRGaussian(T, σs[1]; emit_warning=emit_warning), ()),)

###### Utilities

"""
    kernelfactors(factors::Tuple)

Prepare a factored kernel for filtering. If passed a 2-tuple of
vectors of lengths `m` and `n`, this will return a 2-tuple of
`ReshapedVector`s that are effectively of sizes `m×1` and `1×n`. In
general, each successive `factor` will be reshaped to extend along the
corresponding dimension.

If passed a tuple of general arrays, it is assumed that each is shaped
appropriately along its "leading" dimensions; the dimensionality of each is
"extended" to `N = length(factors)`, appending 1s to the size as needed.
"""
function kernelfactors(factors::NTuple{N,AbstractVector}) where {N}
    total = ntuple(d->true, Val(N))
    _kernelfactors(total, factors)
end

@inline _kernelfactors(total::NTuple{N,Bool}, factors::NTuple{M,Any}) where {N,M} =
    (ReshapedOneD{N,N-M}(factors[1]), _kernelfactors(total, Base.tail(factors))...)
_kernelfactors(::NTuple{N,Bool}, ::Tuple{}) where N = ()

# A variant in which we just need to fill out to N dimensions
kernelfactors(factors::NTuple{N,AbstractArray}) where {N} = map(f->_reshape(f, Val(N)), factors)

function gradfactors(extended::NTuple{N,Bool}, d::Int, k1, k2) where N
    kernelfactors(ntuple(i -> kdim(extended[i], i==d ? k1 : k2), Val(N)))
end

kdim(keep::Bool, k) = keep ? centered(k) : OffsetArray([oneunit(eltype(k))], 0:0)

end


text:
using Plots, LaTeXStrings
pyplot() # gr()
x = LinRange(0,2π,100) # with 100 the number of p
plot(x, sin.(x), c = :black)
savefig("figSin1py.pdf")

max_stars_repo_name: lazarusA/helpfulJulia

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

@inline _adignore_call(f) = f()
@inline _adignore_call_pullback(@nospecialize ΔΩ) = (NoTangent(), NoTangent())
ChainRulesCore.rrule(::typeof(_adignore_call), f) = _adignore_call(f), _adignore_call_pullback

macro _adignore(expr)
    :(_adignore_call(() -> $(esc(expr))))
end

_backing(x::Any) = x
_backing(x::Tangent) = backing(x)
_unpack_tangent(x::Any) = _backing(unthunk(x))

const _ZeroLike = Union{AbstractZero,Nothing}
const _az_tangent(x::AbstractZero) = x
const _az_tangent(::Nothing) = ZeroTangent()  # Still necessary? Return NoTangent() instead?

_notangent_to_zerotangent(x::Any) = x
_notangent_to_zerotangent(x::Union{NoTangent,Nothing}) = ZeroTangent()
_notangent_to_zerotangent(x::Union{Tuple,NamedTuple}) = map(_notangent_to_zerotangent, x)

function _tangent_array(A::AbstractArray{T}) where T
    dA = similar(A, float(T))
    fill!(dA, NaN) # For safety
    return dA
end

function _tangent_array(A::ArrayOfSimilarVectors{T}) where T
    dA = similar(A, Vector{float(T)})
    fill!(flatview(dA), NaN) # For safety
    return dA
end

@inline _keep_zerolike(::Type{T}, x, xs...) where T = T(x, xs...)
@inline _keep_zerolike(::Type{T}, x::_ZeroLike, xs...) where T = _az_tangent(x)
@inline _keep_zerolike(f::F, x, xs...) where F = f(x, xs...)
@inline _keep_zerolike(f::F, x::_ZeroLike, xs...) where F = _az_tangent(x)


text:
using ForwardDiff, OrdinaryDiffEq
using LinearAlgebra: I

include("utils.jl")

struct TestEuclidean{N} <: Manifold end
struct TestEuclideanMetric <: Metric end

@testset "scaled Euclidean metric" begin
    Manifolds.manifold_dimension(::TestEuclidean{N}) where {N} = N
    function Manifolds.local_metric(M::MetricManifold{<:TestEuclidean,<:TestEuclideanMetric}, x)
        return Diagonal(1.0:manifold_dimension(M))
    end

    n = 3
    E = TestEuclidean{n}()
    g = TestEuclideanMetric()
    M = MetricManifold(E, g)
    G = Diagonal(1.0:n)
    invG = inv(G)
    @test manifold_dimension(M) == n
    @test base_manifold(M) === E
    @test metric(M) === g

    @test_throws ErrorException local_metric_jacobian(E, zeros(3))
    @test_throws ErrorException christoffel_symbols_second_jacobian(E, zeros(3))

    for vtype in (Vector, SVector{n}, MVector{n})
        x, v, w = vtype(randn(n)), vtype(randn(n)), vtype(randn(n))

        @test local_metric(M, x) ≈ G
        @test inverse_local_metric(M, x) ≈ invG
        @test det_local_metric(M, x) ≈ *(1.0:n...)
        @test log_local_metric_density(M, x) ≈ sum(log.(1.0:n)) / 2
        @test inner(M, x, v, w) ≈ dot(v, G * w)
        @test norm(M, x, v) ≈ sqrt(dot(v, G * v))

        T = 0:.5:10
        @test exp(M, x, v, T) ≈ [x + t * v for t in T]
        @test geodesic(M, x, v, T) ≈ [x + t * v for t in T]

        @test christoffel_symbols_first(M, x) ≈ zeros(n, n, n)
        @test christoffel_symbols_second(M, x) ≈ zeros(n, n, n)
        @test riemann_tensor(M, x) ≈ zeros(n, n, n, n)
        @test ricci_tensor(M, x) ≈ zeros(n, n)
        @test ricci_curvature(M, x) ≈ 0
        @test gaussian_curvature(M, x) ≈ 0
        @test einstein_tensor(M, x) ≈ zeros(n, n)
    end
end

struct TestSphere{N,T} <: Manifold
    r::T
end

struct TestSphericalMetric <: Metric end

@testset "scaled Sphere metric" begin
    Manifolds.manifold_dimension(::TestSphere{N}) where {N} = N
    function Manifolds.local_metric(M::MetricManifold{<:TestSphere,<:TestSphericalMetric}, x)
        r = base_manifold(M).r
        d = similar(x)
        d[1] = r^2
        d[2] = d[1] * sin(x[1])^2
        return Diagonal(d)
    end
    sph_to_cart(θ, ϕ) = [cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]

    n = 2
    r = 10 * rand()
    θ, ϕ = π * rand(), 2π * rand()
    Sr = TestSphere{n,Float64}(r)
    S = Manifolds.Sphere(n)
    g = TestSphericalMetric()
    M = MetricManifold(Sr, g)

    @test manifold_dimension(M) == n
    @test base_manifold(M) === Sr
    @test metric(M) === g

    for vtype in (Vector, SVector{n}, MVector{n})
        x = vtype([θ, ϕ])
        G = Diagonal(vtype([1, sin(θ)^2])) .* r^2
        invG = Diagonal(vtype([1, 1 / sin(θ)^2])) ./ r^2
        v, w = normalize(randn(n)), normalize(randn(n))

        @test local_metric(M, x) ≈ G
        @test inverse_local_metric(M, x) ≈ invG
        @test det_local_metric(M, x) ≈ r^4 * sin(θ)^2
        @test log_local_metric_density(M, x) ≈ 2log(r) + log(sin(θ))
        @test inner(M, x, v, w) ≈ dot(v, G * w)
        @test norm(M, x, v) ≈ sqrt(dot(v, G * v))

        xcart = sph_to_cart(θ, ϕ)
        vcart = [cos(ϕ)*cos(θ) -sin(ϕ)*sin(θ);
                 sin(ϕ)*cos(θ)  cos(ϕ)*sin(θ);
                 -sin(θ) 0] * v

        if !Sys.iswindows() || Sys.ARCH == :x86_64
            @testset "numerically integrated geodesics for $vtype" begin
                T = 0:.1:1
                @test isapprox([sph_to_cart(yi...) for yi in exp(M, x, v, T)],
                               exp(S, xcart, vcart, T); atol=1e-3, rtol=1e-3)
                @test isapprox([sph_to_cart(yi...) for yi in geodesic(M, x, v, T)],
                               geodesic(S, xcart, vcart, T); atol=1e-3, rtol=1e-3)
            end
        end

        Γ₁ = christoffel_symbols_first(M, x)
        for i=1:n, j=1:n, k=1:n
            if (i,j,k) == (1,2,2) || (i,j,k) == (2,1,2)
                @test Γ₁[i,j,k] ≈ r^2*cos(θ)*sin(θ)
            elseif (i,j,k) == (2,2,1)
                @test Γ₁[i,j,k] ≈ -r^2*cos(θ)*sin(θ)
            else
                @test Γ₁[i,j,k] ≈ 0
            end
        end

        Γ₂ = christoffel_symbols_second(M, x)
        for l=1:n, i=1:n, j=1:n
            if (l,i,j) == (1,2,2)
                @test Γ₂[l,i,j] ≈ -cos(θ)*sin(θ)
            elseif (l,i,j) == (2,1,2) || (l,i,j) == (2,2,1)
                @test Γ₂[l,i,j] ≈ cot(θ)
            else
                @test Γ₂[l,i,j] ≈ 0
            end
        end

        R = riemann_tensor(M, x)
        for l=1:n, i=1:n, j=1:n, k=1:n
            if (l,i,j,k) == (2,1,1,2)
                @test R[l,i,j,k] ≈ -1
            elseif (l,i,j,k) == (2,1,2,1)
                @test R[l,i,j,k] ≈ 1
            elseif (l,i,j,k) == (1,2,1,2)
                @test R[l,i,j,k] ≈ sin(θ)^2
            elseif (l,i,j,k) == (1,2,2,1)
                @test R[l,i,j,k] ≈ -sin(θ)^2
            else
                @test R[l,i,j,k] ≈ 0
            end
        end

        @test ricci_tensor(M, x) ≈ G ./ r^2
        @test ricci_curvature(M, x) ≈ 2 / r^2
        @test gaussian_curvature(M, x) ≈ 1 / r^2
        @test einstein_tensor(M, x) ≈ ricci_tensor(M, x) - gaussian_curvature(M, x)  .* G
    end
end

struct BaseManifold{N} <: Manifold end
struct BaseManifoldMetric{M} <: Metric end
struct DefaultBaseManifoldMetric <: Metric end

@testset "Metric decorator" begin
    Manifolds.manifold_dimension(::BaseManifold{N}) where {N} = N
    Manifolds.inner(::BaseManifold, x, v, w) = 2 * dot(v,w)
    Manifolds.exp!(::BaseManifold, y, x, v) = y .= x + 2 * v
    Manifolds.log!(::BaseManifold, v, x, y) = v .= (y - x) / 2
    Manifolds.project_tangent!(::BaseManifold, w, x, v) = w .= 2 .* v
    Manifolds.project_point!(::BaseManifold, y, x) = (y .= x)
    Manifolds.local_metric(::MetricManifold{BaseManifold{N},BaseManifoldMetric{N}},x) where N = 2*one(x*x')
    Manifolds.exp!(::MetricManifold{BaseManifold{N},BaseManifoldMetric{N}}, y, x, v) where N = exp!(base_manifold(M), y, x, v)
    Manifolds.vector_transport_to!(::BaseManifold, vto, x, v, y, ::ParallelTransport) = (vto .= v)
    Manifolds.is_default_metric(M::BaseManifold,G::DefaultBaseManifoldMetric) = Val(true)
    Manifolds.tangent_orthonormal_basis(M::BaseManifold{N},x,v) where {N} = ( [(Matrix(I, N, N)[:,i]) for i in 1:N], zeros(N))
    Manifolds.projected_distribution(M::BaseManifold, d) = ProjectedPointDistribution(M, d, project_point!, rand(d))
    Manifolds.projected_distribution(M::BaseManifold, d, x) = ProjectedPointDistribution(M, d, project_point!, x)

    function Manifolds.flat!(::BaseManifold, v::FVector{Manifolds.CotangentSpaceType}, x, w::FVector{Manifolds.TangentSpaceType})
        v.data .= 2 .* w.data
        return v
    end
    function Manifolds.sharp!(::BaseManifold, v::FVector{Manifolds.TangentSpaceType}, x, w::FVector{Manifolds.CotangentSpaceType})
        v.data .= w.data ./ 2
        return v
    end

    M = BaseManifold{3}()
    g = BaseManifoldMetric{3}()
    MM = MetricManifold(M, g)
    g2 = DefaultBaseManifoldMetric()
    MM2 = MetricManifold(M,g2)

    @test is_decorator_manifold(MM) == Val{true}()
    @test is_decorator_manifold(MM2) == Val{true}()
    @test is_decorator_manifold(M) == Val{false}()
    @test is_default_metric(MM) == is_default_metric(base_manifold(MM),metric(MM))
    @test is_default_metric(MM2) == is_default_metric(base_manifold(MM2),metric(MM2))
    @test is_default_metric(MM2) == Val{true}()

    @test convert(typeof(MM2),M) == MM2
    @test_throws ErrorException convert(typeof(MM),M)
    x = [0.1 0.2 0.4]
    v = [0.5 0.7 0.11]
    w = [0.13 0.17 0.19]
    y = similar(x)

    # Test fallbacks
    @test_throws ErrorException vee!(M,w,x,v)
    @test_throws ErrorException hat!(M,w,x,v)

    @test inner(M, x, v, w) == 2 * dot(v,w)
    @test inner(MM, x, v, w) === inner(M, x, v, w)
    @test norm(MM, x, v) === norm(M, x, v)
    @test exp(M, x, v) == x + 2 * v
    @test exp(MM2, x, v) == exp(M, x, v)
    @test exp!(MM, y, x, v) === exp!(M, y, x, v)
    @test retract!(MM, y, x, v) === retract!(M, y, x, v)
    @test retract!(MM, y, x, v, 1) === retract!(M, y, x, v, 1)
    # without a definition for the metric from the embedding, no projection possible
    @test_throws ErrorException log!(MM, w, x, y) === project_tangent!(M, w, x, y)
    @test_throws ErrorException project_tangent!(MM, w, x, v) === project_tangent!(M, w, x, v)
    @test_throws ErrorException project_point!(MM, y, x) === project_point!(M, y, x)
    @test_throws ErrorException vector_transport_to!(MM, w, x, v, y) === vector_transport_to!(M, w, x, v, y)
    # without DiffEq, these error
    # @test_throws ErrorException exp(MM,x, v, 1:3)
    # @test_throws ErrorException exp!(MM, y, x, v)
    # these always fall back anyways.
    @test zero_tangent_vector!(MM, v, x) === zero_tangent_vector!(M, v, x)
    @test injectivity_radius(MM, x) === injectivity_radius(M, x)
    @test injectivity_radius(MM) === injectivity_radius(M)
    @test is_manifold_point(MM, x) === is_manifold_point(M, x)
    @test is_tangent_vector(MM, x, v) === is_tangent_vector(M, x, v)

    @test_throws ErrorException local_metric(MM2,x)
    @test_throws ErrorException local_metric_jacobian(MM2,x)
    @test_throws ErrorException christoffel_symbols_second_jacobian(MM2,x)
    # MM falls back to nondefault error
    @test_throws ErrorException projected_distribution(MM,1,x)
    @test_throws ErrorException projected_distribution(MM,1)
    @test_throws ErrorException normal_tvector_distribution(MM,x,0.2)

    @test inner(MM2, x, v, w) === inner(M, x, v, w)
    @test norm(MM2, x, v) === norm(M, x, v)
    @test distance(MM2, x, y) === distance(M, x, y)
    @test exp!(MM2, y, x, v) === exp!(M, y, x, v)
    @test log!(MM2, v, x, y) === log!(M, v, x, y)
    @test retract!(MM2, y, x, v) === retract!(M, y, x, v)
    @test retract!(MM2, y, x, v, 1) === retract!(M, y, x, v, 1)

    @test project_point!(MM2, y, x) === project_point!(M, y, x)
    @test project_tangent!(MM2, w, x, v) === project_tangent!(M, w, x, v)
    @test vector_transport_to!(MM2, w, x, v, y) == vector_transport_to!(M, w, x, v, y)
    @test zero_tangent_vector!(MM2, v, x) === zero_tangent_vector!(M, v, x)
    @test injectivity_radius(MM2, x) === injectivity_radius(M, x)
    @test injectivity_radius(MM2) === injectivity_radius(M)
    @test is_manifold_point(MM2, x) === is_manifold_point(M, x)
    @test is_tangent_vector(MM2, x, v) === is_tangent_vector(M, x, v)

    a = Manifolds.projected_distribution(M, Distributions.MvNormal(zero(zeros(3)), 1.0))
    b = Manifolds.projected_distribution(MM2, Distributions.MvNormal(zero(zeros(3)), 1.0))
    @test isapprox(Matrix(a.d.Σ), Matrix(b.d.Σ))
    @test isapprox(a.d.μ, b.d.μ)
    @test tangent_orthonormal_basis(M,x,v) == tangent_orthonormal_basis(MM2,x,v)
    @test_throws ErrorException tangent_orthonormal_basis(MM,x,v)
    cov = flat(M, x, FVector(TangentSpace, v))
    cow = flat(M, x, FVector(TangentSpace, w))
    @test cov.data ≈ flat(MM, x, FVector(TangentSpace, v)).data
    cotspace = CotangentBundleFibers(M)
    cotspace2 = CotangentBundleFibers(MM)
    @test cov.data ≈ 2 * v
    @test inner(M, x, v, w) ≈ inner(cotspace, x, cov.data, cow.data)
    @test inner(MM, x, v, w) ≈ inner(cotspace, x, cov.data, cow.data)
    @test inner(MM, x, v, w) ≈ inner(cotspace2, x, cov.data, cow.data)
    @test sharp(M, x, cov).data ≈ v
end


text:
mutable struct Buffer
    mDevice::Device
    mCreateInfo::BufferCreateInfo
    mMemory::MemoryChunk
    mStagingBuffer::Ref{Buffer}
    mKeepStagingBuffer::Bool

    mHandle::vk.VkBuffer

    function Buffer(device::Device, createInfo::BufferCreateInfo)
        this = new()
        this.mDevice = device
        this.mCreateInfo = createInfo
        this.mHandle = C_NULL
        this.mKeepStagingBuffer = false
        this.mStagingBuffer = Ref{Buffer}()
        if (handleRef(createInfo)[].size > 0)
            this.mHandle = VkExt.createBuffer(getLogicalDevice(this.mDevice), handleRef(createInfo))
        end
        return this
    end
end

@class Buffer

function handle(this::Buffer)::vk.VkBuffer
    return this.mHandle
end

function destroy(this::Buffer)
    vk.vkDestroyBuffer(getLogicalDevice(this.mDevice), this.mHandle, C_NULL)
    destroy(this.mMemory)
    if isdefined(this.mStagingBuffer, :x)
        destroy(this.mStagingBuffer[])
    end
end

function createBuffer(device::Device, createInfo::BufferCreateInfo)::Buffer
    return Buffer(device, createInfo)
end

function setDataVRAM(this::Buffer, data::Vector, dataType::Type)
    setDataVRAM(this, data, Csize_t(sizeof(dataType) * length(data)))
end

function setDataVRAM(this::Buffer, data::Vector, size::Csize_t)
    # TODO
    #RecordingCommandBuffer::convenienceBufferCheck("Buffer::setDataVRAM()");
    this.initHandle(size)
    if !isdefined(this, :mMemory)
        this.realizeVRAM()
    end

    if isMappable(this.mMemory) # For APUs / integrated GPUs
        mapped = map(this.mMemory)
        memmove(getData(mapped), pointer(data), size)
        unmap(mapped)
    else
        if isdefined(this.mStagingBuffer, :x)
            staging = this.mStagingBuffer[]
        else
            createInfo = copyWithUsage(this.mCreateInfo, vk.VkFlags(vk.VK_BUFFER_USAGE_TRANSFER_SRC_BIT))
            staging = this.mDevice.createBuffer(createInfo)
        end
        staging.setDataRAM(data, size)

        this.copyFrom(staging)

        if this.mKeepStagingBuffer
            this.mStagingBuffer[] = staging
        else
            destroy(staging)
        end
    end
end

function setDataVRAM(this::Buffer, data::Vector, size::Csize_t,
                      cmd::RecordingCommandBuffer)::BufferBarrier
    this.initHandle(size)
    if !isdefined(this, :mMemory)
        this.realizeVRAM()
    end

    if this.mMemory.isMappable() # For APUs / integrated GPUs
        mapped = this.mMemory.map()
        memmove(mapped.getData(), pointer(data), size)
        unmap(mapped)
    else
        return this.executeOnStagingBuffer(staging->begin
            staging.setDataRAM(data, size)
            return this.copyFrom(staging, cmd)
        end)
    end
end

function setDataRAM(this::Buffer, data::Vector, size::Csize_t)
    initHandle(this, size)
    if !isdefined(this, :mMemory)
        realizeRAM(this)
    end

    mapped = map(this.mMemory)
    memmove(getData(mapped), pointer(data), size)
    unmap(mapped)
end

function initHandle(this::Buffer, size::Csize_t)
    if (this.mHandle != C_NULL)
        @assert (size <= handleRef(this.mCreateInfo)[].size) "Buffers in Vulkan cannot be " *
                                                             "enlarged. Create a new one or " *
                                                             "start off with a bigger one."
    else
        this.mCreateInfo = copyWithSize(this.mCreateInfo, size)
        this.mHandle = VkExt.createBuffer(getLogicalDevice(this.mDevice), handleRef(this.mCreateInfo))
    end
end

function pushData(this::Buffer, data::Vector, size::Csize_t)
    cmd = this.mDevice.graphicsQueue().beginCommandBuffer()
    pushData(this, data, size, cmd)
    cmd.endCommandBuffer()
end

function pushData(this::Buffer, data::Vector, size::Csize_t,
                   cmd::RecordingCommandBuffer)
    vk.vkCmdUpdateBuffer(cmd.handle(), this.mHandle, 0, size, pointer(data))
end

function getData(this::Buffer, outData::Vector)
    @assert isdefined(this, :mMemory) "Buffer needs to be realized to get data."

    if this.mMemory.isMappable() # For APUs / integrated GPUs
        mapped = this.mMemory.map()
        memmove(pointer(outData), mapped.getData(), this.mCreateInfo.handleRef()[].size)
        unmap(mapped)
    else
        this.executeOnStagingBuffer(staging->begin
            cmd = this.mDevice.graphicsQueue().beginCommandBuffer()
            region = [vk.VkBufferCopy(
                0, #srcOffset::VkDeviceSize
                0, #dstOffset::VkDeviceSize
                this.mCreateInfo.handleRef()[].size #size::VkDeviceSize
            )]
            vk.vkCmdCopyBuffer(cmd.handle(), this.mHandle, staging.mHandle, 1, pointer(region))
            cmd.endCommandBuffer()
            staging.getData(outData)
        end)
    end
end

function realizeRAM(this::Buffer)
    @assert (this.mHandle != C_NULL) "Cannot realize a Buffer that doesn't have a size yet."
    logicalDevice = getLogicalDevice(this.mDevice)
    req = VkExt.getBufferMemoryRequirements(logicalDevice, this.mHandle)
    this.mMemory = allocate(getSuballocator(this.mDevice), req, RAM)
    @assert (getOffset(this.mMemory) % req.alignment == 0)
    bindToBuffer(this.mMemory, this.mHandle)
end

function realizeVRAM(this::Buffer)
    @assert (this.mHandle != C_NULL) "Cannot realize a Buffer that doesn't have a size yet."
    logicalDevice = getLogicalDevice(this.mDevice)
    req = VkExt.getBufferMemoryRequirements(logicalDevice, this.mHandle)
    this.mMemory = allocate(getSuballocator(this.mDevice), req, VRAM)
    @assert (getOffset(this.mMemory) % req.alignment == 0)
    bindToBuffer(this.mMemory, this.mHandle)
end

function stagesForUsage(usage::vk.VkBufferUsageFlags)::vk.VkPipelineStageFlags
    flags::vk.VkPipelineStageFlags = vk.VkFlags(0)
    check = (use, flag) -> begin
        if ((usage & use) != 0)
            flags |= flag
        end
    end

    check(vk.VK_BUFFER_USAGE_TRANSFER_SRC_BIT,                          vk.VK_PIPELINE_STAGE_TRANSFER_BIT)
    check(vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT,                          vk.VK_PIPELINE_STAGE_TRANSFER_BIT)
    check(vk.VK_BUFFER_USAGE_UNIFORM_TEXEL_BUFFER_BIT,                  vk.VK_PIPELINE_STAGE_ALL_COMMANDS_BIT)
    check(vk.VK_BUFFER_USAGE_STORAGE_TEXEL_BUFFER_BIT,                  vk.VK_PIPELINE_STAGE_ALL_COMMANDS_BIT)
    check(vk.VK_BUFFER_USAGE_UNIFORM_BUFFER_BIT,                        vk.VK_PIPELINE_STAGE_ALL_COMMANDS_BIT)
    check(vk.VK_BUFFER_USAGE_STORAGE_BUFFER_BIT,                        vk.VK_PIPELINE_STAGE_ALL_COMMANDS_BIT)
    check(vk.VK_BUFFER_USAGE_INDEX_BUFFER_BIT,                          vk.VK_PIPELINE_STAGE_VERTEX_INPUT_BIT)
    check(vk.VK_BUFFER_USAGE_VERTEX_BUFFER_BIT,                         vk.VK_PIPELINE_STAGE_VERTEX_INPUT_BIT)
    check(vk.VK_BUFFER_USAGE_INDIRECT_BUFFER_BIT,                       vk.VK_PIPELINE_STAGE_DRAW_INDIRECT_BIT)
    check(vk.VK_BUFFER_USAGE_TRANSFORM_FEEDBACK_BUFFER_BIT_EXT,         vk.VK_PIPELINE_STAGE_TRANSFORM_FEEDBACK_BIT_EXT)
    check(vk.VK_BUFFER_USAGE_TRANSFORM_FEEDBACK_COUNTER_BUFFER_BIT_EXT, vk.VK_PIPELINE_STAGE_TRANSFORM_FEEDBACK_BIT_EXT)
    check(vk.VK_BUFFER_USAGE_CONDITIONAL_RENDERING_BIT_EXT,             vk.VK_PIPELINE_STAGE_CONDITIONAL_RENDERING_BIT_EXT)
    check(vk.VK_BUFFER_USAGE_RAY_TRACING_BIT_NV,                        vk.VK_PIPELINE_STAGE_RAY_TRACING_SHADER_BIT_NV)

    return flags
end

function accessForUsage(usage::vk.VkBufferUsageFlags)::vk.VkAccessFlags
    flags::vk.VkAccessFlags = vk.VkFlags(0)

    check = (use, flag) -> begin
        if ((usage & use) != 0)
            flags |= flag
        end
    end

    check(vk.VK_BUFFER_USAGE_TRANSFER_SRC_BIT,                          vk.VK_ACCESS_TRANSFER_READ_BIT)
    check(vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT,                          vk.VK_ACCESS_TRANSFER_WRITE_BIT)
    check(vk.VK_BUFFER_USAGE_UNIFORM_TEXEL_BUFFER_BIT,                  vk.VK_ACCESS_UNIFORM_READ_BIT)
    check(vk.VK_BUFFER_USAGE_STORAGE_TEXEL_BUFFER_BIT,                  vk.VK_ACCESS_SHADER_READ_BIT |
                                                                        vk.VK_ACCESS_SHADER_WRITE_BIT)
    check(vk.VK_BUFFER_USAGE_UNIFORM_BUFFER_BIT,                        vk.VK_ACCESS_UNIFORM_READ_BIT)
    check(vk.VK_BUFFER_USAGE_STORAGE_BUFFER_BIT,                        vk.VK_ACCESS_SHADER_READ_BIT |
                                                                        vk.VK_ACCESS_SHADER_WRITE_BIT)
    check(vk.VK_BUFFER_USAGE_INDEX_BUFFER_BIT,                          vk.VK_ACCESS_INDEX_READ_BIT)
    check(vk.VK_BUFFER_USAGE_VERTEX_BUFFER_BIT,                         vk.VK_ACCESS_VERTEX_ATTRIBUTE_READ_BIT)
    check(vk.VK_BUFFER_USAGE_INDIRECT_BUFFER_BIT,                       vk.VK_ACCESS_INDIRECT_COMMAND_READ_BIT)
    check(vk.VK_BUFFER_USAGE_TRANSFORM_FEEDBACK_BUFFER_BIT_EXT,         vk.VK_ACCESS_TRANSFORM_FEEDBACK_WRITE_BIT_EXT)
    check(vk.VK_BUFFER_USAGE_TRANSFORM_FEEDBACK_COUNTER_BUFFER_BIT_EXT, vk.VK_ACCESS_TRANSFORM_FEEDBACK_COUNTER_READ_BIT_EXT |
                                                                        vk.VK_ACCESS_TRANSFORM_FEEDBACK_COUNTER_WRITE_BIT_EXT)
    check(vk.VK_BUFFER_USAGE_CONDITIONAL_RENDERING_BIT_EXT,             vk.VK_ACCESS_CONDITIONAL_RENDERING_READ_BIT_EXT)
    check(vk.VK_BUFFER_USAGE_RAY_TRACING_BIT_NV,                        vk.VK_ACCESS_ACCELERATION_STRUCTURE_READ_BIT_NV |
                                                                        vk.VK_ACCESS_ACCELERATION_STRUCTURE_WRITE_BIT_NV)

    return flags
end

function copyFrom(this::Buffer, other::Buffer)
    # TODO
    # RecordingCommandBuffer::convenienceBufferCheck("Buffer::copyFrom()")
    queue = this.mDevice.graphicsQueue()
    cmd = queue.beginCommandBuffer()
    this.copyFrom(other, cmd)
    cmd.endCommandBuffer()
end

function copyFrom(this::Buffer, other::Buffer, cmd::RecordingCommandBuffer)::BufferBarrier
    region = [vk.VkBufferCopy(
        0, #srcOffset::VkDeviceSize
        0, #dstOffset::VkDeviceSize
        min(this.mCreateInfo.handleRef()[].size, other.mCreateInfo.handleRef()[].size) #size::VkDeviceSize
    )]

    vk.vkCmdCopyBuffer(cmd.handle(), other.handle(), this.handle(), 1, pointer(region))
    cmd.attachResource(other)

    barrier = BufferBarrier(cmd, this.handle(),
                          srcStage = vk.VkFlags(vk.VK_PIPELINE_STAGE_TRANSFER_BIT),
                          srcAccessMask = vk.VkFlags(vk.VK_ACCESS_TRANSFER_WRITE_BIT),
                          dstStage = stagesForUsage(this.mCreateInfo.handleRef()[].usage),
                          dstAccessMask = accessForUsage(this.mCreateInfo.handleRef()[].usage))
    barrier.apply()
    return barrier
end

function keepStagingBuffer(this::Buffer, val::Bool = true)
    this.mKeepStagingBuffer = val
end

function executeOnStagingBuffer(this::Buffer, callable)
    staging::Buffer
    if isdefined(this.mStagingBuffer, :x)
        staging = this.mStagingBuffer[]
    else
        createInfo = copyWithUsage(this.mCreateInfo, vk.VkFlags(vk.VK_BUFFER_USAGE_TRANSFER_SRC_BIT |
                                                                vk.VK_BUFFER_USAGE_TRANSFER_DST_BIT))
        staging = this.mDevice.createBuffer(createInfo)
        staging.realizeRAM()
    end

    ret = callable(staging)

    if this.mKeepStagingBuffer
        this.mStagingBuffer[] = staging
    else
        destroy(staging)
    end
    return ret
end


text:
using DRIPs
using Test

@testset "Sims2010.jl" begin
    # Write your own tests here.
    println("Running tests ...")
	    β = 0.9;
		ω = 1.0;
		A = [0.95 0.0; 0.0 0.4];
		Q = [√0.0975 0.0; 0.0 √0.86];
		H = [1.0; 1.0];
	println("Solving Sims (2011) ...")
	    p   = Drip(ω,β,A,Q,H)
	    @test p.ss.err < 1e-4
    println("Solving Sims (2011) with fixed capacity ...")
	    cap = DRIPs.capacity(p)
	    capn= DRIPs.capacity(p; unit = "nat")
	    pf  = Drip(cap,β,A,Q,H;fcap=true)
	    @test pf.ss.err < 1e-4
	    @test pf.ss.Σ_p ≈ p.ss.Σ_p atol = 1e-3
    println("Solving transition dynamics for Sims (2011) ...")
	    S   = DRIPs.Signal([1 0; 0 1],[0 0; 0 0])
	    pt  = Trip(p,S;T=4000)
	    @test pt.err < 1e-4;
		@test pt.con_err < 1e-3;
    println("Checking IRFs ...")
		pirfs     = irfs(p,T = 15)
		psims1    = simulate(p,T = 500, burn = 100, seed = 0)
		psims2    = simulate(p,T = 500, burn = 100, N=100, seed = 0)
		@test psims1.x == psims2.x
		psims     = simulate(p,T = 500, burn = 100, N=1, seed = 0)
		ptss      = Trip(p,p.ss.Σ_1;T = 15)
		ptssirfs  = irfs(ptss,T = 15)
		Sp        = DRIPs.Signal([0 0; 0 0], [1 0; 0 1]);
		ptssirfsp = irfs(p,Sp, T = 30)
		ptssirfsp = irfs(p,Sp, T = 10)
		ptssirfsp = irfs(p,Sp, T = 15)
		ptssirfspn= irfs(p,Sp; reoptimize = false, T = 15)
		@test pirfs.x_hat ≈ ptssirfs.x_hat atol = 1e-3
		@test pirfs.x_hat ≈ ptssirfsp.x_hat atol = 1e-3
		@test pirfs.x_hat ≈ ptssirfspn.x_hat atol = 1e-3
	println("Final tests on internal functions ...")
		x(j) = 0.5^j
		sum  = DRIPs.infinitesum(x)
		@test sum ≈ 2 atol = 1e-4
end


text:
# Load libraries.
using Turing, StatsPlots, Random

# Turn off progress monitor.
Turing.setprogress!(false);

# Set a random seed and use the forward_diff AD mode.
Random.seed!(12345678);

# Define the emission parameter.
y = [
    1.0,
    1.0,
    1.0,
    1.0,
    1.0,
    1.0,
    2.0,
    2.0,
    2.0,
    2.0,
    2.0,
    2.0,
    3.0,
    3.0,
    3.0,
    3.0,
    3.0,
    3.0,
    3.0,
    2.0,
    2.0,
    2.0,
    2.0,
    1.0,
    1.0,
    1.0,
    1.0,
    1.0,
    1.0,
    1.0,
];
N = length(y);
K = 3;

# Plot the data we just made.
plot(y; xlim=(0, 30), ylim=(-1, 5), size=(500, 250))

# Turing model definition.
@model function BayesHmm(y, K)
    # Get observation length.
    N = length(y)

    # State sequence.
    s = tzeros(Int, N)

    # Emission matrix.
    m = Vector(undef, K)

    # Transition matrix.
    T = Vector{Vector}(undef, K)

    # Assign distributions to each element
    # of the transition matrix and the
    # emission matrix.
    for i in 1:K
        T[i] ~ Dirichlet(ones(K) / K)
        m[i] ~ Normal(i, 0.5)
    end

    # Observe each point of the input.
    s[1] ~ Categorical(K)
    y[1] ~ Normal(m[s[1]], 0.1)

    for i in 2:N
        s[i] ~ Categorical(vec(T[s[i - 1]]))
        y[i] ~ Normal(m[s[i]], 0.1)
    end
end;

g = Gibbs(HMC(0.01, 50, :m, :T), PG(120, :s))
chn = sample(BayesHmm(y, 3), g, 1000);

# Extract our m and s parameters from the chain.
m_set = MCMCChains.group(chn, :m).value
s_set = MCMCChains.group(chn, :s).value

# Iterate through the MCMC samples.
Ns = 1:length(chn)

# Make an animation.
animation = @gif for i in Ns
    m = m_set[i, :]
    s = Int.(s_set[i, :])
    emissions = m[s]

    p = plot(
        y;
        chn=:red,
        size=(500, 250),
        xlabel="Time",
        ylabel="State",
        legend=:topright,
        label="True data",
        xlim=(0, 30),
        ylim=(-1, 5),
    )
    plot!(emissions; color=:blue, label="Sample $i")
end every 3

# Index the chain with the persistence probabilities.
subchain = chn[["T[1][1]", "T[2][2]", "T[3][3]"]]

plot(subchain; seriestype=:traceplot, title="Persistence Probability", legend=false)

heideldiag(MCMCChains.group(chn, :T))[1]

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
module TestJSONWiringDiagrams

using Test
import JSON
using Catlab.Doctrines
using Catlab.WiringDiagrams

# Round trip wiring diagrams with dictionary box and wire data.

function roundtrip(f::WiringDiagram)
  json = sprint(JSON.print, generate_json_graph(f), 2)
  parse_json_graph(Dict, Nothing, Dict, JSON.parse(json))
end

ports(n) = Nothing[ nothing for i in 1:n ]
diagram = WiringDiagram(ports(1), ports(1))
f = Box(Dict("name" => "f", "type" => "foo"), ports(1), ports(1))
fv = add_box!(diagram, f)
add_wires!(diagram, [
  Wire(Dict("name" => "w1", "type" => "woo"),
      (input_id(diagram), 1) => (fv, 1)),
  Wire(Dict("name" => "w2", "type" => "woo"),
       (fv, 1) => (output_id(diagram), 1)),
])
@test roundtrip(diagram) == diagram

# Round trip wiring diagrams with symbolic box and port values.

function roundtrip_symbolic(f::WiringDiagram)
  json = sprint(JSON.print, generate_json_graph(f), 2)
  parse_json_graph(Symbol, Symbol, Nothing, JSON.parse(json))
end

A, B, C = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C)
f = to_wiring_diagram(Hom(:f, A, B))
g = to_wiring_diagram(Hom(:g, B, C))

@test roundtrip_symbolic(f) == f
@test roundtrip_symbolic(compose(f,g)) == compose(f,g)
@test roundtrip_symbolic(otimes(f,g)) == otimes(f,g)

# Round trip nested, symbolic wiring diagrams.

f = to_wiring_diagram(Hom(:f, A, B))
diagram = WiringDiagram(:outer, [:A], [:B])
fv = add_box!(diagram, f)
add_wires!(diagram, [
  Wire((input_id(diagram),1) => (fv,1)),
  Wire((fv,1) => (output_id(diagram),1))
])
@test roundtrip_symbolic(diagram) == diagram

end


text:
arg_dict = Dict{Any, Any}(
# specify the physics and SBP operator
"physics" => "Euler",  # specify physics to run
"dimensions" => 2,  # this is a two dimensional mesh

# specify temporal and spatial discretization
"run_type" => 30,  # use LSERK 54
"real_time" => true,  # this is an unsteady problem
"t_max" => 10.0,  # make time
"operator_type" => "SBPOmega",  # specify SBP operator
"order" => 1,  # p = 1 operator
"use_DG" => true,  # use discontinuous galerkin solver
"Flux_name" => "RoeFlux",  # numerical flux function used in face integrals
"CFL" => 0.10,  # CFL number

# specify the problem itself
"IC_name" => "ICUnsteadyVortex",  # initial condtiion
"vortex_x0" => 5.0,  # initial location of vortex
"numBC" => 1,  # number of boundary conditions
"BC1" => [0, 1, 2, 3],  # geometric edges to apply the BC to
"BC1_name" => "unsteadyVortexBC",   # name of boundary condition

# specify mesh
"smb_name" => "../../meshes/vortex_s_1_.smb",
"dmg_name" => ".null",

# misc options
"write_vis" => true,  # write paraview files
"output_freq" => 100,  # do output every this many time steps
"do_postproc" => true,  # calculate error at end of run
"exact_soln_func" => "ICUnsteadyVortex",  # use this function for the exact soltuion (to calculate error)
)


text: min_dist = 0.8;  # minimum separation distance required [m]


text:
function print_exit_message(time_taken::Float64)

        println("Finished!")
        println(" ")
        println("Thank you for choosing KG Solver!")
        println("___________________________________________________________")
        print("Time Taken:  " , time_taken )

        return
end


text:
## ------------------------------------------------------------------
const GLOB_STATE = Dict{Symbol, Any}()

currobj!(o::Union{Nothing, RBObject}) = (setindex!(GLOB_STATE, o, :currobj); o)
currobj!() = currobj!(nothing)
currobj() = get!(GLOB_STATE, :currobj, nothing)

currsec!(s::Union{Nothing, RBSection}) = (setindex!(GLOB_STATE, currobj!(s), :currsec); s)
currsec!() = (currsec!(nothing); currobj!())
currsec() = get!(GLOB_STATE, :currsec, nothing)

currdoc!(d::Union{Nothing, RBDoc}) = (setindex!(GLOB_STATE, currobj!(d), :currdoc); d)
currdoc!() = (currdoc!(nothing); currsec!())
currdoc() = get!(GLOB_STATE, :currdoc, nothing)

currbook!(b::Union{Nothing, RBook}) = (setindex!(GLOB_STATE, currobj!(b), :currbook); b)
currbook!() = (currbook!(nothing); currdoc!())
currbook() = get!(GLOB_STATE, :currbook, nothing)
clearbook!() = currbook!()

## ------------------------------------------------------------------
# book
function _check_currbook()
    book = currbook()
    isnothing(book) && error("No `RBook` selected. See `@openbook` help.")
    return book
end

function _check_currdoc()
    doc = currdoc()
    isnothing(doc) && error("No `RBDoc` selected. See `@new_document!` help.")
    return doc
end

function _check_currsec()
    sec = currsec()
    isnothing(sec) && error("No `RBSec` selected. See `@new_section!` help.")
    return sec
end

function _check_currobj()
    obj = currobj()
    isnothing(obj) && error("No `RBObject` selected. See `currobj` help.")
    return obj
end

## ------------------------------------------------------------------
# Functional Interface
bookbib() = bookbib(bookdir())


text:
"""
  $(TYPEDEF)
  Abstract type for a factorization   with ExtandableSparseMatrix.

  Any such preconditioner should have the following fields
````
  A::ExtendableSparseMatrix
  fact
  phash::UInt64
````
and  provide a method
````
  update!(precon)
````

  The idea is that, depending if the matrix pattern has changed,
  different steps are needed to update the preconditioner.

  Moreover, they have the ExtendableSparseMatrix as a field, ensuring
  consistency after construction.
"""
abstract type AbstractFactorization{Tv, Ti} end

"""
$(TYPEDEF)

Abstract subtype for preconditioners
"""
abstract type AbstractPreconditioner{Tv, Ti} <:AbstractFactorization{Tv, Ti} end

"""
$(TYPEDEF)

Abstract subtype for (full) LU factorizations
"""
abstract type AbstractLUFactorization{Tv, Ti} <:AbstractFactorization{Tv, Ti}  end

"""
```
issolver(factorization)
```

Determine if factorization is a solver or not
"""
issolver(::AbstractLUFactorization)=true
issolver(::AbstractPreconditioner)=false

"""
```
factorize!(factorization, matrix)
```

Update or create factorization, possibly reusing information from the current state.
This method is aware of pattern changes.
"""
function factorize!(p::AbstractFactorization, A::ExtendableSparseMatrix)
    p.A=A
    update!(p)
    p
end

"""
```
lu!(factorization, matrix)
```

Update LU factorization, possibly reusing information from the current state.
This method is aware of pattern changes.

If `nothing` is passed as first parameter, [`factorize`](@ref) is called.
"""
LinearAlgebra.lu!(lufact::AbstractFactorization, A::ExtendableSparseMatrix)=factorize!(lufact,A)

LinearAlgebra.lu(A::ExtendableSparseMatrix)=factorize!(LUFactorization(),A)

"""
```
 lufact\rhs
```

Solve  LU factorization problem.
"""
Base.:\(lufact::AbstractLUFactorization, v::AbstractArray{T,1} where T)=ldiv!(similar(v), lufact,v)

"""
```
update!(factorization)
```
Update factorization after matrix update.
"""
update!(::AbstractFactorization)

"""
```
ldiv!(u,factorization,v)
ldiv!(factorization,v)
```

Solve factorization.
"""
LinearAlgebra.ldiv!(u,fact::AbstractFactorization, v)=ldiv!(u, fact.fact, v)
LinearAlgebra.ldiv!(fact::AbstractFactorization, v)=ldiv!(fact.fact,v)

macro makefrommatrix(fact)
    return quote
        $(esc(fact))(A::ExtendableSparseMatrix{Tv,Ti}; kwargs...) where {Tv,Ti} = factorize!($(esc(fact))(;valuetype=Tv, indextype=Ti,  kwargs...),A)
        $(esc(fact))(A::SparseMatrixCSC{Tv,Ti}; kwargs...) where {Tv,Ti} = $(esc(fact))(ExtendableSparseMatrix(A); kwargs...)
    end
end

include("jacobi.jl")
include("ilu0.jl")
include("parallel_jacobi.jl")
include("umfpack_lu.jl")
include("cholmod_cholesky.jl")

@eval begin

    @makefrommatrix LUFactorization
    @makefrommatrix CholeskyFactorization
    @makefrommatrix ILU0Preconditioner
    @makefrommatrix JacobiPreconditioner
    @makefrommatrix ParallelJacobiPreconditioner

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# displaying type warnings

function warntype_type_printer(io::IO, @nospecialize(ty), used::Bool)
    used || return
    if ty isa Type && (!Base.isdispatchelem(ty) || ty == Core.Box)
        if ty isa Union && Base.is_expected_union(ty)
            Base.emphasize(io, "::$ty", Base.warn_color()) # more mild user notification
        else
            Base.emphasize(io, "::$ty")
        end
    else
        Base.printstyled(io, "::$ty", color=:cyan) # show the "good" type
    end
    nothing
end

"""
    code_warntype([io::IO], f, types; debuginfo=:default)

Prints lowered and type-inferred ASTs for the methods matching the given generic function
and type signature to `io` which defaults to `stdout`. The ASTs are annotated in such a way
as to cause "non-leaf" types to be emphasized (if color is available, displayed in red).
This serves as a warning of potential type instability. Not all non-leaf types are particularly
problematic for performance, so the results need to be used judiciously.
In particular, unions containing either [`missing`](@ref) or [`nothing`](@ref) are displayed in yellow, since
these are often intentional.

Keyword argument `debuginfo` may be one of `:source` or `:none` (default), to specify the verbosity of code comments.

See [`@code_warntype`](@ref man-code-warntype) for more information.
"""
function code_warntype(io::IO, @nospecialize(f), @nospecialize(t); debuginfo::Symbol=:default, optimize::Bool=false)
    debuginfo = Base.IRShow.debuginfo(debuginfo)
    lineprinter = Base.IRShow.__debuginfo[debuginfo]
    for (src, rettype) in code_typed(f, t, optimize=optimize)
        lambda_io::IOContext = io
        if src.slotnames !== nothing
            slotnames = Base.sourceinfo_slotnames(src)
            lambda_io = IOContext(lambda_io, :SOURCE_SLOTNAMES => slotnames)
            println(io, "Variables")
            slottypes = src.slottypes
            for i = 1:length(slotnames)
                print(io, "  ", slotnames[i])
                if isa(slottypes, Vector{Any})
                    warntype_type_printer(io, slottypes[i], true)
                end
                println(io)
            end
            println(io)
        end
        print(io, "Body")
        warntype_type_printer(io, rettype, true)
        println(io)
        # TODO: static parameter values
        Base.IRShow.show_ir(lambda_io, src, lineprinter(src), warntype_type_printer)
    end
    nothing
end
code_warntype(@nospecialize(f), @nospecialize(t); kwargs...) =
    code_warntype(stdout, f, t; kwargs...)

import Base.CodegenParams

# Printing code representations in IR and assembly
function _dump_function(@nospecialize(f), @nospecialize(t), native::Bool, wrapper::Bool,
                        strip_ir_metadata::Bool, dump_module::Bool, syntax::Symbol,
                        optimize::Bool, debuginfo::Symbol=:default,
                        params::CodegenParams=CodegenParams())
    ccall(:jl_is_in_pure_context, Bool, ()) && error("code reflection cannot be used from generated functions")
    if isa(f, Core.Builtin)
        throw(ArgumentError("argument is not a generic function"))
    end
    # get the MethodInstance for the method match
    world = typemax(UInt)
    meth = which(f, t)
    t = to_tuple_type(t)
    tt = signature_type(f, t)
    (ti, env) = ccall(:jl_type_intersection_with_env, Any, (Any, Any), tt, meth.sig)::Core.SimpleVector
    meth = Base.func_for_method_checked(meth, ti, env)
    linfo = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance}, (Any, Any, Any, UInt), meth, ti, env, world)
    # get the code for it
    return _dump_function_linfo(linfo, world, native, wrapper, strip_ir_metadata, dump_module, syntax, optimize, debuginfo, params)
end

function _dump_function_linfo(linfo::Core.MethodInstance, world::UInt, native::Bool, wrapper::Bool,
                              strip_ir_metadata::Bool, dump_module::Bool, syntax::Symbol,
                              optimize::Bool, debuginfo::Symbol=:default,
                              params::CodegenParams=CodegenParams())
    if syntax !== :att && syntax !== :intel
        throw(ArgumentError("'syntax' must be either :intel or :att"))
    end
    if debuginfo === :default
        debuginfo = :source
    elseif debuginfo !== :source && debuginfo !== :none
        throw(ArgumentError("'debuginfo' must be either :source or :none"))
    end
    if native
        llvmf = ccall(:jl_get_llvmf_decl, Ptr{Cvoid}, (Any, UInt, Bool, CodegenParams), linfo, world, wrapper, params)
    else
        llvmf = ccall(:jl_get_llvmf_defn, Ptr{Cvoid}, (Any, UInt, Bool, Bool, CodegenParams), linfo, world, wrapper, optimize, params)
    end
    if llvmf == C_NULL
        error("could not compile the specified method")
    end

    if native
        str = ccall(:jl_dump_function_asm, Ref{String},
                    (Ptr{Cvoid}, Cint, Ptr{UInt8}, Ptr{UInt8}), llvmf, 0, syntax, debuginfo)
    else
        str = ccall(:jl_dump_function_ir, Ref{String},
                    (Ptr{Cvoid}, Bool, Bool, Ptr{UInt8}), llvmf, strip_ir_metadata, dump_module, debuginfo)
    end

    # TODO: use jl_is_cacheable_sig instead of isdispatchtuple
    isdispatchtuple(linfo.specTypes) || (str = "; WARNING: This code may not match what actually runs.\n" * str)
    return str
end

"""
    code_llvm([io=stdout,], f, types; raw=false, dump_module=false, optimize=true, debuginfo=:default)

Prints the LLVM bitcodes generated for running the method matching the given generic
function and type signature to `io`.

If the `optimize` keyword is unset, the code will be shown before LLVM optimizations.
All metadata and dbg.* calls are removed from the printed bitcode. For the full IR, set the `raw` keyword to true.
To dump the entire module that encapsulates the function (with declarations), set the `dump_module` keyword to true.
Keyword argument `debuginfo` may be one of source (default) or none, to specify the verbosity of code comments.
"""
code_llvm(io::IO, @nospecialize(f), @nospecialize(types), raw::Bool, dump_module::Bool=false, optimize::Bool=true, debuginfo::Symbol=:default) =
    print(io, _dump_function(f, types, false, false, !raw, dump_module, :att, optimize, debuginfo))
code_llvm(io::IO, @nospecialize(f), @nospecialize(types=Tuple); raw::Bool=false, dump_module::Bool=false, optimize::Bool=true, debuginfo::Symbol=:default) =
    code_llvm(io, f, types, raw, dump_module, optimize, debuginfo)
code_llvm(@nospecialize(f), @nospecialize(types=Tuple); raw=false, dump_module=false, optimize=true, debuginfo::Symbol=:default) =
    code_llvm(stdout, f, types; raw=raw, dump_module=dump_module, optimize=optimize, debuginfo=debuginfo)

"""
    code_native([io=stdout,], f, types; syntax=:att, debuginfo=:default)

Prints the native assembly instructions generated for running the method matching the given
generic function and type signature to `io`.
Switch assembly syntax using `syntax` symbol parameter set to `:att` for AT&T syntax or `:intel` for Intel syntax.
Keyword argument `debuginfo` may be one of source (default) or none, to specify the verbosity of code comments.
"""
code_native(io::IO, @nospecialize(f), @nospecialize(types=Tuple); syntax::Symbol=:att, debuginfo::Symbol=:default) =
    print(io, _dump_function(f, types, true, false, false, false, syntax, true, debuginfo))
code_native(@nospecialize(f), @nospecialize(types=Tuple); syntax::Symbol=:att, debuginfo::Symbol=:default) =
    code_native(stdout, f, types; syntax=syntax, debuginfo=debuginfo)
code_native(::IO, ::Any, ::Symbol) = error("illegal code_native call") # resolve ambiguous call

max_stars_repo_name: Dictino/julia

max_issues_repo_name: Dictino/julia

max_forks_repo_name: Dictino/julia