text:
# calculateposterior.jl

"""
    calculatedependentfaberschauderfunctions(higestlevel::Int64)

Internal function, not exported.

Calculates which Faber-Schauder functions (ψ_i,ψ_j) have essential overlapping
support, i<j. All other combinations (ψ_i, ψ_j), i<j, have essential
nonoverlapping support and hence int_0^T ψ_i(X_t)ψ_j(X_t)dt=0, which we use in
calculategirsanovmatrix.

Returns a triple (lengthvectors, rowindices, columnindices) where
```julia
lengthvectors == length(rowindices) == length(columnindices),
```
and for i in 1:lengthvectors (ψ_rowindices[i],ψ_columnindices[i]) have
essentially overlapping support. For all i, rowindices[i] < columnindices[i].
"""
function calculatedependentfaberschauderfunctions(higestlevel::Int64)
    lengthvectors = higestlevel*2^(higestlevel+1)+1
    rowindices = Vector{Int}(lengthvectors)
    columnindices = Vector{Int}(lengthvectors)
    index = 1
    # column value need to be strictly greater than row value.
    # psi_1 is dependent with every Faber-Schauder function.
    for i in 2:2^(higestlevel+1)
        rowindices[index] = 1
        columnindices[index] = i
        index += 1
    end
    # psi_{j,k}=psi_{2^j+k} is dependent with
    # psi_{j+d,(k-1)2^d+1},...,psi_{j+d, k2^d}, d\ge1.
    for jone in 0:higestlevel-1
        twotothepowerjone = 2^jone
        for kone in 1:2^jone
            ione = twotothepowerjone+kone
            for  jtwo in jone+1:higestlevel
                d = jtwo - jone
                twotothepowerjtwo = 2^jtwo
                for ktwo in (kone-1)*2^d+1:kone*2^d
                    rowindices[index] = ione
                    columnindices[index] = twotothepowerjtwo+ktwo
                    index += 1
                end
            end
        end
    end
    return (lengthvectors, rowindices, columnindices)
end
#test written

# calculategirsanovmatrixelement methods:

"""
    calculategirsanovmatrixelement(ψ1Xt, ψ2Xt, σ, Δt)
    calculategirsanovmatrixelement(samplevalueindices::BitArray{1},
        ψ1Xt, ψ2Xt, σ, Δt)

Internal function, not exported!

Where ψ1Xt, ψ2Xt are AbstractArray{Float64}, and σ and Δt are either an
AbstractArray{Float64} or a Float64.

Let ψ_a and ψ_b denote two basis elements.
calculategirsanovmatrixelement calculates the (a,b) element of the Girsanov
matrix defined by int_0^T ψ_a(X_t)ψ_b(X_t)/(σ^2(X_t)) dt, where ψ_a is given
by ψ1Xt (already evaluated in X) and ψ_b by ψ2Xt (already evaluated in X).
"""
function calculategirsanovmatrixelement(
        ψ1Xt::S,
        ψ2Xt::T,
        σ::Float64,
        Δt::Float64) where {S<:AbstractArray{Float64},
            T<:AbstractArray{Float64}}
    return Δt * sum(ψ1Xt .* ψ2Xt) / (σ * σ)
end
#test written

function calculategirsanovmatrixelement(
        ψ1Xt::R,
        ψ2Xt::S,
        σ::Float64,
        Δt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{Float64},
            T<:AbstractArray{Float64}}
    return sum(Δt .* ψ1Xt .* ψ2Xt) / (σ * σ)
end
#test written

function calculategirsanovmatrixelement(
        ψ1Xt::R,
        ψ2Xt::S,
        σXt::T,
        Δt::Float64) where {R<:AbstractArray{Float64},
            S<:AbstractArray{Float64}, T<:AbstractArray{Float64}}
    return Δt * sum(ψ1Xt .* ψ2Xt ./ (σXt .* σXt))
end
#test written

function calculategirsanovmatrixelement(
        ψ1Xt::Q,
        ψ2Xt::R,
        σXt::S,
        Δt::T) where {Q<:AbstractArray{Float64}, R<:AbstractArray{Float64},
            S<:AbstractArray{Float64}, T<:AbstractArray{Float64}}
    return sum(Δt .* ψ1Xt .* ψ2Xt ./ (σXt .* σXt))
end
#test written

function calculategirsanovmatrixelement(
        samplevalueindices::BitArray{1},
        ψ1Xt::S,
        ψ2Xt::T,
        σ::Float64,
        Δt::Float64) where {S<:AbstractArray{Float64}, T<:AbstractArray{Float64}}
    filteredψ1Xt = ψ1Xt[samplevalueindices]
    filteredψ2Xt = ψ2Xt[samplevalueindices]
    return Δt * sum(filteredψ1Xt .* filteredψ2Xt) / (σ * σ)
end
#test written

function calculategirsanovmatrixelement(
        samplevalueindices::BitArray{1},
        ψ1Xt::R,
        ψ2Xt::S,
        σ::Float64,
        Δt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{Float64},
            T<:AbstractArray{Float64}}
    filteredΔt = Δt[samplevalueindices]
    filteredψ1Xt = ψ1Xt[samplevalueindices]
    filteredψ2Xt = ψ2Xt[samplevalueindices]
    return sum(filteredΔt .* filteredψ1Xt .* filteredψ2Xt) / (σ * σ)
end
#test written

function calculategirsanovmatrixelement(
        samplevalueindices::BitArray{1},
        ψ1Xt::R,
        ψ2Xt::S,
        σXt::T,
        Δt::Float64) where {R<:AbstractArray{Float64}, S<:AbstractArray{Float64},
            T<:AbstractArray{Float64}}
    filteredψ1Xt = ψ1Xt[samplevalueindices]
    filteredψ2Xt = ψ2Xt[samplevalueindices]
    filteredσXt = σXt[samplevalueindices]
    return Δt * sum(filteredψ1Xt .* filteredψ2Xt ./
        (filteredσXt .* filteredσXt))
end
#test written

function calculategirsanovmatrixelement(
        samplevalueindices::BitArray{1},
        ψ1Xt::Q,
        ψ2Xt::R,
        σXt::S,
        Δt::T) where {Q<:AbstractArray{Float64}, R<:AbstractArray{Float64},
            S<:AbstractArray{Float64}, T<:AbstractArray{Float64}}
    filteredψ1Xt = ψ1Xt[samplevalueindices]
    filteredψ2Xt = ψ2Xt[samplevalueindices]
    filteredσXt = σXt[samplevalueindices]
    filteredΔt = Δt[samplevalueindices]
    return sum(filteredΔt .* filteredψ1Xt .* filteredψ2Xt ./ (
        filteredσXt .* filteredσXt))
end
#test written

"""
    calculateΔt(timeinterval::Range{Float64}) = step(timeinterval)
    calculateΔt(timeinterval::AbstractArray{Float64}) = timeinterval[2:end] -
        timeinterval[1:end-1]

Internal function, not exported!

Calculated Δt, returns a Float64, when timeinterval is a range object, and
wise the increments of the timeinterval vector.
"""
calculateΔt(timeinterval::Range{Float64}) = step(timeinterval)
calculateΔt(timeinterval::AbstractArray{Float64}) = timeinterval[2:end] -
    timeinterval[1:end-1]
# Test written.

# calculategirsanovmatrix methods:

"""
    calculategirsanovmatrix(
        Π::FaberSchauderExpansionWithGaussianCoefficients,
        samplevalueindices::Vector{BitArray{1}},
        X::SamplePath,
        ψXt::Vector{Array{Float64,1}},
        σXt::T) where {T<:Union{S, Float64} where {S<:AbstractArray{Float64}}}
    calculategirsanovmatrix(
        Π::AbstractGaussianProcess,
        X::SamplePath,
        ψXt::S,
        σXt::T) where {S<:AbstractArray{U} where {U<:AbstractArray{Float64}},
            T<:Union{Float64, V} where {V<:AbstractArray{Float64}}}

Internal function, not exported!

Calculates the Girsanov matrix. In case of a
FaberSchauderExpansionWithGaussianCoefficients-object it will make use of the
sparsity structure of the Faber-Schauder basis.
"""
function calculategirsanovmatrix(
        sizesquarematrix::Int,
        timeinterval::R,
        ψXt::S,
        σXt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{U} where U<:AbstractArray{Float64}, T<:Union{Float64, V} where V<:AbstractArray{Float64}}
    girsanovmatrix = Array{Float64}(sizesquarematrix, sizesquarematrix)
    Δt = calculateΔt(timeinterval)
    for k in 1:sizesquarematrix
        for l in k:sizesquarematrix
            girsanovmatrix[k,l] = calculategirsanovmatrixelement(
                ψXt[k], ψXt[l], σXt, Δt)
        end
    end
    return Symmetric(girsanovmatrix)
end

function calculategirsanovmatrix(
        Π::FaberSchauderExpansionWithGaussianCoefficients,
        samplevalueindices::Vector{BitArray{1}},
        timeinterval::S,
        ψXt::Vector{Array{Float64,1}},
        σXt::T) where {S<:AbstractArray{Float64},
            T<:Union{U, Float64} where {U<:AbstractArray{Float64}}}
    lengthvectors, rowindices, columnindices =
        calculatedependentfaberschauderfunctions(Π.higestlevel)
    d = length(Π)
    numberofnonzeroelements = (2*Π.higestlevel+1)*d+2
    V = Vector{Float64}(numberofnonzeroelements)
    Δt = calculateΔt(timeinterval)
    for i in 1:d
        V[i] = calculategirsanovmatrixelement(samplevalueindices[i],
            ψXt[i], ψXt[i], σXt, Δt)
    end
    rowindices, columnindices = vcat(1:d, rowindices, columnindices),
        vcat(1:d, columnindices, rowindices)
    for i in d+1:d+lengthvectors
        V[i+lengthvectors] = V[i] = calculategirsanovmatrixelement(
            samplevalueindices[columnindices[i]], ψXt[rowindices[i]],
            ψXt[columnindices[i]], σXt, Δt)
    end
    return sparse(rowindices, columnindices, V, d, d)
end

"""
    calculategirsanovvectorelement(
        samplevalueindices::BitArray{1},
        ΔXt::R,
        ψXt::S,
        σXt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{Float64},
            T<:AbstractArray{Float64}}
    calculategirsanovvectorelement(
        samplevalueindices::BitArray{1},
        ΔXt::S,
        ψXt::T,
        σ::Float64) where {S<:AbstractArray{Float64}, T<:AbstractArray{Float64}}
    calculategirsanovvectorelement(
        ΔXt::R,
        ψXt::S,
        σXt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{Float64},
            T<:AbstractArray{Float64}}
    calculategirsanovvectorelement(
        ΔXt::S,
        ψXt::T,
        σ::Float64) where {S<:AbstractArray{Float64}, T<:AbstractArray{Float64}}

Internal function, used by calculategirsanovvector, not exported!

Calculates the kth Girsanov vector element int_0^T ψ_k(X_t)/(σ^2(X_t)) dX_t,
where ψ_k is the kth basis function and T is the end time.
"""
function calculategirsanovvectorelement(
    ΔXt::R,
    ψXt::S,
    σXt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{Float64},
        T<:AbstractArray{Float64}}
    return sum(ΔXt .* ψXt ./ (σXt .* σXt))
end
#Test written.

function calculategirsanovvectorelement(
    ΔXt::S,
    ψXt::T,
    σ::Float64) where {S<:AbstractArray{Float64}, T<:AbstractArray{Float64}}
    return sum(ΔXt .* ψXt) / (σ * σ)
end
#Test written.

function calculategirsanovvectorelement(
        samplevalueindices::BitArray{1},
        ΔXt::S,
        ψXt::T,
        σ::Float64) where {S<:AbstractArray{Float64}, T<:AbstractArray{Float64}}
    filteredΔXt = ΔXt[samplevalueindices]
    filteredψXt = ψXt[samplevalueindices]
    return sum(filteredΔXt .* filteredψXt) / (σ*σ)
end
#Test written.

function calculategirsanovvectorelement(
        samplevalueindices::BitArray{1},
        ΔXt::R,
        ψXt::S,
        σXt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{Float64},
            T<:AbstractArray{Float64}}
    filteredΔXt = ΔXt[samplevalueindices]
    filterψXt = ψXt[samplevalueindices]
    filteredσXt = σXt[samplevalueindices]
    return sum(filteredΔXt .* filterψXt ./ (filteredσXt
        .* filteredσXt))
end
#Test written.

"""
    calculategirsanovvector(
        lengthvector::Int,
        samplevalueindices::Vector{BitArray{1}},
        samplevalues::S,
        ψXt::Vector{Array{Float64,1}},
        σXt::T) where {S<:AbstractArray{Float64}, T<:Union{Float64, U} where
            U<:AbstractArray{Float64}}
    calculategirsanovvector(
        lengthvector::Int,
        samplevalues::R,
        ψXt::S,
        σXt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{U} where
            U<:AbstractArray{Float64}, T<:Union{Float64, V} where
            V<:AbstractArray{Float64}}

Internal function, not exported!

Calculates the Girsanov vector.
"""
function calculategirsanovvector(
        lengthvector::Int,
        samplevalues::R,
        ψXt::S,
        σXt::T) where {R<:AbstractArray{Float64}, S<:AbstractArray{U} where
            U<:AbstractArray{Float64}, T<:Union{Float64, V} where
            V<:AbstractArray{Float64}}
    ΔXt = samplevalues[2:end] - samplevalues[1:end-1]
    return [calculategirsanovvectorelement(ΔXt, ψXt[k], σXt) for k in
        1:lengthvector]
end
# Test written.

function calculategirsanovvector(
        lengthvector::Int,
        samplevalueindices::Vector{BitArray{1}},
        samplevalues::S,
        ψXt::Vector{Array{Float64,1}},
        σXt::T) where {S<:AbstractArray{Float64}, T<:Union{Float64, U} where
            U<:AbstractArray{Float64}}
    ΔXt = samplevalues[2:end] - samplevalues[1:end-1]
    return [calculategirsanovvectorelement(samplevalueindices[k], ΔXt, ψXt[k],
        σXt) for k in 1:lengthvector]
end
# Test written.

"""
    calculateσXt(σ::Float64, v::AbstractArray{Float64}) = σ
    calculateσXt(σ::Function, v::AbstractArray{Float64}) = σ.(v)

Internal function, not exported!
"""
calculateσXt(σ::Float64, v::AbstractArray{Float64}) = σ
calculateσXt(σ::Function, v::AbstractArray{Float64}) = σ.(v)
#Tests written.

"""
    calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients,
        X::SamplePath, model::SDEModel)
    calculateposterior(Π::AbstractGaussianProcess,
        X::SamplePath, model::SDEModel)::GaussianProcess

Calculates the posterior distribution Π(⋅∣X), and returns it as a
GaussianProcess object. Uses model to determine the right likelihood.

# Examples

##Example Faber-Schauder expansion
```julia
model = SDEModel(1.0, 0.0, 10000.0, 0.1)
sde = SDE(x->sinpi(2*x), model)
X = rand(sde)
Π = FaberSchauderExpansionWithGaussianCoefficients([2.0^(-j) for j in 0:5])
postΠ = calculateposterior(Π, X, model)
```
##Example Fourier expansion
```julia
model = SDEModel(1.0, 0.0, 10000.0, 0.1)
sde = SDE(x->sinpi(2*x), model)
X = rand(sde)
distribution = GaussianVector(diagm([k^(-1.0) for k in 1.0:50.0]))
Π = GaussianProcess([fourier(k) for k in 1:50], distribution)
postΠ = calculateposterior(Π, X, model)
```
"""
function calculateposterior(Π::AbstractGaussianProcess,
        X::SamplePath, model::SDEModel)::GaussianProcess
    σXt = calculateσXt(model.σ, X.samplevalues[1:end-1])
    ψXt = [f.(X.samplevalues[1:end-1]) for f in Π.basis]
    lengthΠ = length(Π)
    girsanovvector = calculategirsanovvector(lengthΠ, X.samplevalues, ψXt, σXt)
    girsanovmatrix = calculategirsanovmatrix(lengthΠ, X.timeinterval, ψXt, σXt)
    precisionmatrixposterior = full(girsanovmatrix) + inv(full(Π.distribution.var))
    potentialposterior = girsanovvector + full(Π.distribution.mean)
    meanposterior = precisionmatrixposterior \ potentialposterior
    posteriordistribution = GaussianVector(meanposterior,
        precisionmatrixposterior^(-0.5))
    return GaussianProcess(Π.basis, posteriordistribution)
end

function calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients,
        X::SamplePath, model::SDEModel)
    samplevaluesmod1 = mod.(X.samplevalues[1:end-1], 1.0)
    σXt = calculateσXt(model.σ, X.samplevalues[1:end-1])
    ψXt = [f.(X.samplevalues[1:end-1]) for f in Π.basis]
    samplevalueindices = [Π.leftboundssupport[i] .≤
        samplevaluesmod1 .≤ Π.rightboundssupport[i] for i in 1:length(Π)]
    girsanovvector = calculategirsanovvector(length(Π), samplevalueindices,
        X.samplevalues, ψXt, σXt)
    girsanovmatrix = calculategirsanovmatrix(Π, samplevalueindices,
        X.timeinterval, ψXt, σXt)
    precisionmatrixposterior = full(girsanovmatrix) + inv(full(Π.distribution.var))
    potentialposterior = girsanovvector + mean(Π.distribution)
    meanposterior = precisionmatrixposterior \ potentialposterior
    posteriordistribution = GaussianVector(meanposterior,
        precisionmatrixposterior^(-0.5))
    return FaberSchauderExpansionWithGaussianCoefficients(Π.higestlevel, posteriordistribution)
end


text:
function mdl_line_ode!(device_states,
                        output_ode,
                        V_r_from,
                        V_i_from,
                        V_r_to,
                        V_i_to,
                        current_r_from,
                        current_i_from,
                        current_r_to,
                        current_i_to,
                        dv_from,
                        dv_to,
                        sys_f,
                        branch)

    L = branch.x
    R = branch.r
    ω_b = sys_f*2*π
    c_from = branch.b.from
    c_to = branch.b.to

    current_r_from[1] -= (1.0/ω_b)*c_from*dv_from[1] - c_from*V_i_from
    current_i_from[1] -= (1.0/ω_b)*c_from*dv_from[2] + c_from*V_r_from
    current_r_to[1] -= (1.0/ω_b)*c_to*dv_to[1] - c_to*V_i_to
    current_i_to[1] -= (1.0/ω_b)*c_to*dv_to[2] + c_to*V_r_to

    Il_r = device_states[1]
    Il_i = device_states[2]
    output_ode[1] = (ω_b/L)*((V_r_from - V_r_to) - (R*Il_r - L*Il_i))
    output_ode[2] = (ω_b/L)*((V_i_from - V_i_to) - (R*Il_i + L*Il_r))

    current_r_from[1] -= Il_r
    current_i_from[1] -= Il_i
    current_r_to[1] += Il_r
    current_i_to[1] += Il_i

end

ext: jl


text:
# Recursive UNION ALL node.

mutable struct IterateNode <: TabularNode
    over::Union{SQLNode, Nothing}
    iterator::SQLNode

    IterateNode(; over = nothing, iterator) =
        new(over, iterator)
end

IterateNode(iterator; over = nothing) =
    IterateNode(over = over, iterator = iterator)

"""
    Iterate(; over = nothing, iterator)
    Iterate(iterator; over = nothing)

`Iterate` generates the concatenated output of an iterated query.

The `over` query is evaluated first.  Then the `iterator` query is repeatedly
applied: to the output of `over`, then to the output of its previous run, and
so on, until the iterator produces no data.  All these outputs are concatenated
to generate the output of `Iterate`.

The `iterator` query should have an alias specified with [`As`](@ref); it can
refer to the output of the previous iteration using [`From`](@ref) with the same
alias.

```sql
WITH RECURSIVE iterator AS (
  SELECT ...
  FROM \$over
  UNION ALL
  SELECT ...
  FROM \$iterator
)
SELECT ...
FROM iterator
```

# Examples

```jldoctest
julia> q = Define(:n => 1, :f => 1) |>
           Iterate(
               From(:factorial) |>
               Define(:n => Get.n .+ 1) |>
               Define(:f => Get.f .* Get.n) |>
               Where(Get.n .<= 10) |>
               As(:factorial));

julia> print(render(q))
WITH RECURSIVE "factorial_1" ("n", "f") AS (
  SELECT
    1 AS "n",
    1 AS "f"
  UNION ALL
  SELECT
    ("factorial_1"."n" + 1) AS "n",
    ("factorial_1"."f" * ("factorial_1"."n" + 1)) AS "f"
  FROM "factorial_1"
  WHERE (("factorial_1"."n" + 1) <= 10)
)
SELECT
  "factorial_1"."n",
  "factorial_1"."f"
FROM "factorial_1"
```
"""
Iterate(args...; kws...) =
    IterateNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Iterate), pats::Vector{Any}) =
    dissect(scr, IterateNode, pats)

function PrettyPrinting.quoteof(n::IterateNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Iterate))
    if !ctx.limit
        push!(ex.args, quoteof(n.iterator, ctx))
    else
        push!(ex.args, :…)
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::IterateNode) =
    label(n.over)

rebase(n::IterateNode, n′) =
    IterateNode(over = rebase(n.over, n′), iterator = n.iterator)


text:
struct Quaternion{T<:Real} <: Number
    s::T
    v1::T
    v2::T
    v3::T
    norm::Bool
end

(::Type{Quaternion{T}})(x::Real) where {T} = Quaternion{T}(x, 0, 0, 0, false)
(::Type{Quaternion{T}})(q::Quaternion{T}) where {T<:Real} = q
(::Type{Quaternion{T}})(q::Quaternion) where {T<:Real} = Quaternion{T}(q.s, q.v1, q.v2, q.v3, q.norm)
Quaternion(s::Real, v1::Real, v2::Real, v3::Real, n::Bool = false) =
    Quaternion(promote(s, v1, v2, v3)..., n)
Quaternion(x::Real) = Quaternion(x, zero(x), zero(x), zero(x), abs(x) == one(x))
Quaternion(z::Complex) = Quaternion(z.re, z.im, zero(z.re), zero(z.re), abs(z) == one(z.re))
Quaternion(s::Real, a::Vector) = Quaternion(s, a[1], a[2], a[3])
Quaternion(a::Vector) = Quaternion(0, a[1], a[2], a[3])

convert(::Type{Quaternion{T}}, x::Real) where {T} =
    Quaternion(convert(T, x), convert(T, 0), convert(T, 0), convert(T, 0))
convert(::Type{Quaternion{T}}, z::Complex) where {T} =
    Quaternion(convert(T, real(z)), convert(T, imag(z)), convert(T, 0), convert(T, 0))
convert(::Type{Quaternion{T}}, q::Quaternion{T}) where {T <: Real} = q
convert(::Type{Quaternion{T}}, q::Quaternion) where {T} =
    Quaternion(convert(T, q.s), convert(T, q.v1), convert(T, q.v2), convert(T, q.v3), q.norm)

promote_rule(::Type{Quaternion{T}}, ::Type{T}) where {T <: Real} = Quaternion{T}
promote_rule(::Type{Quaternion}, ::Type{T}) where {T <: Real} = Quaternion
promote_rule(::Type{Quaternion{T}}, ::Type{S}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}
promote_rule(::Type{Complex{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}
promote_rule(::Type{Quaternion{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}

quat(p, v1, v2, v3, n = false) = Quaternion(p, v1, v2, v3, n)
quat(x) = Quaternion(x)
quat(s, a) = Quaternion(s, a)

function show(io::IO, q::Quaternion)
    pm(x) = x < 0 ? " - $(-x)" : " + $x"
    print(io, q.s, pm(q.v1), "im", pm(q.v2), "jm", pm(q.v3), "km")
end

real(::Type{Quaternion{T}}) where {T} = T
real(q::Quaternion) = q.s
imag(q::Quaternion) = [q.v1, q.v2, q.v3]

(/)(q::Quaternion, x::Real) = Quaternion(q.s / x, q.v1 / x, q.v2 / x, q.v3 / x)

conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3, q.norm)
abs(q::Quaternion) = sqrt(q.s * q.s + q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3)
float(q::Quaternion{T}) where T = convert(Quaternion{float(T)}, q)
abs_imag(q::Quaternion) = sqrt(q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3)
abs2(q::Quaternion) = q.s * q.s + q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3
inv(q::Quaternion) = q.norm ? conj(q) : conj(q) / abs2(q)

isfinite(q::Quaternion) = q.norm ? true : (isfinite(q.s) && isfinite(q.v1) && isfinite(q.v2) && isfinite(q.v3))

function normalize(q::Quaternion)
    if (q.norm)
        return q
    end
    q = q / abs(q)
    Quaternion(q.s, q.v1, q.v2, q.v3, true)
end

function normalizea(q::Quaternion)
    if (q.norm)
        return (q, one(q.s))
    end
    a = abs(q)
    q = q / a
    (Quaternion(q.s, q.v1, q.v2, q.v3, true), a)
end

function normalizeq(q::Quaternion)
    a = abs(q)
    if a > 0
        q = q / a
        Quaternion(q.s, q.v1, q.v2, q.v3, true)
    else
        Quaternion(0.0, 1.0, 0.0, 0.0, true)
    end
end

(-)(q::Quaternion) = Quaternion(-q.s, -q.v1, -q.v2, -q.v3, q.norm)

(+)(q::Quaternion, w::Quaternion) =
    Quaternion(q.s + w.s, q.v1 + w.v1, q.v2 + w.v2, q.v3 + w.v3)

(-)(q::Quaternion, w::Quaternion) =
    Quaternion(q.s - w.s, q.v1 - w.v1, q.v2 - w.v2, q.v3 - w.v3)

(*)(q::Quaternion, w::Quaternion) = Quaternion(q.s * w.s - q.v1 * w.v1 - q.v2 * w.v2 - q.v3 * w.v3,
                                               q.s * w.v1 + q.v1 * w.s + q.v2 * w.v3 - q.v3 * w.v2,
                                               q.s * w.v2 - q.v1 * w.v3 + q.v2 * w.s + q.v3 * w.v1,
                                               q.s * w.v3 + q.v1 * w.v2 - q.v2 * w.v1 + q.v3 * w.s,
                                               q.norm && w.norm)
(/)(q::Quaternion, w::Quaternion) = q * inv(w)

angleaxis(q::Quaternion) = angle(q), axis(q)

angle(q::Quaternion) = 2 * atan(√(q.v1^2 + q.v2^2 + q.v3^2), q.s)

function axis(q::Quaternion)
    q = normalize(q)
    s = sin(angle(q) / 2)
    abs(s) > 0 ?
        [q.v1, q.v2, q.v3] / s :
        [1.0, 0.0, 0.0]
end

argq(q::Quaternion) = normalizeq(Quaternion(0, q.v1, q.v2, q.v3))

function exp(q::Quaternion)
    s = q.s
    se = exp(s)
    scale = se
    th = abs_imag(q)
    if th > 0
        scale *= sin(th) / th
    end
    Quaternion(se * cos(th), scale * q.v1, scale * q.v2, scale * q.v3, abs(s) < eps(typeof(s)))
end

function log(q::Quaternion)
    q, a = normalizea(q)
    s = q.s
    M = abs_imag(q)
    th = atan(M, s)
    if M > 0
        M = th / M
        return Quaternion(log(a), q.v1 * M, q.v2 * M, q.v3 * M)
    else
        return Quaternion(log(a), th, 0.0, 0.0)
    end
end

function sin(q::Quaternion)
    L = argq(q)
    return (exp(L * q) - exp(-L * q)) / (2 * L)
end

function cos(q::Quaternion)
    L = argq(q)
    return (exp(L * q) + exp(-L * q)) / 2
end

(^)(q::Quaternion, w::Quaternion) = exp(w * log(q))

sqrt(q::Quaternion) = exp(0.5 * log(q))

function linpol(p::Quaternion, q::Quaternion, t::Real)
    p = normalize(p)
    q = normalize(q)
    qm = -q
    if abs(p - q) > abs(p - qm)
        q = qm
    end
    c = p.s * q.s + p.v1 * q.v1 + p.v2 * q.v2 + p.v3 * q.v3
    if c > - 1.0
        if c < 1.0
            o = acos(c)
            s = sin(o)
            sp = sin((1 - t) * o) / s
            sq = sin(t * o) / s
        else
            sp = 1 - t
            sq = t
        end
        Quaternion(sp * p.s  + sq * q.s,
                   sp * p.v1 + sq * q.v1,
                   sp * p.v2 + sq * q.v2,
                   sp * p.v3 + sq * q.v3, true)
    else
        s  =  p.v3
        v1 = -p.v2
        v2 =  p.v1
        v3 = -p.s
        sp = sin((0.5 - t) * pi)
        sq = sin(t * pi)
        Quaternion(s,
                   sp * p.v1 + sq * v1,
                   sp * p.v2 + sq * v2,
                   sp * p.v3 + sq * v3, true)
    end
end

quatrand()  = quat(randn(), randn(), randn(), randn())
nquatrand() = normalize(quatrand())

rand(r::AbstractRNG, ::SamplerType{Quaternion{T}}) where {T<:Real} = quat(rand(r, T), rand(r, T),rand(r,T),rand(r,T))

## Rotations

function qrotation(axis::Vector{T}, theta) where {T <: Real}
    if length(axis) != 3
        error("Must be a 3-vector")
    end
    u = normalize(axis)
    thetaT = convert(eltype(u), theta)
    s = sin(thetaT / 2)
    Quaternion(cos(thetaT / 2), s * u[1], s * u[2], s * u[3], true)
end

# Variant of the above where norm(rotvec) encodes theta
function qrotation(rotvec::Vector{T}) where {T <: Real}
    if length(rotvec) != 3
        error("Must be a 3-vector")
    end
    theta = norm(rotvec)
    if theta > 0
        s = sin(theta / 2) / theta  # divide by theta to make rotvec a unit vector
        return Quaternion(cos(theta / 2), s * rotvec[1], s * rotvec[2], s * rotvec[3], true)
    end
    Quaternion(one(T), zero(T), zero(T), zero(T), true)
end

function qrotation(dcm::Matrix{T}) where {T<:Real}
    # See https://arxiv.org/pdf/math/0701759.pdf
    a2 = 1 + dcm[1,1] + dcm[2,2] + dcm[3,3]
    b2 = 1 + dcm[1,1] - dcm[2,2] - dcm[3,3]
    c2 = 1 - dcm[1,1] + dcm[2,2] - dcm[3,3]
    d2 = 1 - dcm[1,1] - dcm[2,2] + dcm[3,3]

    if a2 >= max(b2, c2, d2)
        a = sqrt(a2)/2
        return Quaternion(a, (dcm[3,2]-dcm[2,3])/4a, (dcm[1,3]-dcm[3,1])/4a, (dcm[2,1]-dcm[1,2])/4a)
    elseif b2 >= max(c2, d2)
        b = sqrt(b2)/2
        return Quaternion((dcm[3,2]-dcm[2,3])/4b, b, (dcm[2,1]+dcm[1,2])/4b, (dcm[1,3]+dcm[3,1])/4b)
    elseif c2 >= d2
        c = sqrt(c2)/2
        return Quaternion((dcm[1,3]-dcm[3,1])/4c, (dcm[2,1]+dcm[1,2])/4c, c, (dcm[3,2]+dcm[2,3])/4c)
    else
        d = sqrt(d2)/2
        return Quaternion((dcm[2,1]-dcm[1,2])/4d, (dcm[1,3]+dcm[3,1])/4d, (dcm[3,2]+dcm[2,3])/4d, d)
    end
end

function qrotation(dcm::Matrix{T}, qa::Quaternion) where {T<:Real}
    q = qrotation(dcm)
    abs(q-qa) < abs(q+qa) ? q : -q
end

rotationmatrix(q::Quaternion) = rotationmatrix_normalized(normalize(q))

function rotationmatrix_normalized(q::Quaternion)
    sx, sy, sz = 2q.s * q.v1, 2q.s * q.v2, 2q.s * q.v3
    xx, xy, xz = 2q.v1^2, 2q.v1 * q.v2, 2q.v1 * q.v3
    yy, yz, zz = 2q.v2^2, 2q.v2 * q.v3, 2q.v3^2
    [1 - (yy + zz)     xy - sz     xz + sy;
        xy + sz   1 - (xx + zz)    yz - sx;
        xz - sy      yz + sx  1 - (xx + yy)]
end

function normalize(v::Vector{T}) where {T}
    nv = norm(v)
    if nv > 0
        return v / nv
    end
    zeros(promote_type(T, typeof(nv)), length(v))
end

function slerp(qa::Quaternion{T}, qb::Quaternion{T}, t::T) where {T}
    # http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/slerp/
    coshalftheta = qa.s * qb.s + qa.v1 * qb.v1 + qa.v2 * qb.v2 + qa.v3 * qb.v3;

    if coshalftheta < 0
        qm = -qb
        coshalftheta = -coshalftheta
    else
        qm = qb
    end
    abs(coshalftheta) >= 1.0 && return qa

    halftheta    = acos(coshalftheta)
    sinhalftheta = sqrt(one(T) - coshalftheta * coshalftheta)

    if abs(sinhalftheta) < 0.001
        return Quaternion(
            T(0.5) * (qa.s  + qb.s),
            T(0.5) * (qa.v1 + qb.v1),
            T(0.5) * (qa.v2 + qb.v2),
            T(0.5) * (qa.v3 + qb.v3),
        )
    end

    ratio_a = sin((one(T) - t) * halftheta) / sinhalftheta
    ratio_b = sin(t * halftheta) / sinhalftheta

    Quaternion(
        qa.s  * ratio_a + qm.s  * ratio_b,
        qa.v1 * ratio_a + qm.v1 * ratio_b,
        qa.v2 * ratio_a + qm.v2 * ratio_b,
        qa.v3 * ratio_a + qm.v3 * ratio_b,
    )
end


text:
"""
    MeshQuadrilateralModule

Module  for generation of meshes composed of quadrilaterals.
"""
module MeshQuadrilateralModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..FESetModule: AbstractFESet, FESetQ4, FESetQ8, bfun, cat, connasarray
import ..FENodeSetModule: FENodeSet, count
import ..MeshModificationModule: mergemeshes
import ..MeshUtilModule: makecontainer, addhyperface!, findhyperface!, linearspace, linearspace
import LinearAlgebra: norm
import Statistics: mean

"""
    Q4annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt)

Mesh of an annulus segment.

Mesh of an annulus segment, centered at the origin, with internal radius `rin`,
and  external radius `rex`, and  development angle `Angl` (in radians). Divided
into elements: nr, nc in the radial and circumferential direction respectively.
"""
function Q4annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt)
    trin=min(rin,rex);
    trex=max(rin,rex);
    fens,fes =Q4block(trex-trin,Angl,nr,nc);
    xy=fens.xyz;
    for i=1:count(fens)
        r=trin+xy[i,1]; a=xy[i,2];
        xy[i,:]=[r*cos(a) r*sin(a)];
    end
    fens.xyz=xy;
    return fens,fes
end

"""
    Q4quadrilateral(xyz::FFltMat, nL::FInt, nW::FInt)

Mesh of a general quadrilateral given by the location of the vertices.
"""
function Q4quadrilateral(xyz::FFltMat, nL::FInt, nW::FInt)
    npts=size(xyz,1);
    if npts==2 # In this case the quadrilateral must be defined in two dimensions
        lo=minimum(xyz, dims = 1);
        hi=maximum(xyz, dims = 1);
        xyz=[[lo[1] lo[2]];
            [hi[1] lo[2]];
            [hi[1] hi[2]];
            [lo[1] hi[2]]];
    elseif npts!=4
        error("Need 2 or 4 points");
    end

    fens,fes = Q4block(2.,2.,nL,nW);

    xyz1=fens.xyz;
    if (size(xyz1,2)<size(xyz,2))
        nxyz1=zeros(FFlt,size(xyz1,1),size(xyz,2));
        nxyz1[:,1:size(xyz1,2)]=xyz1;
        xyz1=nxyz1;
    end

    dummy = FESetQ4(reshape(collect(1:4),1,4))
    pxyz = xyz1;
    for i = 1:count(fens)
        N = bfun(dummy, broadcast(-, pxyz[i,:], 1.0));# shift coordinates by -1
        pxyz[i,:] = N'*xyz;
    end
    fens.xyz = deepcopy(pxyz);
    return fens,fes
end

"""
    Q4elliphole(xradius::FFlt, yradius::FFlt, L::FFlt, H::FFlt,
      nL::FInt, nH::FInt, nW::FInt)

Mesh of one quarter of a rectangular plate with an elliptical hole.

    xradius,yradius = radius of the ellipse,
    L,H= and dimensions of the plate,
    nL,nH= numbers of edges along the side of the plate; this also happens
      to be the number of edges along the circumference of the elliptical
      hole
    nW= number of edges along the remaining straight edge (from the hole
      in the direction of the length),
"""
function Q4elliphole(xradius::FFlt, yradius::FFlt, L::FFlt, H::FFlt,
    nL::FInt, nH::FInt, nW::FInt)
    dA =pi/2/(nL +nH);
    tolerance = (xradius+yradius)/(nL*nH)/100;
    fens= nothing; fes= nothing;
    for i= 1:nH
        xy = [xradius*cos((i-1)*dA) yradius*sin((i-1)*dA);
        L (i-1)/nH*H;
        L (i)/nH*H;
        xradius*cos((i)*dA) yradius*sin((i)*dA)];
        fens1,fes1 = Q4quadrilateral(xy,nW,1);
        if (fens == nothing)
            fens = fens1; fes = fes1;
        else
            fens,fes1,fes2 = mergemeshes(fens1, fes1, fens, fes, tolerance);
            fes = cat(fes1,fes2);
        end
    end
    for i= 1:nL
        xy = [xradius*cos((nH+i-1)*dA)   yradius*sin((nH+i-1)*dA);
        (nL-i+1)/nL*L   H;
        (nL-i)/nL*L  H;
        xradius*cos((nH+i)*dA)   yradius*sin((nH+i)*dA)];
        fens1,fes1 = Q4quadrilateral(xy,nW,1);
        fens,fes1,fes2 = mergemeshes(fens1, fes1, fens, fes, tolerance);
        fes = cat(fes1, fes2);
    end
    return fens,fes
end

"""
    Q4block(Length::FFlt, Width::FFlt, nL::FInt, nW::FInt)

Mesh of a rectangle, Q4 elements.

Divided into elements: nL, nW in the first, second (x,y).
"""
function Q4block(Length::FFlt, Width::FFlt, nL::FInt, nW::FInt)
    return Q4blockx(collect(linearspace(0.0,Length,nL+1)), collect(linearspace(0.0,Width,nW+1)));
end

"""
    Q4blockx(xs::FFltVec, ys::FFltVec)

Graded mesh  of a rectangle, Q4 finite elements.

    Mesh of a 2-D block, Q4 finite elements. The nodes are located at the
    Cartesian product of the two intervals on the input.  This allows for
    construction of graded meshes.

    xs,ys - Locations of the individual planes of nodes.
"""
function Q4blockx(xs::FFltVec, ys::FFltVec)
    nL = length(xs) - 1;
    nW = length(ys) - 1;

    nnodes = (nL+1) * (nW+1);
    ncells = nL * nW;

    # preallocate node locations
    xyz = zeros(FFlt, nnodes, 2);
    k = 1;
    for j = 1:(nW+1)
        for i = 1:(nL+1)
            xyz[k, 1] = xs[i]
            xyz[k, 2] = ys[j]
            k = k + 1;
        end
    end
    # create the nodes
    fens = FENodeSet(xyz);

    #preallocate connectivity matrix
    conn = zeros(FInt, ncells, 4);

    # function  nodenumbers(i,j,nL,nW)
    #     f = (j-1) * (nL+1) + i;
    #     nn = [f, (f+1), f+(nL+1)+1, f+(nL+1)];
    #     return nn
    # end

    k = 1;
    for i = 1:nL
        for j = 1:nW
            f = (j-1) * (nL+1) + i;
            conn[k, 1] = f
            conn[k, 2] = (f+1)
            conn[k, 3] = f+(nL+1)+1
            conn[k, 4] = f+(nL+1)
            k = k + 1;
        end
    end
    # create the cells
    fes = FESetQ4(conn);

    return fens,fes;
end

function Q4blockx(xs::AbstractVector, ys::AbstractVector)
    return Q4blockx(FFltVec(xs), FFltVec(ys))
end

"""
    Q8block(Length::FFlt, Width::FFlt, nL::FInt, nW::FInt)

Mesh of a rectangle of Q8 elements.
"""
function Q8block(Length::FFlt, Width::FFlt, nL::FInt, nW::FInt)
    fens,fes  = Q4block(Length, Width, nL, nW);
    fens,fes = Q4toQ8(fens, fes);
end

"""
    Q4toQ8(fens::FENodeSet, fes::FESetQ4)

Convert a mesh of quadrilateral Q4 to quadrilateral Q8.
"""
function Q4toQ8(fens::FENodeSet, fes::FESetQ4)
    nedges=4;
    ec = [1  2; 2  3; 3  4; 4  1];
    conns = connasarray(fes);
    # Additional node numbers are numbered from here
    newn = count(fens)+1;
    # make a search structure for edges
    edges = makecontainer();
    for i= 1:size(conns,1)
        conn = conns[i,:];
        for J = 1:nedges
            ev = conn[ec[J,:]];
            newn = addhyperface!(edges, ev, newn);
        end
    end
    xyz1 =fens.xyz;             # Pre-existing nodes
    # Allocate for vertex nodes plus edge nodes plus face nodes
    xyz =zeros(FFlt,newn-1,size(xyz1,2));
    xyz[1:size(xyz1,1),:] = xyz1; # existing nodes are copied over
    # calculate the locations of the new nodes
    # and construct the new nodes
    for i in keys(edges)
        C=edges[i];
        for J = 1:length(C)
            ix = vec([item for item in C[J].o])
            push!(ix,  i) # Add the anchor point as well
            xyz[C[J].n, :] = mean(xyz[ix, :], dims = 1);
        end
    end
    # construct new geometry cells
    nconns =zeros(FInt,size(conns,1),8);
    nc=1;
    for i= 1:size(conns,1)
        conn = conns[i,:];
        econn=zeros(FInt,1,nedges);
        for J = 1:nedges
            ev=conn[ec[J,:]];
            h,n = findhyperface!(edges, ev);
            econn[J]=n;
        end
        nconns[nc,:] =vcat(vec(conn), vec(econn));
        nc= nc+ 1;
    end
    fens = FENodeSet(xyz);
    fes = FESetQ8(nconns);
    return fens, fes
end

"""
    Q8blockx(xs::FFltVec, ys::FFltVec)

Graded mesh of a 2-D block of Q8 finite elements.
"""
function Q8blockx(xs::FFltVec, ys::FFltVec)
    fens, fes = Q4blockx(xs, ys);
    fens, fes = Q4toQ8(fens, fes);
end

"""
    Q4refine(fens::FENodeSet, fes::FESetQ4)

Refine a mesh of quadrilaterals by bisection.
"""
function Q4refine(fens::FENodeSet, fes::FESetQ4)
    nedges=4;
    ec = [1  2; 2  3; 3  4; 4  1];
    # make a search structure for edges
    # Additional node numbers are numbered from here
    newn=count(fens)+1;
    # make a search structure for edges
    edges=makecontainer();
    for i= 1:length(fes.conn)
        for J = 1:nedges
            ev=fes.conn[i][ec[J,:]];
            newn = addhyperface!(edges, ev, newn);
        end
    end
    newn=  newn+length(fes.conn) # add the interior nodes to the total
    xyz1 =fens.xyz;             # Pre-existing nodes
    # Allocate for vertex nodes plus edge nodes plus face nodes
    xyz =zeros(FFlt,newn-1,size(xyz1,2));
    xyz[1:size(xyz1,1),:] = xyz1; # existing nodes are copied over
    # calculate the locations of the new nodes
    # and construct the new nodes
    for i in keys(edges)
        C=edges[i];
        for J = 1:length(C)
            ix = vec([item for item in C[J].o])
            push!(ix, i)
            xyz[C[J].n,:] = mean(xyz[ix,:], dims = 1);
        end
    end
    # construct new geometry cells: for new elements out of one old one
    nconn =zeros(FInt,4*length(fes.conn),4);
    nc=1;
    for i= 1:length(fes.conn)
        econn=zeros(FInt,1,nedges);
        for J = 1:nedges
            ev=fes.conn[i][ec[J,:]];
            h,n=findhyperface!(edges, ev);
            econn[J]=n;
        end

        inn=size(xyz,1)-length(fes.conn)+i

        xyz[inn,:]=mean(xyz[[k for k in fes.conn[i]],:], dims = 1); # interior node
        #h,inn=findhyperface!(faces, conn);
        nconn[nc,:] =[fes.conn[i][1] econn[1] inn econn[4]];
        nc= nc+ 1;
        nconn[nc,:] =[fes.conn[i][2] econn[2] inn econn[1]];
        nc= nc+ 1;
        nconn[nc,:] =[fes.conn[i][3] econn[3] inn econn[2]];
        nc= nc+ 1;
        nconn[nc,:] =[fes.conn[i][4] econn[4] inn econn[3]];
        nc= nc+ 1;
    end
    fens =FENodeSet(xyz);
    nfes = FESetQ4(nconn);
    return fens,nfes            # I think I should not be overwriting the input!
end

"""
    Q8annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt)

Mesh of an annulus segment.

Mesh of an annulus segment, centered at the origin, with internal radius
rin, and  external radius rex, and  development angle Angl. Divided into
elements: nr, nc in the radial and circumferential direction
respectively.
"""
function Q8annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt)
    trin=min(rin,rex);
    trex=max(rin,rex);
    fens,fes = Q8block(trex-trin,Angl,nr,nc);
    xy=fens.xyz;
    for i=1:count(fens)
        r=trin+xy[i,1]; a=xy[i,2];
        xy[i,:]=[r*cos(a) r*sin(a)];
    end
    fens.xyz=xy;
    return fens,fes
end

function _ontosphere!(xyz, radius)
    for j = 1:size(xyz, 1)
        xyz[j,:] = xyz[j,:]*radius/norm(xyz[j,:]);
    end
end

# % Surface mesh of 1/8 of a sphere with a given number of elements per circumference.
# %
# % function [fens,fes]=Q4_sphere_n(radius,nperradius,thickness)
# %
# % Create a mesh of 1/8 of the sphere of "radius". The  mesh will consist of
# % 3*(nperradius/2)^2 quadrilateral elements, which corresponds to 4*nperradius
# % element edges per circumference. The parameter nperradius should be an even
# % number; if that isn't so is adjusted to by adding one.
# %
# % The quadrilaterals have thickness "thickness".
# %
# % Examples:
# % [fens,fes]=Q4_sphere_n(77.1,5,1.0);
# % drawmesh({fens,fes},'nodes','fes','facecolor','y', 'linewidth',2); hold on
# %
# % See also: Q4_sphere

"""
    Q4spheren(radius::FFlt, nperradius)

Generate mesh of a sphere.
"""
function Q4spheren(radius::FFlt, nperradius)
    if (mod(nperradius,2) != 0)
        nperradius = nperradius+1;
    end
    nL = Int(nperradius/2); nW = Int(nperradius/2);
    tolerance = radius / nperradius / 100;
    a = sqrt(2.0)/2;
    b = 1/sqrt(3.0);
    c = 0.6*a;
    d = 0.6*b;
    xyz = [1 0 0; 0 1 0; 0 0 1; a a 0; 0 a a; a 0 a; b b b];
    conn = [1 4 7 6; 4 2 5 7; 3 6 7 5];
    fens,fes = Q4quadrilateral(xyz[conn[1,:],:], nL, nW)
    fens1,fes1 = Q4quadrilateral(xyz[conn[2,:],:], nL, nW)
    fens,fes1,fes2 = mergemeshes(fens1, fes1, fens, fes, tolerance);
    fes = cat(fes1, fes2);
    fens1,fes1 = Q4quadrilateral(xyz[conn[3,:],:], nL, nW)
    fens,fes1,fes2 = mergemeshes(fens1, fes1, fens, fes, tolerance);
    fes = cat(fes1, fes2);
    _ontosphere!(fens.xyz, radius);
    return  fens,fes
end

"""
    Q4circlen(radius::FFlt, nperradius)

Mesh of a quarter circle with a given number of elements per radius.

The parameter `nperradius` should be an even
number; if that isn't so is adjusted to by adding one.
"""
function Q4circlen(radius::FFlt, nperradius)
    fens,fes = Q4spheren(radius, nperradius);
    # % apply transformation to project the locations of the nodes into the
    # % plane x-y
    for j=1:count(fens)
    	r = norm(fens.xyz[j,3])
        fens.xyz[j,1:2] = fens.xyz[j,1:2]*((radius-r)+r/2)/radius;
        fens.xyz[j, 3] = 0.0
    end
    return fens,fes
end

end


text:
using StructuredDistances
using Test

include("Distances/Multisets/matching_distances_test.jl")


text:
using JSON
using Pipe: @pipe
using MetaGen
using UnicodePlots

include("helper_function.jl")
include("scripts/useful_functions.jl")

#dict = @pipe "../../metagen-data/data_labelled/data_labelled.json" |> open |> read |> String |> JSON.parse
dict = @pipe "../../scratch_work_07_16_21/09_18/data_labelled.json" |> open |> read |> String |> JSON.parse

num_videos = 50
num_frames = 20
ground_truth_world_states = get_ground_truth(dict, num_videos)

function fit_threshold_NN(num_videos::Int64, num_frames::Int64, dict::Any,
    ground_truth_world_states::Vector{Any}, threshold::Float64)

    params = Video_Params(n_possible_objects = 5)
    receptive_fields = make_receptive_fields(params)
    objects_observed, camera_trajectories = make_observations_office(dict, receptive_fields, num_videos, num_frames, threshold)

    sim_NN = zeros(num_videos)
    for v = 1:num_videos
        for f = 1:num_frames
            #println("v ", v)
            #println("f ", f)
            gt_categories = threeD2twoD_new(ground_truth_world_states[v], v, f, params, camera_trajectories)
            NN_categories = last.(objects_observed[v,f])
            #println(gt_categories)
            #println(NN_categories)
            sim_NN[v] = sim_NN[v] + jaccard_similarity(gt_categories, NN_categories)
            #println(jaccard_similarity(gt_categories, NN_categories)) #printing after doesn't work since running jaccard_similarity actually changes the inputs / outputs
        end
        sim_NN[v] = sim_NN[v]/num_frames
    end
    return sum(sim_NN)/num_videos #averaging across videos
end

thresholds = collect(0:.01:1.0)
#thresholds = collect([0.1])
vs = fill(num_videos, length(thresholds))
fs = fill(num_frames, length(thresholds))
dicts = fill(dict, length(thresholds))
gts = fill(ground_truth_world_states, length(thresholds))
how_good = map(fit_threshold_NN, vs, fs, dicts, gts, thresholds)

println(thresholds[findmax(how_good)[2]])

plt = lineplot(thresholds, how_good)

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using OptionsMod

dispflags = 0
# dispflags = Optim.FINAL | Optim.ITER #| Optim.LINESEARCH | Optim.BRACKET | Optim.BISECT
ops = @options display=dispflags

gtol = 1e-5

# Quadratic objective function
function quadratic(g, x, A, b)
    calc_grad = !(g === nothing)
    v = x'*A*x/2 + b'*x
    if calc_grad
        A_mul_B(g, A, x)
        for i = 1:length(g)
            g[i] += b[i]
        end
    end
    return v[1]
end

# Matlab objective function in objfun.m
# Minimum value is 0, which is tricky
function mobjfun(g, x)
    calc_grad = !(g === nothing)
    x1 = x[1]
    x2 = x[2]
    ex1 = exp(x1)
    v = ex1*(4x1^2 + 2x2^2 + 4x1*x2 + 2x2 + 1)
    if calc_grad
        g[1] = v + ex1*(8x1 + 4x2)
        g[2] = ex1*(4x2 + 4x1 + 2)
    end
    return v
end

# driver1 problem in CG_DESCENT
function driver1{T}(g, x::Array{T})
    calc_grad = !(g === nothing)
    f = zero(T)
    n = length(x)
    for i = 1:n
        t = sqrt(convert(T, i))
        ex = exp(x[i])
        f += ex - t*x[i]
        if calc_grad
            g[i] = ex-t
        end
    end
    return f
end

# Quadratic
N = 8
A = randn(N,N)
A = A'*A
b = randn(N)
x0 = randn(N)
func = (g, x) -> quadratic(g, x, A, b)
x, fval, fcount, converged = cgdescent(func, x0, ops)
@assert converged
@assert all(abs(A*x + b) .< gtol)

l = fill(-2.0, N)
u = fill(2.0, N)
x0 = rand(N)-0.5
x, fval, fcount, converged = fminbox(func, x0, l, u, ops)
@assert converged

# Matlab objective function
x0 = Float64[-1,1]
g = Array(Float64, 2)
mobjfun(g, x0)
fun = x->mobjfun(nothing, x)
d = Calculus.finite_difference(fun, x0)
@assert all(abs(d-g) .< 1e-8)
func = (g, x) -> mobjfun(g, x)
x, fval, fcount, converged = cgdescent(func, x0, ops)
@assert converged

# cgdescent driver1
n = 100
x0 = ones(n)
func = (g,x) -> driver1(g, x)
x, fval, fcount, converged = cgdescent(func, x0, ops)
@assert converged

# Check nnls
A = randn(50,30)
b = randn(size(A, 1))
xnn, fval, fcount, converged = nnls(A, b)
@assert converged


text:
"Convert a structure into a Dict but only with the properties with non-nothing values."
function non_nothing_dict(s, d::Dict)
    for fname ∈ propertynames(s)
        val = getproperty(s, fname)
        if !(val === nothing)
            d[fname] = val
        end
    end
    d
end
non_nothing_dict(s) = non_nothing_dict(s, Dict())

function nothing_or_type(var, type)
    var === nothing || typeassert(var, type)
    true
end


text:
############################################################################################
"""
$(SIGNATURES)
Callable struct to make initializing MCMC sampler easier in sampling library.

# Examples
```julia
```

"""
struct MCMCConstructor{M,S<:Union{Symbol,NTuple{k,Symbol} where k},D<:MCMCDefault} <: AbstractConstructor
    "Valid MCMC kernel."
    kernel::M
    "Parmeter to be tagged in MCMC sampler."
    sym::S
    "MCMC Default Arguments"
    default::D
    function MCMCConstructor(
        kernel::Type{M}, sym::S, default::D
    ) where {M<:MCMCKernel,S<:Union{Symbol,NTuple{k,Symbol} where k},D<:MCMCDefault}
        tup = BaytesCore.to_Tuple(sym)
        return new{typeof(kernel),typeof(tup),D}(kernel, tup, default)
    end
end
function (constructor::MCMCConstructor)(
    _rng::Random.AbstractRNG,
    model::ModelWrapper,
    data::D,
    temperature::F,
    info::BaytesCore.SampleDefault
) where {D, F<:AbstractFloat}
    return MCMC(
        _rng,
        constructor.kernel,
        Objective(model, data, Tagged(model, constructor.sym), temperature),
        constructor.default,
        info
    )
end
function MCMC(
    kernel::Type{M}, sym::S; kwargs...
) where {M<:MCMCKernel,S<:Union{Symbol,NTuple{k,Symbol} where k}}
    return MCMCConstructor(kernel, sym, MCMCDefault(; kwargs...))
end

############################################################################################
function infer(diagnostics::Type{AbstractDiagnostics}, kernel::MCMCKernel)
    return println("No known diagnostics for given kernel")
end

"""
$(SIGNATURES)
Infer MCMC diagnostics type.

# Examples
```julia
```

"""
function infer(
    _rng::Random.AbstractRNG,
    diagnostics::Type{AbstractDiagnostics},
    mcmc::MCMC,
    model::ModelWrapper,
    data::D,
) where {D}
#    Tlik = eltype(mcmc.kernel.result.θᵤ)
    TKernel = infer(_rng, diagnostics, mcmc.kernel, model, data)
    TPrediction = infer(_rng, mcmc, model, data)
    TGenerated = infer_generated(_rng, mcmc, model, data)
    return MCMCDiagnostics{TPrediction,TKernel,TGenerated}
end

"""
$(SIGNATURES)
Infer type of predictions of MCMC sampler.

# Examples
```julia
```

"""
function infer(_rng::Random.AbstractRNG, mcmc::MCMC, model::ModelWrapper, data::D) where {D}
    objective = Objective(model, data, mcmc.tune.tagged)
    return typeof(predict(_rng, mcmc, objective))
end

"""
$(SIGNATURES)
Infer type of generated quantities of MCMC sampler.

# Examples
```julia
```

"""
function infer_generated(
    _rng::Random.AbstractRNG, mcmc::MCMC, model::ModelWrapper, data::D
) where {D}
    objective = Objective(model, data, mcmc.tune.tagged)
    return typeof(generate(_rng, objective, mcmc.tune.generated))
end

############################################################################################
"""
$(SIGNATURES)
Print result for a single trace.

# Examples
```julia
```

"""
function results(
    diagnosticsᵛ::AbstractVector{M}, mcmc::MCMC, Ndigits::Integer, quantiles::Vector{T}
) where {T<:Real,M<:MCMCDiagnostics}
    println(
        "### ",
        Base.nameof(typeof(mcmc.kernel)),
        " parameter target: ",
        keys(mcmc.tune.tagged.parameter),
    )
    println(
        "Sampler finished after ",
        size(diagnosticsᵛ, 1),
        " iterations with acceptance rates of ",
        round.(
            mean(diagnosticsᵛ[iter].accept.accepted for iter in eachindex(diagnosticsᵛ)) *
            100;
            digits=Ndigits,
        ),
        "%.",
    )
    println(
        "Avg. initial ℓposterior: ",
        round(mean(diagnosticsᵛ[begin].base.ℓobjective); digits=Ndigits),
        ", Avg. final ℓposterior: ",
        round(mean(diagnosticsᵛ[end].base.ℓobjective); digits=Ndigits),
        ".",
    )
    ## Print kernel specific diagnostics
    results(
        [diagnosticsᵛ[iter].kernel for iter in eachindex(diagnosticsᵛ)], Ndigits, quantiles
    )
    ## Print Divergences
    print_divergences(diagnosticsᵛ, mcmc.tune.phase)
    return nothing
end

############################################################################################
function result!(mcmc::MCMC, result::L) where {L<:ℓObjectiveResult}
    mcmc.kernel.result = result
    return nothing
end

function get_result(mcmc::MCMC)
    return mcmc.kernel.result
end

function predict(_rng::Random.AbstractRNG, mcmc::MCMC, objective::Objective)
    return predict(_rng, objective)
end

############################################################################################
# Export
export MCMCConstructor, infer, predict


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

# Floating point numbers - basic tests
@test 4.00000000000001 ≈ 4.0
@test 5.0 ≈ 4.999999999999993
@test 4.000000002 ≉ 4.00300002

# Other tolerance levels
@test isapprox(4.32, 4.3; rtol=0.1, atol=0.01)
@test isapprox(1.001, 1.002; rtol=0.001, atol=0.0001)
@test !isapprox(4.5, 4.9; rtol=0.001, atol=0.001)

# Complex numbers
@test 1.0 + 1.0im ≈ 1.0 + 1.00000000000001im
@test 0.9999999999999 + 1.0im ≈ 1.0 + 1.000000000000001im
@test isapprox(0.9999 + 1.0im, 1.0 + 1.1im; rtol = 0.0001, atol=1.1)

# Complex <-> reals
@test 1.0 + 0im ≈ 1.0000000000001
@test isapprox(0.9999999999999, 1.0 + 0im)
@test !isapprox(1.0+1im, 1.000000000000001)

# Comparing NaNs
@test 4.0 ≉ NaN
@test NaN ≉ 4.0
@test complex(2.3,NaN) ≉ complex(NaN,2.3)
@test NaN ≉ NaN
@test complex(NaN,NaN) ≉ complex(NaN,NaN)
@test complex(NaN,2.3) ≉ complex(NaN,2.3)
@test complex(2.3,NaN) ≉ complex(2.3,NaN)

# Comparing Infs
@test Inf ≈ Inf
@test Inf ≉ 1
@test Inf ≉ -Inf
@test complex(0.0,Inf) ≈ complex(0.0,Inf)
@test complex(0.0,Inf) ≉ complex(0.0,-Inf)

# Tests for integers and rationals
@test isapprox(4,4)
@test isapprox(4,5; atol=1)
@test !isapprox(4,6; atol=1)
@test !isapprox(4,5)

@test isapprox(1//2+3im, 1//2+3im)
@test isapprox(1//3, 0.33333333333333333)
@test isapprox(1//3, 0.3333; rtol=0.0001, atol=0.0001)
@test isapprox(1+1im, 1im+1)

# Notably missing from this test suite at the moment
# * Tests for other magnitudes of numbers - very small, very big, and combinations of small and big
# * Tests for various odd combinations of types, e.g. isapprox(x::Integer, y::Rational)

# issue #12375:
@test !isapprox(1e17, 1)

# Tests for arrays:
@test [1,2,3] ≈ [1,2,3+1e-9]
@test [0,1] ≈ [1e-9, 1]
@test [0,Inf] ≈ [0,Inf]
@test [0,Inf] ≉ [0,-Inf]

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using Makie, GeometryTypes, Distributions
import Makie: convert_arguments
import Bridge: mcsvd3, visualize_uncertainty

convert_arguments(P::Type{<:Union{Lines,Scatter}}, X::SamplePath{<:AbstractVector}) = convert_arguments(P, X.yy)
convert_arguments(P::Type{<:Union{Lines,Scatter}}, X::SamplePath{<:Real}) = convert_arguments(P, X.tt, X.yy)

function band!(scene, x, ylower, yupper; nargs...)
    n = length(x)
    coordinates = [x ylower; x yupper]
    ns = 1:n-1
    ns2 = n+1:2n-1
    connectivity = [ns ns .+ 1 ns2;
                    ns2 ns2 .+ 1 ns .+ 1]
    mesh!(scene, coordinates, connectivity; shading = false, nargs...)
end

"""
    mcsvd3(mcstates) -> mean, q, sv

Compute mean and decompose variance into singular values `sv` and
3d rotation quaternion `q`.

Returns vectors of appropriate `GeometryTypes`.
"""
function mcsvd3(states)
    Xmean = Point3f0[]
    Xscal = Point3f0[]
    Xrot = Vec4f0[]

    for i in 1:length(states)
        xx, vv = Bridge.mcstats(states[i])
        svds = svd.(vv)
        append!(Xmean, xx)
        append!(Xscal, map(svd->sqrt.(svd[2]), svds))
        append!(Xrot, Bridge.quaternion.(first.(svds)))
        if i != length(states)
            pop!(Xmean); pop!(Xscal); pop!(Xrot)
        end
    end
    Xmean, Xrot, Xscal
end

function visualize_uncertainty(scene, X, skip = 10, qu = 0.95; args...)
    c = sqrt(quantile(Chisq(3), qu))
    Xmean, Xrot, Xscal = X
    sphere = Sphere(Point3f0(0,0,0), 1.0f0)
    viz = visualize(
        (sphere, (Xmean[1:skip:end])),
        scale = c*(Xscal[1:skip:end]),
        rotation = (Xrot[1:skip:end]);
        args...
    ).children[]
    Makie.insert_scene!(scene, :mesh, viz, Dict(:show=>true, :camera=>:auto))
end

viridis(X, alpha = 0.9f0, maxviri = 200) = map(x->RGBA(Float32.(x)..., alpha), Bridge._viridis[round.(Int, range(1, stop=maxviri, length=length(X)))])
viridis(n::Integer, alpha = 0.9f0, maxviri = 200) = map(x->RGBA(Float32.(x)..., alpha), Bridge._viridis[round.(Int, range(1, stop=n>1 ? maxviri : 1, length=n))])

#set_perspective!(scene, perspective) = (push!(Makie.getscreen(scene).cameras[:perspective].view, perspective); scene)
#get_perspective(scene) = Makie.getscreen(scene).cameras[:perspective].view.value


text:
"""
nonparametric regression and quantile regression
* uses pre-generated kernel weights
* can be local constant, local linear, or local quadratic

y is n X p

x is n X dimx

xeval is neval X dimx

weights is n X neval (different weights for each eval. point)

weights should sum to one by columns, they are typical
nonparametric weights from a kernel.

keyword order=0,1 or 2 for local consant, local linear, or local quadratic
(default local linear)

keyword do_median (default false)
keyword do_ci to compute .05 and 0.95 quantiles

execute npreg() for an example
"""

using Random, Econometrics, Plots
function npreg(bandwidth=-1)
    println("npreg(), with no arguments, runs a simple example")
    println("npreg(bw) will run the example with your chosen bandwidth (Float64)")
    println("execute edit(npreg,()) to see the code")
    k = 1 # number of regressors
    Random.seed!(1) # set seed to
    n = 1000
    bandwidth == -1 ? bandwidth = 0.25*n^(-1.0/(4 + k)) : nothing
    neval = 100
    x = rand(n)*pi*2.0 # from 0 to 2π
    xeval = collect(range(pi/2., stop=pi*1.5, length=neval)) # from π/2 to 3π/2
    y = cos.(x) .+ 0.25 .* cos.(3.0.*x) + 0.1*randn(n)
    # npreg wants args to be matrices, not vectors, so convert them
    y = reshape(y,n,1)
    x = reshape(x,n,1)
    xeval = reshape(xeval,neval,1)
    ytrue = cos.(xeval) .+ 0.25.*cos.(3.0*xeval)
    weights = kernelweights(x, xeval, bandwidth, true, "knngaussian", 200)
    yhat, y50, y05, y95 = npreg(y, x, xeval, weights, order=1, do_median=true, do_ci=true)
    labels = ["true" "mean" "median" "0.05 quantile" "0.95 quantile"]
    title = "Kernel regression and quantiles"
    p = Plots.plot(xeval, [ytrue yhat y50 y05 y95], labels=labels, title=title)
    display(p)
    nothing
end

function npreg(y, x, xeval, weights; order=1, do_median=false, do_ci=false)
    weights = sqrt.(weights)
    neval, dimx = size(xeval)
    n, dimy = size(y)
    # local constant
    if order==0
        X = ones(n,1)
        Xeval = ones(neval,1)
    elseif order==1
    # local linear
        X = [ones(n,1) x]
        Xeval = [ones(neval,1) xeval]
    else
    # local quadratic
        stacked = [x; xeval]
        # cross products
        CP = zeros(n+neval, Int((dimx-1)*dimx/2))
        cpind = 0
        for i = 1:dimx-1
            for j = (i+1):dimx
                cpind += 1
                CP[:,cpind] = stacked[:,[i]].*stacked[:,[j]]
            end
        end
        ZZ = [ones(n+neval,1) stacked CP]
        X = view(ZZ,1:n,:)
        Xeval = view(ZZ,(n+1):n+neval,:)
    end
    # do the fit
    yhat = zeros(neval, dimy)
    y50 = nothing
    y05 = nothing
    y95 = nothing
    if do_median y50 = zeros(neval, dimy) end
    if do_ci
        y05 = zeros(neval, dimy)
        y95 = zeros(neval, dimy)
    end
    for i = 1:neval
        WX = weights[:,i] .* X
        Wy = weights[:,i] .* y
        b = WX\Wy
        yhat[i,:] = Xeval[i:i,:]*b
        if do_median
            for j = 1:dimy
                y50[i,j] = (Xeval[i:i,:]*qreg_coef(Wy[:,j], WX, 0.5))[1,1]
            end
        end
        if do_ci
            for j = 1:dimy
                y05[i,j] = (Xeval[i:i,:]*qreg_coef(Wy[:,j], WX, 0.05))[1,1]
                y95[i,j] = (Xeval[i:i,:]*qreg_coef(Wy[:,j], WX, 0.95))[1,1]
            end
        end
    end
    if do_median==false
        return yhat
    elseif do_ci==false
        return yhat, y50
    else
        return yhat, y50, y05, y95
    end
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using Test
using QuantumPropagators
using LinearAlgebra
using UnicodePlots

@testset "TLS Rabi Cycling" begin

    # We use Rabi cycling in a TLS as a test case, which allows to to compare
    # the propagation with the known analytic solution.

    SHOWPLOT = false

    Ψ0 = ComplexF64[1, 0]
    Ĥ = ComplexF64[0  0.5; 0.5 0]
    tlist = collect(range(0, 1.5π, length=101)) # 3π/2 pulse

    genfunc(tlist, i; kwargs...) = Ĥ

    storage = init_storage(Ψ0, tlist)
    @test isa(storage, Matrix)
    @test eltype(storage) == ComplexF64

    Ψ_out = propagate(Ψ0, genfunc, tlist; storage=storage)
    Ψ_expected = ComplexF64[-1/√2, -1im/√2]  # note the phases

    pop0 = abs.(storage[1,:]).^2
    SHOWPLOT && println(
        lineplot(tlist ./ π, pop0, ylim=[0, 1], title="fw prop")
    )

    @test norm(Ψ_out - Ψ_expected) < 1e-12
    @test pop0[end] ≈ 0.5

    # Propagating backwards in time should exactly reverse the dyanmics (since
    # they are unitary). Thus, we should end up back at the initial state, and
    # the stored states (being filled back-to-front) should exactly match the
    # stored states from the forward propagation.

    storage_bw = init_storage(Ψ0, tlist)
    Ψ_out_bw = propagate(Ψ_out, genfunc, tlist; backwards=true,
                         storage=storage_bw)

    pop0_bw = abs.(storage_bw[1,:]).^2
    SHOWPLOT && println(
        lineplot(tlist ./ π, pop0_bw, ylim=[0, 1], title="bw prop")
    )

    @test norm(Ψ_out_bw - Ψ0) < 1e-12
    @test pop0_bw[1] ≈ 1.0

    @test norm(storage - storage_bw) < 1e-12

end

@testset "Optomech propagation" begin
    include("optomech.jl")
    Ψ0 = ket(0, N_cav) ⊗ ket(2, N_mech)
    H = H_cav + H_mech + H_int
    tlist = [0:0.2:50;]
    Ψ1 = propagate(Ψ0, (tlist, i; kwargs...)->H, tlist, method=:newton)
    @test (norm(Ψ1) - 1.0) < 1e-12
    Ψ2 = propagate(Ψ0, (tlist, i; kwargs...)->H, tlist, method=:cheby)
    @test norm(Ψ1 - Ψ2) < 1e-10
    Ψ3 = propagate(Ψ0, (tlist, i; kwargs...)->H, tlist, method=:auto)
    @test norm(Ψ1 - Ψ3) < 1e-10
end


text:
module SBTestGradient

using SpatialBoundaries
using NeutralLandscapes
using Test

landscape = rand(PlanarGradient(), 40, 50)
x = collect(LinRange(0.2, 1.8, size(landscape, 1)))
y = collect(LinRange(0.2, 0.8, size(landscape, 2)))

W = wombling(x, y, landscape)
@test isa(W, LatticeWomble)
@test size(W.m) == size(landscape) .- 1
@test length(W.x) == length(x) - 1
@test length(W.y) == length(y) - 1

end


text:
struct AssimilationProblem
    """The solution initial condition, i.e. `u(tspan[1])=u₀` """
    u₀
    """The solution `u(t)` will be computed for `tspan[1] ≤ t ≤ tspan[2]`. """
    tspan
    """The underlying ODE is `du = f(u,p,t)` for out-of-place and `f(du,u,p,t)` for in-place."""
    f
    """Constant parameters to be supplied as the second argument of f"""
    p
    """Analysis time step """
    Δt
    """Process noise covariance matrix"""
    Q
    """Initial model error covariance matrix"""
    B₀
    """Observation error covariance matrix"""
    R
    """Observation time step """
    Δt_obs
    """Observation function `h:u\to w` """
    h
    """Jacobian of the observation function, i.e. `Dh = ∂ h_i / ∂ u_j` """
    Dh
    """Observation data"""
    w
end

export AssimilationProblem


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SourcePath(sources, batchsize=10^3)

Traverse a spatial object from `sources` and outwards.
Optionally pass a `batchsize` for KD-tree evaluations.
"""
struct SourcePath <: AbstractPath
  sources::Vector{Int}
  batchsize::Int
end

SourcePath(sources) = SourcePath(sources, 10^3)

function traverse(object, path::SourcePath)
  sources   = path.sources
  batchsize = path.batchsize
  @assert allunique(sources) "non-unique sources"
  @assert all(1 .≤ sources .≤ nelms(object)) "sources must be valid locations"
  @assert length(sources) ≤ nelms(object) "more sources than points in object"

  # coordinate matrix for source points
  S = coordinates(object, sources)

  # fit search tree
  kdtree = KDTree(S)

  # other locations that are not sources
  others = setdiff(1:nelms(object), sources)

  # process points in batches
  batches = Iterators.partition(others, batchsize)

  # pre-allocate memory for coordinates
  X = Matrix{coordtype(object)}(undef, ncoords(object), batchsize)

  # compute distances to sources
  dists = []
  for batch in batches
    coordinates!(X, object, batch)
    _, ds = knn(kdtree, view(X,:,1:length(batch)), length(sources), true)
    append!(dists, ds)
  end

  [sources; view(others, sortperm(dists))]
end


text:
using SparseArrays

@testset "SZigZag (Iter)" begin
    d = 8
    S = 1.3I + 0.5sprandn(d, d, 0.1)
    Γ = S*S'

    ∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x) # sparse computation
    ∇ϕmoving(t, x, θ, i, t′, F, Γ) = ZigZagBoomerang.idot_moving!(Γ, i, t, x, θ, t′, F) # sparse computation

    t0 = 0.0
    x0 = rand(d)
    θ0 = rand([-1.0,-0.5,0.5,1.0], d)

    c = .7*[norm(Γ[:, i], 2) for i in 1:d]

    Z = ZigZag(0.9Γ, x0*0)
    T = 1000.0

    sampler = FactSampler(∇ϕ, t0 =>(x0, θ0), c, Z, Γ)

    tr = trace(sampler, T)
    dt = 0.5
    @test 0.1/sqrt(T) < mean(abs.(mean(tr))) < 2/sqrt(T)

    ts, xs = sep(collect(discretize(tr, dt)))

    @test mean(abs.(mean(xs))) < 2/sqrt(T)
    @test mean(abs.(cov(xs) - inv(Matrix(Γ)))) < 2.5/sqrt(T)
end


text:
function test_hinge_loss_layer(backend::Backend, T, eps)
  println("-- Testing HingeLossLayer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################
  tensor_dim = abs(rand(Int)) % 4 + 2
  dims = tuple((abs(rand(Int,tensor_dim)) % 6 + 6)...)
  println("    > $dims")
  preds = rand(T, dims)*4-2
  labels = round(rand(T, dims))*2-1
  errs_mask = preds.*labels .< one(T)

  ############################################################
  # Setup
  ############################################################
  layer  = HingeLossLayer(; bottoms=[:predictions, :labels])
  pred_blob  = make_blob(backend, T, dims)
  label_blob = make_blob(backend, T, dims)
  diff_blob  = make_blob(backend, T, dims)
#   diff_blob2 = make_blob(backend, T, dims)
  copy!(pred_blob, preds)
  copy!(label_blob, labels)
  inputs = Blob[pred_blob, label_blob]
  diffs  = Blob[diff_blob, NullBlob()]
#   diffs2 = Blob[diff_blob, diff_blob2]

  state  = setup(backend, layer, inputs, diffs)

  forward(backend, state, inputs)

  loss = sum(max(one(T) .- preds.*labels, zero(T))) / dims[end]
  @test -eps < loss-state.loss < eps

  backward(backend, state, inputs, diffs)
  grad  = -labels .* errs_mask / dims[end]
#   grad2 = -preds  .* errs_mask / dims[end]
  diff  = similar(grad)
  copy!(diff, diffs[1])

  @test all(-eps .< grad - diff .< eps)

  # check accumulated/reset statistics
  forward(backend, state, Blob[label_blob, label_blob])
  @test 0 == state.loss
  @test -eps < 0.5loss-state.loss_accum < eps
  reset_statistics(state)
  forward(backend, state, inputs)
  @test -eps < loss-state.loss_accum < eps

  shutdown(backend, state)
end

function test_hinge_loss_layer(backend::Backend)
  test_hinge_loss_layer(backend, Float32, 1e-2)
  test_hinge_loss_layer(backend, Float64, 1e-8)
end

if test_cpu
  test_hinge_loss_layer(backend_cpu)
end
if test_gpu
  test_hinge_loss_layer(backend_gpu)
end


text:
function ispangram(input)
    # https://docs.julialang.org/en/v1/base/strings/#Base.Unicode.lowercase
    # https://docs.julialang.org/en/v1/base/strings/#Base.join
    # https://docs.julialang.org/en/v1/base/collections/#Base.all-Tuple{Any}
    input = lowercase(input)
    alphabet = join('a':'z')

    return all(c in input for c in alphabet)
end

max_stars_repo_name: joaopalmeiro/learning-julia

max_issues_repo_name: joaopalmeiro/learning-julia

max_forks_repo_name: joaopalmeiro/learning-julia


text:
# sparse Nauty graphs

const sg_weight = Cint

struct SparseNautyGraph <: AbstractGraph{Int}
    nv::Cint # number of vertices
    nde::Cint # number of directed edges (1 for loops)
    v::Ptr{Cint} # for each vertex, its neighbour list is at e[v[i]:v[i]+d[i]-1]
    d::Ptr{Cint} # (out)degrees of vertices
    e::Ptr{Cint} # list of neighourhoods, see above
    w::Ptr{sg_weight} # unused, should be NULL
    vlen::Int # allocated length of v
    dlen::Int # ... d
    elen::Int # ... e
    wlen::Int # ... w; should be 0
end


text:
include("test_ascii.jl")
include("test_utf8.jl")


text:
Base.split(dlm) = s -> split(s, dlm, keepempty = false)

shorten(a::Vector{SubString{String}}) = Symbol(a[1][1], a[2][1], a[3][1])

extend_rotation(sym::Symbol) = join([c * "rotation" for c ∈ string(sym)[1:3]], ' ')
extend_position(sym::Symbol) = join([c * "position" for c ∈ string(sym)[1:3]], ' ')

constructor(sym::Symbol) = getfield(BVHFiles, Symbol("Rot", sym))

degrees(R::Rotation) = [getfield(R, θ) |> rad2deg for θ ∈ [Symbol("theta", n) for n ∈ 1:3]]

radians(v::Vector{Float64}, sym::Symbol, constr::Type = constructor(sym)) = constr([deg2rad(θ) for θ ∈ v]...)

rotation(g::BVHGraph, v::Integer, sym::Symbol, f::Integer) = radians(rotations(g, v, f), sym)
rotation(g::BVHGraph, v::Integer, f::Integer) = radians(rotations(g, v, f), sequence(g, v))

function rotation!(g::BVHGraph, v::Integer, f::Integer, rot)
    sym = sequence(g, v)
    constr = constructor(sym)
    rotations(g, v)[f, :] = rot |> constr |> degrees
    return nothing
end

Rx(ψ) = [1.0 0.0 0.0; 0.0 cos(ψ) -sin(ψ); 0.0 sin(ψ) cos(ψ)]
Ry(θ) = [cos(θ) 0.0 sin(θ); 0.0 1.0 0.0; -sin(θ) 0.0 cos(θ)]
Rz(φ) = [cos(φ) -sin(φ) 0.0; sin(φ) cos(φ) 0.0; 0.0 0.0 1.0]

for name in ("Rxyz", "Rxyx", "Rxzy", "Rxzx", "Ryxz", "Ryzx", "Ryxy", "Ryzy", "Rzxy", "Rzyx", "Rzxz", "Rzyz")
    s = name * "(ψ, θ, φ) = " * "R" * name[2] * "(ψ) * " * "R" * name[3] * "(θ) * " * "R" * name[4] * "(φ)"
    eval(Meta.parse(s))
end

function rot(g::BVHGraph, v::Integer, ψ, θ, φ)
    sym = sequence(g, v)
    f = getfield(BVHFiles, Symbol("R", string(sym) |> lowercase))
    return f(ψ, θ, φ)
end

rot(g::BVHGraph, v::Integer, vec) = rot(g, v, vec...)

function rot(sym::Symbol, ψ, θ, φ)
    f = getfield(BVHFiles, Symbol("R", string(sym) |> lowercase))
    return f(ψ, θ, φ)
end

rot(sym::Symbol, vec) = rot(sym, vec...)


text:
@non_differentiable !(::Any)
@non_differentiable *(::Union{Regex, AbstractChar, AbstractString}...)
@non_differentiable -(::AbstractChar, ::AbstractChar)
@non_differentiable <=(::Any)
@non_differentiable <=(::Any, ::Any)
@non_differentiable <(::Any)
@non_differentiable <(::Any, ::Any)
@non_differentiable >=(::Any)
@non_differentiable >=(::Any, ::Any)
@non_differentiable >(::Any)
@non_differentiable >(::Any, ::Any)
@non_differentiable ==(::Any)
@non_differentiable ==(::Any, ::Any)
@non_differentiable ===(::Any, ::Any)

@non_differentiable Bool(::Any)

@non_differentiable accumulate(::Any, ::AbstractArray{Bool})
@non_differentiable accumulate!(::Any, ::Any, ::AbstractArray{Bool})
@non_differentiable cat(::AbstractArray{Bool}...)
@non_differentiable circcopy!(::Any, ::AbstractArray{Bool})
@non_differentiable circshift(::AbstractArray{Bool}, ::Any)
@non_differentiable circshift!(::Any, ::AbstractArray{Bool}, ::Any)
@non_differentiable circshift!(::AbstractArray{Bool}, ::Any)
@non_differentiable conj(::AbstractArray{Bool})
@non_differentiable conj!(::AbstractArray{Bool})
@non_differentiable cumprod(::AbstractArray{Bool})
@non_differentiable cumprod!(::Any, ::AbstractArray{Bool})
@non_differentiable cumsum(::AbstractArray{Bool})
@non_differentiable cumsum!(::Any, ::AbstractArray{Bool})
@non_differentiable diff(::AbstractArray{Bool})
@non_differentiable dropdims(::AbstractArray{Bool})
@non_differentiable eachcol(::AbstractArray{Bool})
@non_differentiable eachrow(::AbstractArray{Bool})
@non_differentiable eachslice(::AbstractArray{Bool})
@non_differentiable fill(::AbstractArray{Bool})
@non_differentiable fill!(::AbstractArray, ::Any)
@non_differentiable findprev(::AbstractArray{Bool}, ::Any)
@non_differentiable findprev(::Any, ::AbstractArray{Bool}, ::Any)
@non_differentiable getindex(::AbstractArray{Bool}, ::Any...)
@non_differentiable hcat(::AbstractArray{Bool}...)
@non_differentiable hvcat(::Any, ::AbstractArray{Bool}...)
@non_differentiable isperm(::AbstractArray{Bool})
@non_differentiable invperm(::AbstractArray{Bool})
@non_differentiable invpermute!(::Any, ::AbstractArray{Bool})
@non_differentiable mapslices(::Any, ::AbstractArray{Bool})
@non_differentiable Matrix(::AbstractArray{Bool})
@non_differentiable ndims(::AbstractArray{Bool})
@non_differentiable parent(::AbstractArray{Bool})
@non_differentiable parentindices(::AbstractArray{Bool})
@non_differentiable permute!(::AbstractArray{Bool}, ::Any)
@non_differentiable permutedims(::AbstractArray{Bool})
@non_differentiable permutedims(::AbstractArray{Bool}, ::Any)
@non_differentiable promote_shape(::AbstractArray{Bool}, ::Any)
@non_differentiable repeat(::AbstractArray{Bool}, ::Any...)
@non_differentiable reshape(::AbstractArray{Bool}, ::Any...)
@non_differentiable reverse(::AbstractArray{Bool}, ::Any...)
@non_differentiable reverse!(::AbstractArray{Bool}, ::Any...)
@non_differentiable reverseind(::AbstractArray{Bool}, ::Any)
@non_differentiable rot180(::AbstractArray{Bool})
@non_differentiable rot180(::AbstractArray{Bool}, ::Any)
@non_differentiable rotl90(::AbstractArray{Bool})
@non_differentiable rotl90(::AbstractArray{Bool}, ::Any)
@non_differentiable rotr90(::AbstractArray{Bool})
@non_differentiable rotr90(::AbstractArray{Bool}, ::Any)
@non_differentiable selectdim(::AbstractArray{Bool}, ::Any, ::Any)
@non_differentiable setindex!(::AbstractArray{Bool}, ::Any, ::Any)
@non_differentiable similar(::AbstractArray{Bool}, ::Any...)
@non_differentiable stride(::AbstractArray{Bool}, ::Any)
@non_differentiable strides(::AbstractArray{Bool})
@non_differentiable vcat(::AbstractArray{Bool}...)
@non_differentiable vec(::AbstractArray{Bool})
@non_differentiable Vector(::AbstractArray{Bool})

@non_differentiable Array(::AbstractArray{Bool})
@non_differentiable IndexStyle(::AbstractArray{Bool})

@non_differentiable abspath(::AbstractString...)
@non_differentiable all(::Any)
@non_differentiable all(::Any, ::Any)
@non_differentiable any(::Any)
@non_differentiable any(::Any, ::Any)
@non_differentiable ascii(::AbstractString)
@non_differentiable axes(::Any)
@non_differentiable axes(::Any, ::Any)

@non_differentiable basename(::AbstractString)
@non_differentiable bytesavailable(::Any)

@non_differentiable cd(::AbstractString)
@non_differentiable chomp(::AbstractString)
@non_differentiable chop(::AbstractString)
@non_differentiable close(::Any)
@non_differentiable cmp(::AbstractString, ::AbstractString)
@non_differentiable cmp(::Symbol, ::Symbol)
@non_differentiable codepoint(::Char)
@non_differentiable codeunit(::AbstractString)
@non_differentiable codeunit(::SubString)
@non_differentiable codeunits(::AbstractString)
@non_differentiable copy(::Base.GenericIOBuffer)
@non_differentiable count(::Union{Regex, AbstractString}, ::AbstractString)
@non_differentiable countlines(::AbstractString)
@non_differentiable countlines(::IO)
@non_differentiable cp(::AbstractString, ::AbstractString)

@non_differentiable detach(::Cmd)
@non_differentiable dirname(::AbstractString)
@non_differentiable displayable(::AbstractDisplay, ::AbstractString)
@non_differentiable displayable(::AbstractString)
@non_differentiable displaysize(::IO)
@non_differentiable download(::AbstractString)
@non_differentiable download(::AbstractString, ::AbstractString)

@non_differentiable eachindex(::Any...)
@non_differentiable eachline(::AbstractString)
@non_differentiable eachline(::IO)
@non_differentiable eachmatch(::Regex, ::AbstractString)
@non_differentiable endswith(::AbstractString, ::AbstractString)
@non_differentiable endswith(::AbstractString, ::Regex)
@non_differentiable eof(::Any)
@non_differentiable error(::AbstractString)
@non_differentiable escape_string(::AbstractString)
@non_differentiable escape_string(::IO, ::AbstractString)
@non_differentiable evalfile(::AbstractString)
@non_differentiable evalfile(::AbstractString, ::Array{String,1})
@non_differentiable expanduser(::AbstractString)

@non_differentiable fd(::Base.Filesystem.File)
@non_differentiable fd(::IOStream)
@non_differentiable fieldtype(T, ::Union{Symbol, Integer})
VERSION >= v"1.1" && @non_differentiable fieldtypes(T)
@non_differentiable fieldname(T, ::Integer)
@non_differentiable fieldnames(T)

@non_differentiable findall(::Union{Regex, AbstractString, Function}, ::AbstractString)
@non_differentiable findall(::Function, ::AbstractArray)
@non_differentiable findall(::AbstractArray)
@non_differentiable findfirst(::Union{Regex, AbstractString, Function, AbstractChar}, ::AbstractString)
@non_differentiable findfirst(::Function, ::Any)
@non_differentiable findfirst(::AbstractArray)
@non_differentiable findlast(::Union{AbstractString, Function, AbstractChar}, ::AbstractString)
@non_differentiable findlast(::Function, ::Any)
@non_differentiable findlast(::AbstractArray)
@non_differentiable first(::AbstractChar)
@non_differentiable first(::Cmd)
@non_differentiable firstindex(::AbstractChar)
@non_differentiable firstindex(::AbstractString)
@non_differentiable firstindex(::Cmd)
@non_differentiable flush(::Any)

@non_differentiable gensym(::Symbol)
@non_differentiable gensym(::String...)
@non_differentiable getindex(::AbstractChar)
@non_differentiable getindex(::AbstractString, ::Colon)
@non_differentiable getindex(::Base.TTY, ::Symbol)
@non_differentiable getindex(::RegexMatch, ::AbstractString)
@non_differentiable getindex(::RegexMatch, ::Symbol)
@non_differentiable getpid(::Base.Process)

@non_differentiable haskey(::Any, ::Any)
@non_differentiable hash(::Any)
@non_differentiable hash(::Any, ::UInt)
@non_differentiable hex2bytes(::AbstractString)

@non_differentiable ignorestatus(::Cmd)
@non_differentiable in(::Any, ::Any)
@non_differentiable include_dependency(::AbstractString)
@non_differentiable isa(::Any, ::Any)
@non_differentiable isabspath(::AbstractString)
@non_differentiable isabstracttype(::Any)
@non_differentiable isapprox(::Any, ::Any)
@non_differentiable isascii(::AbstractChar)
@non_differentiable isascii(::AbstractString)
@non_differentiable isassigned(::Any...)
@non_differentiable isbits(::Any)
@non_differentiable isbitstype(::Any)
@non_differentiable isblockdev(::Any...)
@non_differentiable ischardev(::Any...)
@non_differentiable iscntrl(::AbstractChar)
@non_differentiable isconcretetype(::Any)
@non_differentiable isconst(::Module, ::Symbol)
@non_differentiable isdigit(::AbstractChar)
@non_differentiable isdir(::Any...)
@non_differentiable isdirpath(::AbstractString)
VERSION >= v"1.5" && @non_differentiable isdisjoint(::Any, ::Any)
@non_differentiable isdispatchtuple(::Any)
@non_differentiable isempty(::Any)
@non_differentiable isequal(::Any)
@non_differentiable isequal(::Any, ::Any)
@non_differentiable iseven(::Any)
@non_differentiable isfifo(::Any...)
@non_differentiable isfile(::Any...)
@non_differentiable isfinite(::Any)
@non_differentiable isimmutable(::Any)
@non_differentiable isinf(::Any)
@non_differentiable isinteger(::Any)
@non_differentiable isless(::Any, ::Any)
@non_differentiable isletter(::AbstractChar)
@non_differentiable islink(::Any...)
@non_differentiable islowercase(::AbstractChar)
@non_differentiable ismarked(::IO)
@non_differentiable ismissing(::Any)
@non_differentiable ismount(::Any...)
VERSION >= v"1.5" && @non_differentiable ismutable(::Any)
@non_differentiable isnan(::Any)
VERSION >= v"1.1" && @non_differentiable isnothing(::Any)
@non_differentiable isnumeric(::AbstractChar)
@non_differentiable isodd(::Any)
@non_differentiable isone(::Any)
@non_differentiable isopen(::Any)
@non_differentiable ispath(::Any...)
@non_differentiable isperm(::Any)
@non_differentiable ispow2(::Any)
@non_differentiable isprimitivetype(::Any)
@non_differentiable isprint(::AbstractChar)
@non_differentiable ispunct(::AbstractChar)
@non_differentiable ispunct(::Any)
@non_differentiable isreadable(::Any)
@non_differentiable isreadonly(::Any)
@non_differentiable isready(::Any)
@non_differentiable isreal(::Any)
@non_differentiable issetequal(::Any, ::Any)
@non_differentiable issetgid(::Any)
@non_differentiable issetuid(::Any)
@non_differentiable issocket(::Any)
@non_differentiable issorted(::Any)
@non_differentiable issorted(::Any, ::Base.Order.Ordering)
@non_differentiable isspace(::AbstractChar)
@non_differentiable issticky(::Any)
@non_differentiable isstructtype(::Any)
@non_differentiable issubnormal(::Any)
@non_differentiable issubset(::Any, ::Any)
@non_differentiable istaskdone(::Task)
VERSION >= v"1.3" && @non_differentiable istaskfailed(::Task)
@non_differentiable istaskstarted(::Task)
@non_differentiable istextmime(::AbstractString)
@non_differentiable isuppercase(::AbstractChar)
@non_differentiable isvalid(::AbstractChar)
@non_differentiable isvalid(::String)
@non_differentiable isvalid(::SubString{String})
@non_differentiable iswritable(::Any)
@non_differentiable isxdigit(::AbstractChar)
@non_differentiable iszero(::Any)
@non_differentiable iterate(::AbstractChar)
@non_differentiable iterate(::Cmd)
@non_differentiable iterate(::AbstractString)

@non_differentiable join(::Any)
@non_differentiable join(::Any, ::Any)
@non_differentiable join(::Any, ::Any, ::Any)
@non_differentiable join(::IO, ::Any, ::Any, ::Any)
@non_differentiable joinpath(::AbstractString, ::AbstractString...)

@non_differentiable keys(::AbstractString)
@non_differentiable keys(::Cmd)
@non_differentiable kill(::Array{Base.Process,1})
@non_differentiable kill(::Base.Process)
@non_differentiable kill(::Base.ProcessChain)

@non_differentiable last(::AbstractChar)
@non_differentiable last(::Cmd)
@non_differentiable lastindex(::AbstractChar)
@non_differentiable lastindex(::AbstractString)
@non_differentiable lastindex(::Cmd)
@non_differentiable length(::Any)
@non_differentiable lowercase(::AbstractString)
@non_differentiable lowercase(::AbstractChar)
@non_differentiable lowercasefirst(::AbstractString)
@non_differentiable lstat(::AbstractString)
@non_differentiable lstrip(::AbstractString)

@non_differentiable mark(::IO)
@non_differentiable match(::Regex, ::AbstractString)
@non_differentiable mkdir(::AbstractString)
@non_differentiable mkpath(::AbstractString)
@non_differentiable mktemp(::AbstractString)
@non_differentiable mktempdir(::AbstractString)
@non_differentiable mv(::AbstractString, ::AbstractString)

@non_differentiable ncodeunits(::AbstractChar)
@non_differentiable ncodeunits(::AbstractString)
@non_differentiable ndims(::AbstractChar)
@non_differentiable nfields(::Any)
@non_differentiable normpath(::AbstractString...)

@non_differentiable occursin(::Regex, ::AbstractString)
@non_differentiable occursin(::Union{AbstractChar, AbstractString}, ::AbstractString)
VERSION >= v"1.4" && @non_differentiable only(::Char)
@non_differentiable open(::Any)

@non_differentiable pointer(::Any)
@non_differentiable popdisplay(::AbstractDisplay)
@non_differentiable position(::Base.Filesystem.File)
@non_differentiable position(::Base.GenericIOBuffer)
@non_differentiable position(::Base.SecretBuffer)
@non_differentiable position(::IOStream)
@non_differentiable print(::Any...)
@non_differentiable println(::Any...)
@non_differentiable process_exited(::Any)
@non_differentiable process_running(::Any)
@non_differentiable pushdisplay(::AbstractDisplay)

@non_differentiable read(::IO)
@non_differentiable readavailable(::Any)
@non_differentiable readdir(::AbstractString)
@non_differentiable readline(::AbstractString)
@non_differentiable readline(::IO)
@non_differentiable readlines(::AbstractString)
@non_differentiable readlink(::AbstractString)
@non_differentiable readuntil(::IO, ::AbstractChar)
@non_differentiable readuntil(::IO, ::AbstractString)
@non_differentiable realpath(::AbstractString)
if isdefined(Base, :redirect_stdio)
    @non_differentiable (::Base.redirect_stdio)(
        ::Union{IOStream, Base.LibuvStream, Base.DevNull, Base.AbstractPipe},
    )
else
    @non_differentiable redirect_stderr(
        ::Union{IOStream, Base.LibuvStream, Base.DevNull, IOContext},
    )
    @non_differentiable redirect_stdin(
        ::Union{IOStream, Base.LibuvStream, Base.DevNull, IOContext},
    )
    @non_differentiable redirect_stdout(
        ::Union{IOStream, Base.LibuvStream, Base.DevNull, IOContext},
    )
end
@non_differentiable relpath(::AbstractString, ::AbstractString)
@non_differentiable relpath(::String)
@non_differentiable repr(::Any)
@non_differentiable repr(::MIME, ::Any)
@non_differentiable repr(::AbstractString, ::Any)
@non_differentiable reset(::IO)
@non_differentiable reverse(::AbstractString)
@non_differentiable rm(::AbstractString)
@non_differentiable rsplit(::AbstractString)
@non_differentiable rsplit(::AbstractString, ::AbstractChar)
@non_differentiable rstrip(::AbstractString)

@non_differentiable schedule(::Task)
@non_differentiable schedule(::Task, ::Any)
@non_differentiable seekend(::Any)
@non_differentiable seekstart(::Any)
@non_differentiable setenv(::Cmd)
@non_differentiable show(::Any)
@non_differentiable show(::Any, ::Any)
@non_differentiable showerror(::IO, ::Exception)
@non_differentiable size(::Any)
@non_differentiable size(::Any, ::Any)
@non_differentiable sizeof(::Any)
@non_differentiable sleep(::Any)
@non_differentiable split(::AbstractString)
@non_differentiable split(::AbstractString, ::AbstractChar)
@non_differentiable splitdir(::AbstractString)
@non_differentiable splitdrive(::AbstractString)
@non_differentiable splitext(::AbstractString)
@non_differentiable startswith(::AbstractString, ::AbstractString)
VERSION >= v"1.1" && @non_differentiable splitpath(::AbstractString)
@non_differentiable startswith(::AbstractString, ::Regex)
@non_differentiable stat(::AbstractString)
@non_differentiable stat(::Base.Filesystem.File)
@non_differentiable stat(::IOStream)
@non_differentiable stat(::RawFD)
@non_differentiable string(::Any...)
@non_differentiable strip(::AbstractString)
@non_differentiable success(::Array{Base.Process,1})
@non_differentiable success(::Base.Process)
@non_differentiable success(::Base.ProcessChain)
@non_differentiable supertype(::Any)
@non_differentiable symlink(::AbstractString, ::AbstractString)

@non_differentiable take!(::Base.GenericIOBuffer)
@non_differentiable take!(::IOStream)
@non_differentiable tempname(::AbstractString)
@non_differentiable textwidth(::AbstractChar)
@non_differentiable textwidth(::AbstractString)
@non_differentiable titlecase(::AbstractString)
@non_differentiable titlecase(::AbstractChar)
@non_differentiable touch(::AbstractString)
@non_differentiable typemin(::String)

@non_differentiable unescape_string(::AbstractString)
@non_differentiable unescape_string(::IO, ::AbstractString)
@non_differentiable unmark(::IO)
@non_differentiable unsafe_string(::Cstring)
@non_differentiable uppercase(::AbstractString)
@non_differentiable uppercase(::AbstractChar)
@non_differentiable uppercasefirst(::AbstractString)

@non_differentiable wait(::Base.Process)
@non_differentiable wait(::Base.ProcessChain)
@non_differentiable wait(::RawFD)
@non_differentiable write(::IO, ::Any...)

@non_differentiable xor(::Any...)

# Non-public Base
@non_differentiable Base.gc_num()
@non_differentiable Base.time_ns()
@non_differentiable Base.typename(::Any)

@non_differentiable Broadcast.combine_styles(::Any...)
@non_differentiable Broadcast.result_style(::Any)
@non_differentiable Broadcast.result_style(::Any, ::Any)

@non_differentiable Libc.free(::Any)
@non_differentiable Libc.getpid()
@non_differentiable Libc.strptime(::AbstractString)
@non_differentiable Libc.strptime(::AbstractString, ::AbstractString)

@non_differentiable Meta.parse(str::AbstractString)
@non_differentiable Meta.parse(str::AbstractString, pos::Integer)

@non_differentiable Sys.cpu_summary(::IO)
@non_differentiable Sys.isapple(::Symbol)
@non_differentiable Sys.isbsd(::Symbol)
VERSION >= v"1.1" && @non_differentiable Sys.isdragonfly(::Symbol)
@non_differentiable Sys.isexecutable(::AbstractString)
VERSION >= v"1.1" && @non_differentiable Sys.isfreebsd(::Symbol)
VERSION >= v"1.2" && @non_differentiable Sys.isjsvm(::Symbol)
@non_differentiable Sys.islinux(::Symbol)
VERSION >= v"1.1" && @non_differentiable Sys.isnetbsd(::Symbol)
VERSION >= v"1.1" && @non_differentiable Sys.isopenbsd(::Symbol)
@non_differentiable Sys.isunix(::Symbol)
@non_differentiable Sys.iswindows(::Symbol)
@non_differentiable Sys.which(::AbstractString)

@non_differentiable Threads.nthreads()
@non_differentiable Threads.threadid()
@non_differentiable Threads.threadid(::Task)


text:
using Compat.Test
using Tries

import Compat: Nothing

tests = ["trie"]

if length(ARGS) > 0
    tests = ARGS
end

@testset "Tries" begin

for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    println("$fp ...")
    include(fp)
end

end # @testset


text:
abstract type Message
end

struct DZip{S,T}
  X::S
  Y::T
end

struct Start
end
struct Close{T}
    state::T
end
struct Postprocess{T,S} <: Message
    state::T
    signal::S
end
Base.iterate(M::Message) = getfield(M, 1), 1
Base.iterate(M::Message, Any) = getfield(M, 2), nothing
Base.iterate(M::Message, ::Nothing) = nothing

struct DRange{T}
    start::T
    stop::T
end
dyniterate(r::DRange, ::Start) = (r.start, r.start)
dyniterate(r::DRange, s) = (s == r.stop) ? nothing : (s+1, s+1)
dyniterate(r::DRange, E::Close) = nothing

Base.iterate(Z::DZip) = dyniterate(Z, Start())
Base.iterate(Z::DZip, state) = dyniterate(Z, state)

dyniterate(Z::DZip, ::Start) = dyniterate(Z::DZip, (Start(), Start()))
dyniterate(Z::DZip, (start, signal)::Postprocess{Start}) = dyniterate(Z::DZip, Postprocess((Start(), Start()), signal))

function dyniterate(Z::DZip, (p, q))
    ϕ = dyniterate(Z.X, p)
    ϕ === nothing && return nothing
    x, p = ϕ

    ϕ = dyniterate(Z.Y, q)
    ϕ === nothing && return nothing
    y, q = ϕ

    (x, y), (p, q)
end

function dyniterate(Z::DZip, ((p,q), signal)::Postprocess)
    ϕ = dyniterate(Z.X, p)
    ϕ === nothing && return dyniterate(Z.Y, signal(q))
    x, p = ϕ

    ϕ = dyniterate(Z.Y, q)
    ϕ === nothing && return dyniterate(Z.X, signal(p))
    y, q = ϕ

    (x, y), Postprocess((p, q), signal)
end

using BenchmarkTools

function f(n)
    Z = zip(1:n, 1:n)
    i = 0
    for (x, y) in Z
        i += x + 2y % 100
    end
    i
end

function g(n)
    Z = DZip(DRange(1,n), DRange(1,n))
    i = 0

    for (x, y) in Z
        i += x + 2y % 100
    end
    i
end

function h(n)
    Z = DZip(DRange(1,n), DRange(1,n))
    i = 0
    ϕ = dyniterate(Z, Postprocess(Start(), x->Close(x)))
    while !(ϕ === nothing)
        (x, y), state = ϕ
        i += x + 2y % 100
        ϕ = dyniterate(Z, state)
        ϕ === nothing && return i
    end
    i
end

@btime f(10000)
@btime g(10000)
@btime h(10000)


text:
function weno5_js(epsilon)
       return Weno5(-2, 2, epsilon)
end

function flux_positive_weno5(scheme::Weno5,
                             vector_in::Array{<:AbstractFloat, 1})
    is_w1 = 13. / 12
    is_w2 = 0.25
    tmp1 = [1.; -2.;  1.]
    tmp2 = [1.; -4.;  3.]
    tmp3 = [1.;  0.; -1.]
    tmp4 = [3.; -4.;  1.]
    is_0 = is_w1 * (dot(tmp1, vector_in[1:3])) ^ 2 +
           is_w2 * (dot(tmp2, vector_in[1:3])) ^ 2
    is_1 = is_w1 * (dot(tmp1, vector_in[2:4])) ^ 2 +
           is_w2 * (dot(tmp3, vector_in[2:4])) ^ 2
    is_2 = is_w1 * (dot(tmp1, vector_in[3:5])) ^ 2 +
           is_w2 * (dot(tmp4, vector_in[3:5])) ^ 2
    a_0 = 0.1 / (scheme.epsilon + is_0)
    a_1 = 0.6 / (scheme.epsilon + is_1)
    a_2 = 0.3 / (scheme.epsilon + is_2)
    f_0 = (vector_in[1] * 2 - vector_in[2] * 7 + vector_in[3] * 11) / 6
    f_1 = (- vector_in[2] + vector_in[3] * 5 + vector_in[4] * 2) / 6
    f_2 = (vector_in[3] * 2 + vector_in[4] * 5 - vector_in[5]) / 6
    return (a_0 * f_0 + a_1 * f_1 + a_2 * f_2) / (a_0 + a_1 + a_2)
end # flux_positive_weno5


text:
# Types for molecular models, their particles and trajectories

struct MolecularModelHeader <: Header
	D::Dict{Symbol,Any}

	MolecularModelHeader() = new(Dict{Symbol,Any}())
end

function MolecularModelHeader(header::MolecularModelHeader)
	header2 = MolecularModelHeader()
	merge!(header2, header)
end

Headers.headerval(::MolecularModelHeader, ::Val{:title}, v::AbstractString) =
		String(v)

Headers.headerval(::MolecularModelHeader, ::Val{:step}, v::Integer) = Int(v)

Headers.headerval(::MolecularModelHeader, ::Val{:time}, v::Real) = Float64(v)

Headers.headerval(::MolecularModelHeader, ::Val{:lambda}, v::Real) = Float64(v)

Headers.headerval(::MolecularModelHeader, ::Val{:cell}, v::TriclinicPBC) =
		TriclinicCell(v)

Headers.headerval(::MolecularModelHeader, ::Val{:pressure},
		v::AbstractMatrix{<:Real}) = Basis3D(v)

Headers.headerval(::MolecularModelHeader, ::Val{:virial},
		v::AbstractMatrix{<:Real}) = Basis3D(v)

mutable struct MolecularModel <: Multicollection
	header::MolecularModelHeader
	n::Int
	D::Dict{Symbol,Any}

	MolecularModel(n::Integer = 0) =
			new(MolecularModelHeader(), n, Dict{Symbol,Any}())
end

const MolecularModelView = MulticollectionView{MolecularModel}

const ParticleCollection = Union{MolecularModel,MolecularModelView}

const Particle = MulticollectionItem{MolecularModel}

abstract type MolecularTrajectory <: FormattedStream{MolecularModel} end

function MolecularModel(model::ParticleCollection)
	model2 = MolecularModel(length(model))
	for (key, val) in pairs(model)
		get!(model2, key, val)
	end
	merge!(model2.header, model.header)
	model2
end

MolecularModel(p::Particle) = MolecularModel(view(parent(p), parentindices(p)))

Base.similar(model::ParticleCollection, n::Integer) = MolecularModel(n)

function Multicollections.collval(::MolecularModel, ::Val{:topology},
		n::Integer, v::AbstractGraph)
	@boundscheck length(v) == n || error("expected $(n)-index graph")
	merge!(Graph(n), v)
end

Multicollections.collval(::MolecularModel, ::Val{:topology}, n::Integer,
		V::AbstractVector) = pair!(Graph(n), V...)

Multicollections.collval(::MolecularModel, ::Val{:ids}, n::Integer, v) =
		(Vector{Int}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:ids}, n::Integer,
		::UndefInitializer) = Vector{Int}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:names}, n::Integer, v) =
		(Vector{String}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:names}, n::Integer,
		::UndefInitializer) = Vector{String}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:resids}, n::Integer, v) =
		(Vector{Int}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:resids}, n::Integer,
		::UndefInitializer) = Vector{Int}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:resnames}, n::Integer, v) =
		(Vector{String}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:resnames}, n::Integer,
		::UndefInitializer) = Vector{String}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:chainids}, n::Integer, v) =
		(Vector{String}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:chainids}, n::Integer,
		::UndefInitializer) = Vector{String}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:elements}, n::Integer, v) =
		(Vector{String}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:elements}, n::Integer,
		::UndefInitializer) = Vector{String}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:R}, n::Integer, v) =
		(Vector{Vector3D}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:R}, n::Integer,
		::UndefInitializer) = Vector{Vector3D}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:V}, n::Integer, v) =
		(Vector{Vector3D}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:V}, n::Integer,
		::UndefInitializer) = Vector{Vector3D}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:F}, n::Integer, v) =
		(Vector{Vector3D}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:F}, n::Integer,
		::UndefInitializer) = Vector{Vector3D}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:masses}, n::Integer, v) =
		(Vector{Float64}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:masses}, n::Integer,
		::UndefInitializer) = Vector{Float64}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:charges}, n::Integer, v) =
		(Vector{Float64}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:charges}, n::Integer,
		::UndefInitializer) = Vector{Float64}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:bfactors}, n::Integer, v) =
		(Vector{Float64}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:bfactors}, n::Integer,
		::UndefInitializer) = Vector{Float64}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:occupancies}, n::Integer, v) =
		(Vector{Float64}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:occupancies}, n::Integer,
		::UndefInitializer) = Vector{Float64}(undef, n)

Multicollections.collval(::MolecularModel, ::Val{:SS}, n::Integer, v) =
		(Vector{String}(undef, n) .= v)

Multicollections.collval(::MolecularModel, ::Val{:SS}, n::Integer,
		::UndefInitializer) = Vector{String}(undef, n)

Multicollections.itemtocollprop(::Particle, ::Val{:id}) = :ids
Multicollections.colltoitemprop(::Particle, ::Val{:ids}) = :id

Multicollections.itemtocollprop(::Particle, ::Val{:name}) = :names
Multicollections.colltoitemprop(::Particle, ::Val{:names}) = :name

Multicollections.itemtocollprop(::Particle, ::Val{:resid}) = :resids
Multicollections.colltoitemprop(::Particle, ::Val{:resids}) = :resid

Multicollections.itemtocollprop(::Particle, ::Val{:resname}) = :resnames
Multicollections.colltoitemprop(::Particle, ::Val{:resnames}) = :resname

Multicollections.itemtocollprop(::Particle, ::Val{:chainid}) = :chainids
Multicollections.colltoitemprop(::Particle, ::Val{:chainids}) = :chainid

Multicollections.itemtocollprop(::Particle, ::Val{:element}) = :elements
Multicollections.colltoitemprop(::Particle, ::Val{:elements}) = :element

Multicollections.itemtocollprop(::Particle, ::Val{:mass}) = :masses
Multicollections.colltoitemprop(::Particle, ::Val{:masses}) = :mass

Multicollections.itemtocollprop(::Particle, ::Val{:charge}) = :charges
Multicollections.colltoitemprop(::Particle, ::Val{:charges}) = :charge

Multicollections.itemtocollprop(::Particle, ::Val{:bfactor}) = :bfactors
Multicollections.colltoitemprop(::Particle, ::Val{:bfactors}) = :bfactor

Multicollections.itemtocollprop(::Particle, ::Val{:occupancy}) = :occupancies
Multicollections.colltoitemprop(::Particle, ::Val{:occupancies}) = :occupancy

Multicollections.itemtocollprop(::Particle, ::Val{:ss}) = :SS
Multicollections.colltoitemprop(::Particle, ::Val{:SS}) = :ss

Multicollections.collitemname(::MolecularModel) = "particle"

function pbcstrategy(model::ParticleCollection)
	cell = pbccell(get(model.header, :cell, nothing))
	R = model.R
	Rw, Kw = pbcpos(R, cell)
	cell, Rw, Kw, R
end

function pbcstrategy!(Rw::AbstractVector{Vector3D},
		Kw::AbstractVector{Vector3D}, model::ParticleCollection)
	cell = pbccell(get(model.header, :cell, nothing))
	R = model.R
	pbcpos!(Rw, Kw, R, cell)
	cell, Rw, Kw, R
end


text:
push!(LOAD_PATH, "/mnt/Data/Documents/Agents")
#push!(LOAD_PATH, "D:\\projects\\Agents.jl\\")

using Documenter, Agents

makedocs(
  sitename="Agents.jl",
    pages = [
        "Introduction" => "index.md",
        "Tutorial" => "tutorial.md",
        "Examples" => [
          "Boltzmann wealth distribution" => "boltzmann_example01.md",
          "Forest fire" => "forest_fire.md",
          "Cellular Automata" => "CA.md"
        ],
        "Built-in funtions" => "builtin_functions.md",
        "Comparison against Mesa" => "mesa.md"
    ]

)

deploydocs(
	    repo = "github.com/kavir1698/Agents.jl.git",
   )