text:
# Store the rank with the value, necessary for collecting values in order
struct pval{T}
    rank :: Int
    errorstatus :: Bool
    value :: T
end
pval{T}(p::pval{T}) where {T} = p
pval{T}(p::pval) where {T} = pval{T}(p.rank, p.errorstatus, convert(T, p.value))

errorpval(rank) = pval(rank, true, nothing)

errorstatus(p::pval) = p.errorstatus

# Function to obtain the value of pval types
value(p::pval) = p.value
value(p::Any) = p

struct Product{I}
    iterators :: I
end

getiterators(p::Product) = p.iterators

Base.length(p::Product) = length(Iterators.product(p.iterators...))
Base.iterate(p::Product, st...) = iterate(Iterators.product(p.iterators...), st...)

function product(iter...)
    any(x -> x isa Product, iter) && throw(ArgumentError("the iterators should not be Products"))
    Product(iter)
end

struct Hold{T}
    iterators :: T
end

getiterators(h::Hold) = getiterators(h.iterators)

Base.length(h::Hold) = length(h.iterators)

function check_knownsize(iterator)
    itsz = Base.IteratorSize(iterator)
    itsz isa Base.HasLength || itsz isa Base.HasShape
end

struct ZipSplit{Z, I}
    z :: Z
    it :: I
    skip :: Int
    N :: Int
end

# This constructor differs from zipsplit, as it uses skipped and retained elements
# and not p and np. This type is added to increase compatibility with SplittablesBase
function ZipSplit(itzip, skipped_elements::Integer, elements_on_proc::Integer)
    it = Iterators.take(Iterators.drop(itzip, skipped_elements), elements_on_proc)
    ZipSplit{typeof(itzip), typeof(it)}(itzip, it, skipped_elements, elements_on_proc)
end

Base.length(zs::ZipSplit) = length(zs.it)
Base.eltype(zs::ZipSplit) = eltype(zs.it)
Base.iterate(z::ZipSplit, i...) = iterate(takedrop(z), i...)
takedrop(zs::ZipSplit) = zs.it

function SplittablesBase.halve(zs::ZipSplit)
    nleft = zs.N ÷ 2
    ZipSplit(zs.z, zs.skip, nleft), ZipSplit(zs.z, zs.skip + nleft, zs.N - nleft)
end

zipsplit(iterators::Tuple, np::Integer, p::Integer) = zipsplit(zip(iterators...), np, p)

function zipsplit(itzip::Iterators.Zip, np::Integer, p::Integer)
    check_knownsize(itzip)
    d,r = divrem(length(itzip), np)
    skipped_elements = d*(p-1) + min(r,p-1)
    lastind = d*p + min(r,p)
    elements_on_proc = lastind - skipped_elements
    ZipSplit(itzip, skipped_elements, elements_on_proc)
end

_split_iterators(iterators, np, p) = (zipsplit(iterators, np, p),)
function _split_iterators(iterators::Tuple{Hold{<:Product}}, np, p)
    it_hold = first(iterators)
    (ProductSplit(getiterators(it_hold), np, p), )
end

############################################################################################
# Local mapreduce
############################################################################################

struct NoSplat <: Function
    f :: Function
end
NoSplat(u::NoSplat) = u

_maybesplat(f) = Base.splat(f)
_maybesplat(f::NoSplat) = f

_mapreduce(f, op, iterators...; reducekw...) = mapreduce(f, op, iterators...; reducekw...)
function _mapreduce(fun::NoSplat, op, iterators...; reducekw...)
    mapval = fun.f(iterators...)
    reduce(op, (mapval,); reducekw...)
end

function mapreducenode(f, op, rank, pipe::BranchChannel, selfoutchannel, iterators...; reducekw...)
    # Evaluate the function
    # No communication with other nodes happens here
    try
        fmap = _maybesplat(f)
        if rank == 1
            res = _mapreduce(fmap, op, iterators...; reducekw...)
        else
            # init should only be used once on the first rank
            # remove it from the kwargs on other workers
            kwdict = Dict(reducekw)
            pop!(kwdict, :init, nothing)
            res = _mapreduce(fmap, op, iterators...; kwdict...)
        end
        val = pval(rank, false, res)
        put!(selfoutchannel, val)
    catch
        put!(selfoutchannel, errorpval(rank))
        rethrow()
    end
end

############################################################################################
# Reduction across workers
############################################################################################

abstract type ReductionNode end
struct TopTreeNode <: ReductionNode
    rank :: Int
end
struct SubTreeNode <: ReductionNode
    rank :: Int
end

_maybesort(op::Commutative, vals) = vals
_maybesort(op, vals) = sort!(vals, by = pv -> pv.rank)

function reducechannel(op, c, N; reducekw...)
    vals = [take!(c) for i = 1:N]
    vals = _maybesort(op, vals)
    v = [value(v) for v in vals]
    reduce(op, v; reducekw...)
end

seterrorflag(c, val) = put!(c, take!(c) | val)

function reducedvalue(op, node::SubTreeNode, pipe::BranchChannel, selfoutchannel; reducekw...)
    rank = node.rank

    N = nchildren(pipe) + 1
    err_ch = Channel{Bool}(1)
    put!(err_ch, false)

    self_pval = take!(selfoutchannel)
    if errorstatus(self_pval)
        return errorpval(rank)
    else
        put!(selfoutchannel, self_pval)
    end

    @sync for i = 1:nchildren(pipe)
        @async begin
            child_pval = take!(pipe.childrenchannel)
            if errorstatus(child_pval)
                seterrorflag(err_ch, true)
            else
                put!(selfoutchannel, child_pval)
                seterrorflag(err_ch, false)
            end
        end
    end

    take!(err_ch) && return errorpval(rank)

    redval = reducechannel(op, selfoutchannel, N; reducekw...)

    return pval(rank, false, redval)
end
function reducedvalue(op, node::TopTreeNode, pipe::BranchChannel, ::Any; reducekw...)
    rank = node.rank

    N = nchildren(pipe)
    c = Channel(N)
    err_ch = Channel{Bool}(1)
    put!(err_ch, false)

    @sync for i in 1:N
        @async begin
            child_pval = take!(pipe.childrenchannel)
            if errorstatus(child_pval)
                seterrorflag(err_ch, true)
            else
                put!(c, child_pval)
                seterrorflag(err_ch, false)
            end
        end
    end

    take!(err_ch) && return errorpval(rank)

    redval = reducechannel(op, c, N; reducekw...)

    return pval(rank, false, redval)
end

function reducenode(op, node::ReductionNode, pipe::BranchChannel, selfoutchannel = nothing; kwargs...)
    # This function that communicates with the parent and children
    rank = node.rank
    try
        kwdict = Dict(kwargs)
        pop!(kwdict, :init, nothing)
        res = reducedvalue(op, node, pipe, selfoutchannel; kwdict...)
        put!(pipe.parentchannel, res)
    catch
        put!(pipe.parentchannel, errorpval(rank))
        rethrow()
    finally
        GC.gc()
    end

    return nothing
end

function pmapreduceworkers(f, op, tree_branches, iterators; reducekw...)

    tree, branches = tree_branches

    nworkerstree = nworkers(tree)
    extrareducenodes = length(tree) - nworkerstree

    @sync for (ind, mypipe) in enumerate(branches)
        p = mypipe.p
        ind_reduced = ind - extrareducenodes
        rank = ind_reduced
        if ind_reduced > 0
            iterable_on_proc = _split_iterators(iterators, nworkerstree, rank)
            @spawnat p begin
                selfoutchannel = Channel(nchildren(mypipe) + 1)
                @sync begin
                    @async mapreducenode(f, op, rank, mypipe, selfoutchannel, iterable_on_proc...; reducekw...)
                    @async reducenode(op, SubTreeNode(rank), mypipe, selfoutchannel; reducekw...)
                end
            end
        else
            @spawnat p reducenode(op, TopTreeNode(rank), mypipe; reducekw...)
        end
    end

    tb = topbranch(tree, branches)
    value(take!(tb.parentchannel))
end

"""
    pmapreduce(f, op, [pool::AbstractWorkerPool], iterators...; reducekw...)

Evaluate a parallel `mapreduce` over the elements from `iterators`.
For multiple iterators, apply `f` elementwise.

The keyword arguments `reducekw` are passed on to the reduction.

See also: [`pmapreduce_productsplit`](@ref)
"""
function pmapreduce(f, op, pool::AbstractWorkerPool, iterators...; reducekw...)
    N = length(zip(iterators...))

    if N <= 1 || nworkers(pool) == 1
        iterable_on_proc = _split_iterators(iterators, 1, 1)
        fmap = _maybesplat(f)
        if nprocs() == 1 # no workers added
            return _mapreduce(fmap, op, iterable_on_proc...; reducekw...)
        else # one worker or single-valued iterator
            return @fetchfrom workers(pool)[1] _mapreduce(fmap, op, iterable_on_proc...; reducekw...)
        end
    end

    tree_branches = createbranchchannels(pool, N)
    pmapreduceworkers(f, op, tree_branches, iterators; reducekw...)
end

function pmapreduce(f, op, iterators...; reducekw...)
    N = length(zip(iterators...))
    pool = maybetrimmedworkerpool(workers(), N)
    pmapreduce(f, op, pool, iterators...; reducekw...)
end

"""
    pmapreduce_productsplit(f, op, [pool::AbstractWorkerPool], iterators...; reducekw...)

Evaluate a parallel mapreduce over the outer product of elements from `iterators`.
The product of `iterators` is split over the workers available, and each worker is assigned a section
of the product. The function `f` should accept a single argument that is a collection of `Tuple`s.

The keyword arguments `reducekw` are passed on to the reduction.

See also: [`pmapreduce`](@ref)
"""
pmapreduce_productsplit(f, op, pool::AbstractWorkerPool, iterators...; reducekw...) =
    pmapreduce(NoSplat(f), op, pool, Hold(product(iterators...)); reducekw...)

function pmapreduce_productsplit(f, op, iterators...; reducekw...)
    N = length(product(iterators...))
    pool = maybetrimmedworkerpool(workers(), N)
    pmapreduce_productsplit(f, op, pool, iterators...; reducekw...)
end

"""
    pmapbatch(f, [pool::AbstractWorkerPool], iterators...)

Carry out a `pmap` with the `iterators` divided evenly among the available workers.

See also: [`pmapreduce`](@ref)
"""
function pmapbatch(f, pool::AbstractWorkerPool, iterators...)
    pmapreduce((x...) -> [f(x...)], vcat, pool, iterators...)
end

function pmapbatch(f, iterators...)
    N = length(zip(iterators...))
    pool = maybetrimmedworkerpool(workers(), N)
    pmapbatch(f, pool, iterators...)
end

"""
    pmapbatch_productsplit(f, [pool::AbstractWorkerPool], iterators...)

Carry out a `pmap` with the outer product of `iterators` divided evenly among the available workers.
The function `f` must accept a collection of `Tuple`s.

See also: [`pmapbatch`](@ref), [`pmapreduce_productsplit`](@ref)
"""
function pmapbatch_productsplit(f, pool::AbstractWorkerPool, iterators...)
    pmapreduce_productsplit(x -> [f(x)], vcat, pool, iterators...)
end

function pmapbatch_productsplit(f, iterators...)
    N = length(product(iterators...))
    pool = maybetrimmedworkerpool(workers(), N)
    pmapbatch_productsplit(f, pool, iterators...)
end


text:
using Statistics, LinearAlgebra
burninfrac = 0.9
temperaturenum = 6
x_ft = transD_GP.assembleTat1(opt, :x_ftrain, burninfrac=burninfrac, temperaturenum=temperaturenum)
n = transD_GP.assembleTat1(opt, :nodes, burninfrac=burninfrac, temperaturenum=temperaturenum)
nmodels = length(n)
i = nmodels
ftest, = transD_GP.GP.GPfit(opt.K, x_ft[i][2,1:n[i]]', x_ft[i][1,1:n[i]]',
    opt.xall, opt.λ², opt.δ, nogetvars=true, demean=opt.demean, p=2,
    my=mean(opt.fbounds, dims=2))
transD_GP.plot_posterior(csem, opt, burninfrac=burninfrac,
            pdfnormalize=true, cmappdf="inferno_r", temperaturenum=temperaturenum)
##
ax[2].plot(ftest,opt.xall[:])

ext: jl

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
export Conversion

abstract type Conversion{T}<:Operator{T} end

struct ConcreteConversion{S<:Space,V<:Space,T} <: Conversion{T}
    domainspace::S
    rangespace::V
end

ConcreteConversion(a::Space,b::Space)=ConcreteConversion{typeof(a),typeof(b),
        promote_type(rangetype(a),rangetype(b))}(a,b)

function convert(::Type{Operator{T}},C::ConcreteConversion{S,V}) where {T,S,V}
    if T==eltype(C)
        C
    else
        ConcreteConversion{S,V,T}(C.domainspace,C.rangespace)
    end
end

domainspace(C::ConcreteConversion)=C.domainspace
rangespace(C::ConcreteConversion)=C.rangespace

function defaultConversion(a::Space,b::Space)
    if a==b
        Conversion(a)
    elseif conversion_type(a,b)==NoSpace()
        sp=canonicalspace(a)
        if typeof(sp) == typeof(a)
            error("Implement Conversion from " * string(typeof(sp)) * " to " * string(typeof(b)))
        elseif typeof(sp) == typeof(b)
            error("Implement Conversion from " * string(typeof(a)) * " to " * string(typeof(sp)))
        else
            Conversion(a,sp,b)
        end
    else
        error("Implement Conversion from " * string(typeof(a)) * " to " * string(typeof(b)))
    end
end

Conversion(a::Space,b::Space) = defaultConversion(a,b)
Conversion(a::Space) = ConversionWrapper(Operator(I,a))
Conversion() = ConversionWrapper(Operator(I,UnsetSpace()))

## Wrapper
# this allows for a Derivative implementation to return another operator, use a SpaceOperator containing
# the domain and range space
# but continue to know its a derivative

struct ConversionWrapper{S<:Operator,T} <: Conversion{T}
    op::S
end

@wrapper ConversionWrapper

ConversionWrapper(::Type{T},op) where {T} = ConversionWrapper{typeof(op),T}(op)
ConversionWrapper(B::Operator) =
    ConversionWrapper{typeof(B),eltype(B)}(B)
Conversion(A::Space,B::Space,C::Space) =
    ConversionWrapper(Conversion(B,C)*Conversion(A,B))
Conversion(A::Space,B::Space,C::Space,D::Space...) =
    ConversionWrapper(Conversion(C,D...)*Conversion(B,C)*Conversion(A,B))

==(A::ConversionWrapper,B::ConversionWrapper) = A.op==B.op

function convert(::Type{Operator{T}},D::ConversionWrapper) where T
    if T==eltype(D)
        D
    else
        BO=convert(Operator{T},D.op)
        ConversionWrapper{typeof(BO),T}(BO)
    end
end

#promotedomainspace(P::Conversion,sp::Space)=ConversionWrapper(Operator(I, sp))


text:
using LinearAlgebra, DiffEqOperators, Random, Test, BandedMatrices, SparseArrays

# 4th derivative, 4th order
function fourth_deriv_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:8] = [3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3]
    A[2,1:8] = [2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6]

    A[N-1,N-5:end] = reverse([2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6], dims=2)
    A[N,N-5:end] = reverse([3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3], dims=2)

    for i in 3:N-2
        A[i,i-2:i+4] = [-1/6 2.0 -13/2 28/3 -13/2 2.0 -1/6]
    end
    return A
end

function second_deriv_fourth_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:6] = [5/6 -15/12 -1/3 7/6 -6/12 5/60]
    A[N,N-3:end] = [1/12 -6/12 14/12 -4/12 -15/12 10/12]
    for i in 2:N-1
        A[i,i-1:i+3] = [-1/12 4/3 -5/2 4/3 -1/12]
    end
    return A
end

function second_derivative_stencil(N)
  A = zeros(N,N+2)
  for i in 1:N, j in 1:N+2
      (j-i==0 || j-i==2) && (A[i,j]=1)
      j-i==1 && (A[i,j]=-2)
  end
  A
end

@testset "Test Constructor, Multiplication, and Concretization" begin
    # Generate random parameters
    al = rand()
    bl = rand()
    cl = rand()

    ar = rand()
    br = rand()
    cr = rand()
    dx = rand()

    Q = RobinBC((al, bl, cl), (ar, br, cr), dx)
    N = 20
    L = CenteredDifference(4, 4, dx, N)
    L2 = CenteredDifference(2,4, dx, N)
    L1 = UpwindDifference(1, 1, dx, N, 1.)

    function coeff_func(du,u,p,t)
        du .= u
    end

    cL = coeff_func*L
    coeffs = rand(N)
    DiffEqOperators.update_coefficients!(cL,coeffs,nothing,0.0)
    cA = cL*Q
    DiffEqOperators.update_coefficients!(cA,coeffs,nothing,0.0)
    @test cL.coefficients == coeffs == cA.L.coefficients

    cL1 = coeff_func*L1
    coeffs = rand(N)
    DiffEqOperators.update_coefficients!(cL1,coeffs,nothing,0.0)
    cA1 = cL1*Q
    DiffEqOperators.update_coefficients!(cA1,coeffs,nothing,0.0)
    @test cL1.coefficients == coeffs == cA1.L.coefficients

    A = L*Q
    c2A = coeff_func*A
    DiffEqOperators.update_coefficients!(c2A,coeffs,nothing,0.0)
    @test coeffs == c2A.L.coefficients

    # Test GhostDerivativeOperator constructor by *
    u = rand(N)
    A = L*Q
    A1 = L1*Q

    # Test for consistency of GhostDerivativeOperator*u with L*(Q*u)
    @test A == L * Q
    @test A*u ≈ L*(Q*u)
    @test A1 == L1 * Q
    @test A1*u ≈ L1*(Q*u)

    # Test for consistency of c*GhostDerivativeOperator*u with alternative methods
    c = 2.1
    @test c * A == (c * L) * Q == c * (L * Q)
    @test c * A * u ≈ (c * L) * (Q * u) ≈ c * (L * Q) * u
    @test c * A1 == (c * L1) * Q == c * (L1 * Q)
    @test c * A1 * u ≈ (c * L1) * (Q * u) ≈ c * (L1 * Q) * u

    # check A + B, where A and B are GhostDerivativeOperators
    B = c * A
    B1 = c * A1
    @test (A + B) * u == (A + c * A) * u == (B + A) * u == (c * A + A) * u
    @test (A + B) * u == A * u + B * u
    @test (A1 + B1) * u == (A1 + c * A1) * u == (B1 + A1) * u == (c * A1 + A1) * u
    @test (A1 + B1) * u == A1 * u + B1 * u

    # Check (L + L) * Q works
    LLQ = (L + L) * Q
    LLQ1 = (L1 + L1) * Q
    @test LLQ * u == A * u + A * u == (A + A) * u
    @test LLQ1 * u == A1 * u + A1 * u == (A1 + A1) * u

    # Test for consistency of c*GhostDerivativeOperator*u when c is a vector
    c = rand(N)
    L = CenteredDifference(4, 4, dx, N)
    L1 = UpwindDifference(1, 1, 1., N, 1.)
    A1 = L1 * Q
    cA = c * A
    cL = c * L
    cLQ = (c * L) * Q
    cA1 = c * A1
    cL1 = c * L1
    cLQ1 = (c * L1) * Q
    @test cA.L == cL
    @test cA.L.coefficients == cL.coefficients == cLQ.L.coefficients
    @test c * A == (c * L) * Q == c * (L * Q)
    @test c * A * u ≈ (c * L) * (Q * u) ≈ c * (L * Q) * u
    @test cA1.L.coefficients == cL1.coefficients == cLQ1.L.coefficients
    @test c * A1 == (c * L1) * Q == c * (L1 * Q)
    @test c * A1 * u ≈ (c * L1) * (Q * u) ≈ c * (L1 * Q) * u

    # check A + B, where A and B are GhostDerivativeOperators
    B = c * A
    B1 = c * A1
    # @test (A + B) == A + c * A == B + A == c * A + A # uncomment if we want to implement and test for equality of linear combinations of operators
    @test (A + B) * u == (A + c * A) * u == (B + A) * u == (c * A + A) * u
    @test (A + B) * u == A * u + B * u
    # @test (A1 + B1) == A + c * A1 == B1 + A1 == c * A1 + A1
    @test (A1 + B1) * u == (A1 + c * A1) * u == (B1 + A1) * u == (c * A1 + A1) * u
    @test (A1 + B1) * u == A1 * u + B1 * u

    # Test for consistency of GhostDerivativeOperator*M with L*(Q*M)
    M = rand(N,10)
    LQM = zeros(N,10)
    for i in 1:10
        mul!(view(LQM,:,i), L, Q*M[:,i])
    end
    ghost_LQM = A*M
    @test ghost_LQM ≈ LQM

    u = rand(N + 2)
    @test (L + L2) * u ≈ convert(AbstractMatrix,L + L2) * u ≈ (BandedMatrix(L) + BandedMatrix(L2)) * u

    # Test concretization
    @test Array(A)[1] ≈ (Array(L)*Array(Q,N)[1], Array(L)*Array(Q,N)[2])[1]
    @test Array(A)[2] ≈ (Array(L)*Array(Q,N)[1], Array(L)*Array(Q,N)[2])[2]
    @test SparseMatrixCSC(A)[1] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[1]
    @test SparseMatrixCSC(A)[2] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[2]
    @test sparse(A)[1] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[1]
    @test sparse(A)[2] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[2]
    # BandedMatrix not implemented for boundary operator
    @test_broken BandedMatrix(A)[1] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[1]
    @test_broken BandedMatrix(A)[2] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[2]

    # Test that concretization works with multiplication
    u = rand(N)
    @test Array(A)[1]*u + Array(A)[2] ≈ L*(Q*u) ≈ A*u
    @test sparse(A)[1]*u + sparse(A)[2] ≈ L*(Q*u) ≈ A*u
end

@testset "Test Constructor, Multiplication, and Concretization (Non-uniform grid)" begin
    # Generate random parameters
    al = rand()
    bl = rand()
    cl = rand()

    ar = rand()
    br = rand()
    cr = rand()
    N = 20
    dx = rand(N + 2)

    Q = RobinBC((al, bl, cl), (ar, br, cr), dx)
    L = UpwindDifference(1, 1, dx, N + 1, 1.)
    L2 = CenteredDifference(2, 2, dx, N)
    L4 = CenteredDifference(4, 4, dx, N - 2)

    function coeff_func(du,u,p,t)
        du .= u
    end

    cL = coeff_func*L
    coeffs = rand(N + 1)
    DiffEqOperators.update_coefficients!(cL,coeffs,nothing,0.0)
    cA = cL*Q
    DiffEqOperators.update_coefficients!(cA,coeffs,nothing,0.0)
    @test cL.coefficients == coeffs == cA.L.coefficients

    cL2 = coeff_func*L2
    coeffs = rand(N)
    DiffEqOperators.update_coefficients!(cL2,coeffs,nothing,0.0)
    cA2 = cL2*Q
    DiffEqOperators.update_coefficients!(cA2,coeffs,nothing,0.0)
    @test cL2.coefficients == coeffs == cA2.L.coefficients

    cL4 = coeff_func*L4
    coeffs = rand(N - 2)
    DiffEqOperators.update_coefficients!(cL4,coeffs,nothing,0.0)
    cA4 = cL4*Q
    DiffEqOperators.update_coefficients!(cA4,coeffs,nothing,0.0)
    @test cL4.coefficients == coeffs == cA4.L.coefficients

    # Test GhostDerivativeOperator constructor by *
    u = rand(N + 1)
    u2 = rand(N)
    u4 = rand(N - 2)
    A = L*Q
    A2 = L2*Q
    A4 = L4*Q

    # Test for consistency of GhostDerivativeOperator*u with L*(Q*u)
    @test A == L * Q
    @test A*u ≈ L*(Q*u)
    @test A2 == L2 * Q
    @test A2*u2 ≈ L2*(Q*u2)
    @test A4 == L4 * Q
    @test A4*u4 ≈ L4*(Q*u4)

    # Test for consistency of c*GhostDerivativeOperator*u with alternative methods
    c = 2.1
    cA = c * A
    cL = c * L
    @test c * A == (c * L) * Q == c * (L * Q)
    @test c * A * u ≈ (c * L) * (Q * u) ≈ c * (L * Q) * u
    @test c * A2 == (c * L2) * Q == c * (L2 * Q)
    @test c * A2 * u2 ≈ (c * L2) * (Q * u2) ≈ c * (L2 * Q) * u2
    @test c * A4 == (c * L4) * Q == c * (L4 * Q)
    @test c * A4 * u4 ≈ (c * L4) * (Q * u4) ≈ c * (L4 * Q) * u4

    # check A + B, where A and B are GhostDerivativeOperators
    B = c * A
    B2 = c * A2
    B4 = c * A4
    @test (A + B) * u == (A + c * A) * u == (B + A) * u == (c * A + A) * u
    @test (A + B) * u == A * u + B * u
    @test (A2 + B2) * u2 == (A2 + c * A2) * u2 == (B2 + A2) * u2 == (c * A2 + A2) * u2
    @test (A2 + B2) * u2 == A2 * u2 + B2 * u2
    @test (A4 + B4) * u4 == (A4 + c * A4) * u4 == (B4 + A4) * u4 == (c * A4 + A4) * u4
    @test (A4 + B4) * u4 == A4 * u4 + B4 * u4

    # Check (L + L) * Q works
    LLQ = (L + L) * Q
    LLQ2 = (L2 + L2) * Q
    LLQ4 = (L4 + L4) * Q
    @test LLQ * u == A * u + A * u == (A + A) * u
    @test LLQ2 * u2 == A2 * u2 + A2 * u2 == (A2 + A2) * u2 # this comparison fails, even though the operators seem alike
    @test LLQ4 * u4 == A4 * u4 + A4 * u4 == (A4 + A4) * u4

    # Test for consistency of c*GhostDerivativeOperator*u when c is a vector
    c = rand(N + 1)
    L1 = L
    A1 = L1 * Q
    cA = c * A
    cL = c * L
    cLQ = (c * L) * Q
    cA1 = c * A1
    cL1 = c * L1
    cLQ1 = (c * L1) * Q
    @test cA.L == cL
    @test cA.L.coefficients == cL.coefficients == cLQ.L.coefficients
    @test c * A == (c * L) * Q == c * (L * Q)
    @test c * A * u ≈ (c * L) * (Q * u) ≈ c * (L * Q) * u
    @test cA1.L.coefficients == cL1.coefficients == cLQ1.L.coefficients
    @test c * A1 == (c * L1) * Q == c * (L1 * Q)
    @test c * A1 * u ≈ (c * L1) * (Q * u) ≈ c * (L1 * Q) * u

    # check A + B, where A and B are GhostDerivativeOperators
    B = c * A
    # @test (A + B) == A + c * A == B + A == c * A + A # uncomment if implement equality of combinations of linear operators
    @test (A + B) * u == (A + c * A) * u == (B + A) * u == (c * A + A) * u
    @test (A + B) * u == A * u + B * u

    # Test for consistency of GhostDerivativeOperator*M with L*(Q*M)
    M = rand(N+1,10)
    LQM = zeros(N+1,10)
    for i in 1:10
        mul!(view(LQM,:,i), L, Q*M[:,i])
    end
    ghost_LQM = A*M
    @test ghost_LQM ≈ LQM

    # Test concretization for UpwindDifference
    @test Array(A)[1] ≈ (Array(L)*Array(Q,N + 1)[1], Array(L)*Array(Q,N + 1)[2])[1]
    @test Array(A)[2] ≈ (Array(L)*Array(Q,N + 1)[1], Array(L)*Array(Q,N + 1)[2])[2]
    @test SparseMatrixCSC(A)[1] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N + 1)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N + 1)[2])[1]
    @test SparseMatrixCSC(A)[2] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N + 1)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N + 1)[2])[2]
    @test sparse(A)[1] ≈ (sparse(L)*sparse(Q,N + 1)[1], sparse(L)*sparse(Q,N + 1)[2])[1]
    @test sparse(A)[2] ≈ (sparse(L)*sparse(Q,N + 1)[1], sparse(L)*sparse(Q,N + 1)[2])[2]
    # BandedMatrix not implemeted for boundary operator
    @test_broken BandedMatrix(A)[1] ≈ (BandedMatrix(L)*BandedMatrix(Q,N + 1)[1], BandedMatrix(L)*BandedMatrix(Q,N + 1)[2])[1]
    @test_broken BandedMatrix(A)[2] ≈ (BandedMatrix(L)*BandedMatrix(Q,N + 1)[1], BandedMatrix(L)*BandedMatrix(Q,N + 1)[2])[2]

    # Test concretization for CenteredDifference
    A2 = L2 * Q
    @test Array(A2)[1] ≈ (Array(L2)*Array(Q,N)[1], Array(L2)*Array(Q,N)[2])[1]
    @test Array(A2)[2] ≈ (Array(L2)*Array(Q,N)[1], Array(L2)*Array(Q,N)[2])[2]
    @test SparseMatrixCSC(A2)[1] ≈ (SparseMatrixCSC(L2)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L2)*SparseMatrixCSC(Q,N)[2])[1]
    @test SparseMatrixCSC(A2)[2] ≈ (SparseMatrixCSC(L2)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L2)*SparseMatrixCSC(Q,N)[2])[2]
    @test sparse(A2)[1] ≈ (sparse(L2)*sparse(Q,N)[1], sparse(L2)*sparse(Q,N)[2])[1]
    @test sparse(A2)[2] ≈ (sparse(L2)*sparse(Q,N)[1], sparse(L2)*sparse(Q,N)[2])[2]
    # BandedMatrix not implemeted for boundary operator
    @test_broken BandedMatrix(A2)[1] ≈ (BandedMatrix(L2)*BandedMatrix(Q,N)[1], BandedMatrix(L2)*BandedMatrix(Q,N)[2])[1]
    @test_broken BandedMatrix(A2)[2] ≈ (BandedMatrix(L2)*BandedMatrix(Q,N)[1], BandedMatrix(L2)*BandedMatrix(Q,N)[2])[2]

    # Test that concretization works with multiplication, UpwindDifference
    u = rand(N + 1)
    @test Array(A)[1]*u + Array(A)[2] ≈ L*(Q*u) ≈ A*u
    @test sparse(A)[1]*u + sparse(A)[2] ≈ L*(Q*u) ≈ A*u

    # Test that concretization works with multiplication, CenteredDifference
    u = rand(N)
    @test Array(A2)[1]*u + Array(A2)[2] ≈ L2*(Q*u) ≈ A2*u
    @test sparse(A2)[1]*u + sparse(A2)[2] ≈ L2*(Q*u) ≈ A2*u
end

@testset "Test Left Division L2 (second order)" begin
    # Test \ Homogenous BC
    # f(x) = -x^2 + x
    # f''(x) = -2.0
    # f'(0) = f'(1) = 0
    f(x) = -x^2 + x
    f2(x) = -2.0

    dx = 0.01
    x = dx:dx:1.0-dx
    N = length(x)

    L = CenteredDifference(2, 2, dx, N)
    Q = RobinBC((1.0, 0.0, 0.0), (1.0, 0.0, 0.0), dx)
    A = L*Q

    analytic_L = second_derivative_stencil(N) ./ dx^2
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_AL = analytic_L*analytic_QL

    # No affine component to the this system
    analytic_f = analytic_AL \ f2.(x)
    ghost_f = A \ f2.(x)

    # Check that A\f2.(x) is consistent with analytic_AL \ f2.(x)
    @test analytic_f ≈ ghost_f

    # Additionally test that A\f2.(x) ≈ f.(x)
    @test f.(x) ≈ ghost_f ≈ analytic_f

    # Check ldiv!
    f_temp = zeros(N)
    ldiv!(f_temp, A, f2.(x))
    @test f_temp ≈ ghost_f ≈ analytic_f

    # Check that left division with matrices works
    ghost_fM = A \ [f2.(x) f2.(x)]
    analytic_fM = analytic_AL \ [f2.(x) f2.(x)]
    @test ghost_fM ≈ analytic_fM

    fM_temp = zeros(N,2)
    ldiv!(fM_temp, A, [f2.(x) f2.(x)])
    @test fM_temp ≈ analytic_fM

    # Test \ Inhomogenous BC
    # f(x) = -x^2 + x + 4.0
    # f''(x) = -2.0
    # f'(0) = f'(1) = 0
    f(x) = -x^2 + x + 4.0
    f2(x) = -2.0

    dx = 0.01
    x = dx:dx:1.0-dx
    N = length(x)

    L = CenteredDifference(2, 2, dx, N)
    Q = RobinBC((1.0, 0.0, 4.0), (1.0, 0.0, 4.0), dx)
    A = L*Q

    analytic_L = second_derivative_stencil(N) ./ dx^2
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_Qb = [4.0; zeros(N); 4.0]
    analytic_AL = analytic_L*analytic_QL
    analytic_Ab = analytic_L*analytic_Qb

    analytic_f = analytic_AL \ (f2.(x) - analytic_Ab)
    ghost_f = A \ f2.(x)

    # Check that A\f2.(x) is consistent with analytic_AL \ f2.(x)
    @test analytic_f ≈ ghost_f

    # Additionally test that A\f2.(x) ≈ f.(x)
    @test f.(x) ≈ ghost_f ≈ analytic_f

    # Check ldiv!
    f_temp = zeros(N)
    ldiv!(f_temp, A, f2.(x))
    @test f_temp ≈ ghost_f ≈ analytic_f

    # Check \ for Matrix
    M2 = [f2.(x) 2.0*f2.(x) 10.0*f2.(x)]
    analytic_M = analytic_AL \ (M2 .- analytic_Ab)
    ghost_M = A \ M2
    @test analytic_M ≈ ghost_M

    # Check ldiv! for Matrix
    M_temp = zeros(N,3)
    ldiv!(M_temp, A, M2)
    @test M_temp ≈ analytic_M ≈ ghost_M

    # Additionally test that A\M2 ≈ [f, 2.0(f-4.0)+4.0, 10.0(f-4.0)+4.0]
    M = [f.(x) 2.0*(f.(x) .- 4.0).+4.0 10.0*(f.(x) .- 4.0).+4.0]
    @test M ≈ M_temp ≈ analytic_M ≈ ghost_M
end

@testset "Test Left Division L4 (fourth order)" begin

    # Test \ homogenous and inhomogenous BC
    dx = 0.01
    x = 0.01:dx:0.2
    N = length(x)
    u = sin.(x)

    L = CenteredDifference(4, 4, dx, N)
    Q = RobinBC((1.0, 0.0, sin(0.0)), (1.0, 0.0, sin(0.2+dx)), dx)
    A = L*Q

    analytic_L = fourth_deriv_approx_stencil(N) ./ dx^4
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_AL = analytic_L*analytic_QL
    analytic_Qb = [zeros(N+1); sin(0.2+dx)]
    analytic_Ab = analytic_L*analytic_Qb

    analytic_u = analytic_AL \ (u - analytic_Ab)
    ghost_u = A \ u

    # Check that A\u.(x) is consistent with analytic_AL \ u.(x)
    @test analytic_u ≈ ghost_u

    # Check ldiv!
    u_temp = zeros(N)
    ldiv!(u_temp, A, u)
    @test u_temp ≈ ghost_u ≈ analytic_u

    # Check \ for Matrix
    M2 = [u 2.0*u 10.0*u]
    analytic_M = analytic_AL \ (M2 .- analytic_Ab)
    ghost_M = A \ M2
    @test analytic_M ≈ ghost_M

    # Check ldiv! for Matrix
    M_temp = zeros(N,3)
    ldiv!(M_temp, A, M2)
    @test M_temp ≈ ghost_M ≈ analytic_M
end


text:
#!/usr/bin/env julia -

using Luxor

using Test

function svgstring_test()
    @testset "svg string test" begin
        Drawing(400, 800, :svg)
        origin()
        juliacircles()
        @test finish() == true
        svgsource = svgstring()
        matches = collect(eachmatch(r"<.*?>", svgsource))
        @test first(matches).match == """<?xml version="1.0" encoding="UTF-8"?>"""
        @test last(matches).match == "</svg>"
    end
    println("...finished svgstring_test")
end

function place_svgtest(fnamein, fnameout)
    svgin = readsvg(fnamein)
    @testset "readsvg" begin
        @test svgin.xpos == 400
        @test svgin.ypos == 400
        @test svgin.width == 400
        @test svgin.height == 400
    end
    @testset "placeimage" begin
        Drawing(400, 800, fnameout)
        origin()
        background("green")
        placeimage(svgin, O, centered = true)
        @test finish() == true
        svgin2 = readsvg(fnameout)
        @test svgin2.xpos == 400
        @test svgin2.ypos == 800
        @test svgin2.width == 400
        @test svgin2.height == 800
    end
    println("...finished test: output in $(fnameout)")
end

svgstring_test()
place_svgtest("polysample.svg", "place-svg.svg")


text:
module GeneSets

# Write your package code here.
include("gset_struct.jl")
export GeneSet, AbstractGeneSet, get_msigdb

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Test conversion from osculating to mean elements using J2 osc. orbit
#   propagator.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/orbit/mean_elements/j2osc.jl
# ========================================

# Function rv_to_mean_elements_j2osc
# ----------------------------------

@testset "Function rv_to_mean_elements_j2osc" begin
    # We just need to run the J2 osc. orbit propagator, obtain the osculating
    # elements, convert them to mean elements, and compare with the original
    # ones.
    jd₀ = date_to_jd(2021, 6, 19, 19, 35, 35)
    orb_input = KeplerianElements(
        jd₀,
        7130.982e3,
        0.001111,
        deg2rad(98.405),
        compute_RAAN_lt(jd₀, 22.5),
        deg2rad(90),
        deg2rad(10)
    )

    # Generate the osculating elements.
    orbp   = init_orbit_propagator(Val(:J2osc), orb_input)
    vr, vv = propagate!(orbp, 0:10:12000)
    vjd    = get_epoch(orbp) .+ (0:10:12000) ./ 86400

    # Obtain the mean elements.
    orb, ~ = rv_to_mean_elements_j2osc(
        vjd,
        vr,
        vv;
        mean_elements_epoch = :begin,
        print_debug         = false,
        max_iterations      = 1000,
        atol                = 1e-10,
        rtol                = 1e-10
    )

    # Compare.
    @test orb.t == orb_input.t
    @test orb.a ≈  orb_input.a atol = 1e-4
    @test orb.e ≈  orb_input.e atol = 1e-11
    @test orb.i ≈  orb_input.i atol = 1e-10
    @test orb.Ω ≈  orb_input.Ω atol = 1e-10
    @test orb.ω ≈  orb_input.ω atol = 1e-7
    @test orb.f ≈  orb_input.f atol = 1e-7
end


text:
export Poisson, poisson_aux
export Gauss, gauss_aux
export ScaledGauss, scaled_gauss_aux

"""
    poisson_aux(μ, meas)

Calculates the Poisson loss for `μ` and `meas`.
`μ` can be of larger size than `meas`. In that case
we extract a centered region from `μ` of the same size as `meas`.
"""
function poisson_aux(μ, meas)
    μ = center_extract(μ, size(meas))

    μ[μ .<= 1f-8] .= 1f-8
    return sum(μ .- meas .* log.(μ))
end

 # define custom gradient for speed-up
 # ChainRulesCore offers the possibility to define a backward AD rule
 # which can be used by several different AD systems
function ChainRulesCore.rrule(::typeof(poisson_aux), μ, meas)
    Y = poisson_aux(μ, meas)

    function poisson_aux_pullback(xbar)
        meas_new = copy(μ)
        meas_new = center_set!(meas_new, meas)

        ∇ = xbar .* (one(eltype(μ)) .- meas_new ./ μ)
        ∇[μ .< 1f-8] .= 0
        return zero(eltype(μ)), ∇, zero(eltype(μ))
    end

    return Y, poisson_aux_pullback
end

"""
    Poisson()

Returns a function to calculate Poisson loss
Check the help of `poisson_aux`.
"""
function Poisson()
      return poisson_aux
end

"""
    gauss_aux(μ, meas)

Calculates the Gauss loss for `μ` and `meas`.
`μ` can be of larger size than `meas`. In that case
we extract a centered region from `μ` of the same size as `meas`.
"""
function gauss_aux(μ, meas)
    μ = center_extract(μ, size(meas))
    return sum(abs2.(μ - meas))
end

 # define custom gradient for speed-up
function ChainRulesCore.rrule(::typeof(gauss_aux), μ, meas)
    Y = gauss_aux(μ, meas)
    function gauss_aux_pullback(xbar)
        meas_new = copy(μ)
        meas_new = center_set!(meas_new, meas)
        return zero(eltype(μ)), 2 .* (μ - meas_new), zero(eltype(μ))
    end
    return Y, gauss_aux_pullback
end

"""
    Gauss()

Returns a function to calculate Gauss loss.
Check the help of `gauss_aux`.
"""
function Gauss()
    return gauss_aux
end

"""
    scaled_gauss_aux(μ, meas)
Calculates the scaled Gauss loss for `μ` and `meas`.
`μ` can be of larger size than `meas`. In that case
we extract a centered region from `μ` of the same size as `meas`.
"""
function scaled_gauss_aux(μ, meas)
    μ = center_extract(μ, size(meas))
    μ[μ .<= 1f-8] .= 1f-8
    return sum(1/2f0 .* log.(μ) .+ (meas .- μ).^2 ./ (2 .* μ))
end

 # define custom gradient for speed-up
function ChainRulesCore.rrule(::typeof(scaled_gauss_aux), μ, meas)
    Y = scaled_gauss_aux(μ, meas)
    function scaled_gauss_aux_pullback(xbar)
        meas_new = copy(μ)
        meas_new = center_set!(meas_new, meas)
        ∇ = (μ + μ.^2 - meas_new.^2)./(2 .* μ.^2)
        ∇[μ .<= 1f-8] .= 0
        return zero(eltype(μ)), ∇, zero(eltype(μ))
    end
    return Y, scaled_gauss_aux_pullback
end

"""
    ScaledGauss()

Returns a function to calculate scaled Gauss loss.
Check the help of `scaled_gauss_aux`.
"""
function ScaledGauss()
    return scaled_gauss_aux
end


text:
@reexport module LogicBase

export ¬, ∧, ∨, ⊻, ⇒, ⇔, ⋁, ⋀

"""
    ¬(x::Bool) → Bool

A logical "not". Return `false` if the given argument is `true`, and return
`true` if the given argument is `false`.
"""
¬(x::Bool) = ~x

"""
    ∧(x::Bool, y::Bool) → Bool

A logical "and" that does not short-circuit. Return `true` if both arguments are
`true`, and `false` otherwise.
"""
∧(x::Bool, y::Bool) = x & y

"""
    ∨(x::Bool, y::Bool) → Bool

A logical "or" that does not short-circuit. Return `true` if at least one
argument is `true`, and `false` otherwise.
"""
∨(x::Bool, y::Bool) = x | y

"""
    ⊻(x::Bool, y::Bool) → Bool

A logical "exclusive or". Return `true` if exactly one argument is `true`, and
`false` otherwise.
"""
⊻(x::Bool, y::Bool) = x $ y

"""
    ⇒(x::Bool, y::Bool) → Bool

A logical "implies". Return `true` if the first argument is `false`, or if the
second argument is `true`. Return `false` if the first argument is `true` and
the second argument is `false`.
"""
⇒(x::Bool, y::Bool) = ~x | y

"""
    ⇔(x::Bool, y::Bool) → Bool

A logical "if and only if". Return `true` if the first argument is the same as
the second argument.
"""
⇔(x::Bool, y::Bool) = x ≡ y

"""
    ⋁(xs) → Bool

Return `true` if at least one of the `xs` is true.
"""
⋁ = any

"""
    ⋀(xs) → Bool

Return `true` if all of the `xs` are true.
"""
⋀ = all

end  # module

max_stars_repo_name: TotalVerb/FirstOrderLogic

max_issues_repo_name: TotalVerb/FirstOrderLogic

max_forks_repo_name: TotalVerb/FirstOrderLogic


text:
"""
    AbstractMultirouteFlowAlgorithm

Abstract type that allows users to pass in their preferred algorithm.
"""
abstract type AbstractMultirouteFlowAlgorithm end

"""
    KishimotoAlgorithm

Used to specify the Kishimoto algorithm.
"""
struct KishimotoAlgorithm <: AbstractMultirouteFlowAlgorithm end

"""
    ExtendedMultirouteFlowAlgorithm

Used to specify the Extended Multiroute Flow algorithm.
"""
struct ExtendedMultirouteFlowAlgorithm <: AbstractMultirouteFlowAlgorithm end

# Methods when the number of routes is more than the connectivity
# 1) When using Boykov-Kolmogorov as a flow subroutine
# 2) Other flow algorithm
function empty_flow(
    capacity_matrix::AbstractMatrix{T},     # edge flow capacities
    flow_algorithm::BoykovKolmogorovAlgorithm # keyword argument for algorithm
    ) where T<:Real
    n = size(capacity_matrix, 1)
    return zero(T), zeros(T, n, n), zeros(T, n)
end
# 2) Other flow algorithm
function empty_flow(
    capacity_matrix::AbstractMatrix{T},     # edge flow capacities
    flow_algorithm::AbstractFlowAlgorithm     # keyword argument for algorithm
    ) where T<:Real
    n = size(capacity_matrix, 1)
    return zero(T), zeros(T, n, n)
end

# Method for Kishimoto algorithm
@traitfn function multiroute_flow(
    flow_graph::::lg.IsDirected,                   # the input graph
    source::Integer,                       # the source vertex
    target::Integer,                       # the target vertex
    capacity_matrix::AbstractMatrix,  # edge flow capacities
    flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
    mrf_algorithm::KishimotoAlgorithm,     # keyword argument for algorithm
    routes::Int                            # keyword argument for routes
    )
    return kishimoto(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)
end

## Methods for Extended Multiroute Flow Algorithm
#1 When the breaking points are not already known
@traitfn function multiroute_flow(
    flow_graph::::lg.IsDirected,                       # the input graph
    source::Integer,                                # the source vertex
    target::Integer,                                # the target vertex
    capacity_matrix::AbstractMatrix,           # edge flow capacities
    flow_algorithm::AbstractFlowAlgorithm,          # keyword argument for algorithm
    mrf_algorithm::ExtendedMultirouteFlowAlgorithm, # keyword argument for algorithm
    routes::Real                                       # keyword argument for routes
    )
    return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)
end
#2 When the breaking points are already known
#2-a Output: flow value (paired with the associated restriction)
multiroute_flow(
    breakingpoints::Vector{Tuple{T,T,Int}},       # vector of breaking points
    routes::R                                       # keyword argument for routes
    ) where T<:Real where R<:Real =
    intersection(breakingpoints, routes)

#2-b Output: flow value, flows(, labels)
function multiroute_flow(
    breakingpoints::AbstractVector{Tuple{T1,T1,Int}}, # vector of breaking points
    routes::R,                                # keyword argument for routes
    flow_graph::lg.AbstractGraph,                      # the input graph
    source::Integer,                          # the source vertex
    target::Integer,                          # the target vertex
    capacity_matrix::AbstractMatrix{T2} =   # edge flow capacities
    DefaultCapacity(flow_graph);
    flow_algorithm::AbstractFlowAlgorithm  =  # keyword argument for algorithm
    PushRelabelAlgorithm()
    ) where T2 where T1<:Real where R<:Real
    x, f = intersection(breakingpoints, routes)
    # For other cases, capacities need to be Floats
    if !(T2<:AbstractFloat)
        capacity_matrix = convert(AbstractMatrix{Float64}, capacity_matrix)
    end

    return maximum_flow(flow_graph, source, target, capacity_matrix,
    algorithm = flow_algorithm, restriction = x)
end

### TODO: CLEAN UP THIS FUNCTION AND DOCUMENTATION. THERE SHOULD BE NO NEED TO
### HAVE A TYPE-UNSTABLE FUNCTION HERE. (sbromberger 2017-03-26)
"""
    multiroute_flow(flow_graph, source, target[, DefaultCapacity][, flow_algorithm][, mrf_algorithm][, routes])

The generic multiroute_flow function.

The output will vary depending on the input:

- When the number of `route`s is `0`, return the set of breaking points of
the multiroute flow.
- When the number of `route`s is `1`, return a flow with a set of 1-disjoint paths
(this is the classical max-flow implementation).
- When the input is limited to a set of breaking points and a route value `k`,
return only the k-route flow.
- Otherwise, a tuple with 1) the maximum flow and 2) the flow matrix. When the
max-flow subroutine is the Boykov-Kolmogorov algorithm, the associated mincut is
returned as a third output.

When the input is a network, it requires the following arguments:

- `flow_graph`: the input graph
- `source`: the source vertex
- `target`: the target vertex
- `capacity_matrix`: matrix of edge flow capacities
- `flow_algorithm`: keyword argument for flow algorithm
- `mrf_algorithm`: keyword argument for multiroute flow algorithm
- `routes`: keyword argument for the number of routes

When the input is only the set of (breaking) points and the number of route,
it requires the following arguments:

- `breakingpoints`: vector of breaking points
- `routes`: number of routes

When the input is the set of (breaking) points, the number of routes,
and the network descriptors, it requires the following arguments:

- `breakingpoints`: vector of breaking points
- `routes`: number of routes
- `flow_graph`: the input graph
- `source`: the source vertex
- `target`: the target vertex
- `capacity_matrix`: matrix of edge flow capacities
- `flow_algorithm`: keyword argument for flow algorithm

The function defaults to the Push-relabel (classical flow) and Kishimoto
(multiroute) algorithms. Alternatively, the algorithms to be used can also
be specified through  keyword arguments. A default capacity of `1` is assumed
for each link if no capacity matrix is provided.

The `mrf_algorithm` keyword is inforced to Extended Multiroute Flow
in the following cases:

- The number of routes is non-integer
- The number of routes is 0 or non-specified

### Usage Example :
(please consult the  [`maximum_flow`](@ref) section for options about flow_algorithm
and capacity_matrix)

```julia
julia> flow_graph = lg.DiGraph(8) # Create a flow graph

julia> flow_edges = [
(1, 2, 10), (1, 3, 5),  (1, 4, 15), (2, 3, 4),  (2, 5, 9),
(2, 6, 15), (3, 4, 4),  (3, 6, 8),  (4, 7, 16), (5, 6, 15),
(5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6),  (7, 8, 10)
]

julia> capacity_matrix = zeros(Int, 8, 8) # Create a capacity matrix

julia> for e in flow_edges
    u, v, f = e
    add_edge!(flow_graph, u, v)
    capacity_matrix[u, v] = f
end

julia> f, F = multiroute_flow(flow_graph, 1, 8, capacity_matrix, routes = 2) # Run default multiroute_flow with an integer number of routes = 2

julia> f, F = multiroute_flow(flow_graph, 1, 8, capacity_matrix, routes = 1.5) # Run default multiroute_flow with a noninteger number of routes = 1.5

julia> points = multiroute_flow(flow_graph, 1, 8, capacity_matrix) # Run default multiroute_flow for all the breaking points values

julia> f, F = multiroute_flow(points, 1.5) # Then run multiroute flow algorithm for any positive number of routes

julia> f = multiroute_flow(points, 1.5, valueonly = true)

julia> f, F, labels = multiroute_flow(flow_graph, 1, 8, capacity_matrix, algorithm = BoykovKolmogorovAlgorithm(), routes = 2) # Run multiroute flow algorithm using Boykov-Kolmogorov algorithm as maximum_flow routine

```
"""
function multiroute_flow(
        flow_graph::lg.AbstractGraph,                    # the input graph
        source::Integer,                                 # the source vertex
        target::Integer,                                 # the target vertex
        capacity_matrix::AbstractMatrix{T} =             # edge flow capacities
            DefaultCapacity(flow_graph);
        flow_algorithm::AbstractFlowAlgorithm =          # keyword argument for algorithm
        PushRelabelAlgorithm(),
        mrf_algorithm::AbstractMultirouteFlowAlgorithm = # keyword argument for algorithm
        KishimotoAlgorithm(),
        routes::R = 0                                    # keyword argument for number of routes (0 = all values)
    ) where T where R <: Real

    # a flow with a set of 1-disjoint paths is a classical max-flow
    (routes == 1) &&
    return maximum_flow(flow_graph, source, target, capacity_matrix, flow_algorithm)

    # routes > λ (connectivity) → f = 0
    λ = maximum_flow(flow_graph, source, target, DefaultCapacity(flow_graph),
    algorithm = flow_algorithm)[1]
    (routes > λ) && return empty_flow(capacity_matrix, flow_algorithm)

    # For other cases, capacities need to be Floats
    if !(T<:AbstractFloat)
        capacity_matrix = convert(AbstractMatrix{Float64}, capacity_matrix)
    end

    # Ask for all possible values (breaking points)
    (routes == 0) &&
    return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm)
    # The number of routes is a float → EMRF
    (R <: AbstractFloat) &&
    return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)

    # Other calls
    return multiroute_flow(flow_graph, source, target, capacity_matrix,
    flow_algorithm, mrf_algorithm, routes)
end


text:
# Probability that
dot(x::T, y::T) where {T<:BitArray} = x.&y
dot(x::T, y::T) where {T<:Bool} = x&y

dist(x::T, y::T) where {T<:Bool} = dot(x, y)

# TODO sum needs to customized
# TODO dot needs to be renamed to something else
# TODO dist should be probability that x != y given function outputs

dist(x::T, y::T) where {T<:BitArray} = sum(dot(x, y))/length(x)

bloss(x, y) = dist(x, y)


text:
# 2019.10.06 (dahong): this was an attempt to implement the iterative algorithm
# using iterables to make the code more modular, but this seemed to be less
# memory efficient for some reason, so we are just using loops for now.
# todos: figure out memory inefficiency, make types parametric as appropriate

using IterTools

import Base: iterate, IteratorSize, IsInfinite, SizeUnknown, eltype, tail

struct FilterHaltIterable{T}
    it::T
    tol
end
IteratorSize(::Type{<:FilterHaltIterable}) = SizeUnknown()
eltype(::Type{FilterHaltIterable{I}}) where {I} = eltype(I)

function iterate(fh::FilterHaltIterable{T},state=(false,copy(fh.it.H0),)) where{T}
    halt, Hprev, itstate = state[1], state[2], tail(tail(state))
    halt && return nothing

    itnext = iterate(fh.it,itstate...)
    itnext === nothing && return nothing

    H = itnext[2].H
    haltnext = (normdiff(H,Hprev)/norm(H) <= fh.tol)::Bool

    copyto!(Hprev,H)
    return itnext[1],(haltnext,Hprev,itnext[2])
end

struct CAOLIterable{fmatT}
    x
    H0::fmatT
    R
    λ

    CAOLIterable{fmatT}(x,H0,R,λ) where {fmatT} =
        !(H0'H0 ≈ (1/prod(R))*I) ?
            error("Initial filters not orthonormal.") : new(x,H0,R,λ)
end
CAOLIterable(x,H0::fmatT,R,λ) where {fmatT} = CAOLIterable{fmatT}(x,H0,R,λ)
IteratorSize(::Type{<:CAOLIterable}) = IsInfinite()
eltype(::Type{CAOLIterable{fmatT}}) where {fmatT} = Tuple{fmatT,eltype(fmatT)}

struct CAOLState
    xpad   # padded images
    H      # vectorized form
    h      # natural form view

    # Temporary variables
    zlk
    ΨZ
    ψz
    ψztemp
    HΨZ
    UVt
end
function CAOLState(it::CAOLIterable)   # Form initial state from CAOLIterable
    K = size(it.H0,2)

    # Padded images
    xpad = [padarray(xl,Pad(:circular,ntuple(_->0,ndims(xl)),it.R)) for xl in it.x]

    # Initial filters
    H = copy(it.H0)
    h = [reshape(view(H,:,k),map(n->1:n,it.R)) for k in 1:K]

    # Temporary variables
    zlk = similar(first(it.x),map(n->0:n-1,size(first(it.x))))
    ΨZ = similar(H)
    ψz = [reshape(view(ΨZ,:,k),map(n->1:n,it.R)) for k in 1:K]
    ψztemp = similar(first(ψz))
    HΨZ = similar(H,K,K)
    UVt = HΨZ

    return CAOLState(xpad,H,h,zlk,ΨZ,ψz,ψztemp,HΨZ,UVt)
end

_obj(zlk,λ) = sum(z -> (abs(z) < sqrt(2λ)) ? abs2(z)/2 : λ, zlk)
function iterate(it::CAOLIterable{fmatT},s::CAOLState=CAOLState(it)) where {fmatT}
    # Compute objective and ΨZ
    obj = zero(eltype(it.H0))
    fill!(s.ΨZ,zero(eltype(s.ΨZ)))
    for xpadl in s.xpad, k in 1:length(s.h)
        imfilter!(s.zlk,xpadl,(s.h[k],),NoPad(),Algorithm.FIR())
        obj += _obj(s.zlk,it.λ)
        s.zlk .= hard.(s.zlk,sqrt(2*it.λ))
        imfilter!(s.ψztemp,xpadl,(s.zlk,),NoPad(),Algorithm.FIR())
        s.ψz[k] .+= s.ψztemp
    end

    # Update filter via polar factorization
    mul!(s.HΨZ,it.H0',s.ΨZ)
    F = svd!(s.HΨZ)
    mul!(s.UVt,F.U,F.Vt)
    mul!(s.H,it.H0,s.UVt)

    return (s.H,obj)::eltype(it),s
end

function CAOL(x,h0::Vector,λ,niters,tol)
    R, K = size(h0[1]), length(h0)

    H0 = similar(h0[1],prod(R),K) # vectorized form
    for k in 1:K
        H0[:,k] = vec(h0[k])
    end

    return CAOL(x,H0,R,λ,niters,tol)
end
function CAOL(x,H0,R,λ,niters,tol)
    outs = collect(imap(deepcopy,Iterators.take(FilterHaltIterable(CAOLIterable(x,H0,R,λ),tol),niters)))

    H_trace = getindex.(outs,1)
    H_convergence = [normdiff(H_trace[t],t == 1 ? H0 : H_trace[t-1])/norm(H_trace[t]) for t in 1:length(H_trace)]
    obj_fnc_vals = getindex.(outs,2)
    niter = length(outs)

    h = [reshape(view(H_trace[end],:,k),map(n->1:n,R)) for k in 1:size(H_trace[end],2)]

    return (h,niter,obj_fnc_vals,H_trace,H_convergence)
end

function CAOLfast(x,h0::Vector,λ,niters,tol)
    R, K = size(h0[1]), length(h0)

    H0 = similar(h0[1],prod(R),K) # vectorized form
    for k in 1:K
        H0[:,k] = vec(h0[k])
    end

    return CAOLfast(x,H0,R,λ,niters,tol)
end

function CAOLfast(x,H0,R,λ,niters,tol)
    H = H0
    for outs in Iterators.take(FilterHaltIterable(CAOLIterable(x,H0,R,λ),tol),niters)
        H = outs[1]
    end

    return H
end


text: abstract type FCS end


text:
import Base.ndims, Base.inv
using AffineTransforms

abstract CoordSystem

immutable AbsoluteCoordSystem <: CoordSystem
    name :: Symbol
end

type RelativeCoordSystem <: CoordSystem
    name :: Symbol
end

absolute_coord = Set([:MNI, :Talairach])

coordsys(name :: Symbol, absolute :: Bool) = absolute ? AbsoluteCoordSystem(name) : RelativeCoordSystem(name)
coordsys(name :: Symbol) = coordsys(name, name in absolute_coord)

mni() = coordsys(:MNI, true)
talairach() = coordsys(:Talairach, true)
voxel() = coordsys(:voxel, false)
scanner() = coordsys(:scanner, false)

immutable Coordmap
    domain :: CoordSystem
    range :: CoordSystem
    transformation
end

function *(coord1 :: Coordmap, coord2 :: Coordmap)
    if coord1.domain != coord2.range
        error("Coordinate maps can not be composed because of inconsistent input, output coordinate systems")
    end
    Coordmap(coord2.domain, coord1.range, coord1.transformation * coord2.transformation)
end

inv(coord :: Coordmap) = Coordmap(coord.range, coord.domain, inv(coord.transformation))

ndims(coord :: Coordmap) = ndims(transformation)

typealias CoordmapSet Set{Coordmap}

function find_transformation(coord_set :: CoordmapSet, initial_coordsys :: CoordSystem, target_coordsys :: CoordSystem, ndims :: Integer)
    found_coordsys = Dict{CoordSystem, Vector{Coordmap}}([initial_coordsys], {[Coordmap(initial_coordsys, initial_coordsys, tformeye(ndims))]})
    previous_run = 0
    while found_coordsys.count > previous_run
        previous_run = found_coordsys.count
        for coordmap in coord_set
            if (coordmap.domain in keys(found_coordsys)) && ~(coordmap.range in keys(found_coordsys))
                found_coordsys[coordmap.range] = [[coordmap], found_coordsys[coordmap.domain]]
            elseif coordmap.range in keys(found_coordsys) && ~(coordmap.domain in keys(found_coordsys))
                found_coordsys[coordmap.domain] = [[inv(coordmap)], found_coordsys[coordmap.range]]
            end
            if target_coordsys in keys(found_coordsys)
                return prod(found_coordsys[target_coordsys]).transformation
            end
        end
        if target_coordsys in keys(found_coordsys)
            return prod(found_coordsys[target_coordsys]).transformation
        end
    end
    error("No transformation to target coordinate system has been found")
end


text:
"""
運動学
"""
module Kinematics

using LinearAlgebra

export q_neutral
export q_max
export q_min
export calc_all
export cpoints_local

"""DHパラメータ（修正？）"""
mutable struct DHparam{T}
    α::T
    a::T
    d::T
    θ::T
end

# パラメータ
const L = 278e-3
const h = 64e-3
const H = 1104e-3
const L0 = 270.35e-3
const L1 = 69e-3
const L2 = 364.35e-3
const L3 = 69e-3
const L4 = 374.29e-3
const L5 = 10e-3
const L6 = 368.3e-3

const q_neutral = [0.0, -31.0, 0.0, 43.0, 0.0, 72.0, 0.0] * pi/180  # ニュートラルの姿勢
const q_max = [51.0, 60.0, 173.0, 150.0, 175.0, 120.0, 175.0] * pi/180
const q_min = [-141.0, -123.0, -173.0, -3.0, -175.0, -90.0, -175.0] * pi/180

const DHparams_neutral = [
    DHparam(0.0, 0.0, 0.0, q_neutral[1])
    DHparam(-pi/2, L1, 0.0, q_neutral[2]+pi/2)
    DHparam(pi/2, 0.0, L2, q_neutral[3])
    DHparam(-pi/2, L3, 0.0, q_neutral[4])
    DHparam(pi/2, 0.0, L4, q_neutral[5])
    DHparam(-pi/2, L5, 0.0, q_neutral[6])
    DHparam(pi/2, 0.0, 0.0, q_neutral[7])
]

# 制御点
const cpoints_local = (
    (
        [0.0; L1/2; -L0/2; 1.0],
        # [0.0; -L1/2; -L0/2; 1.0],
        # [L1/2; 0.0; -L0/2; 1.0],
        # [-L1/2; 0.0; -L0/2; 1.0],
    ),  # 1
    (
        # [0.0; 0.0; L3/2; 1.0],
        [0.0; 0.0; -L3/2; 1.0],
    ),  # 2
    (
        [0.0; L3/2; -L2*2/3; 1.0],
        [0.0; -L3/2; -L2*2/3; 1.0],
        [L3/2; 0.0; -L2*2/3; 1.0],
        [-L3/2; 0.0; -L2*2/3; 1.0],
        [0.0; L3/2; -L2*1/3; 1.0],
        [0.0; -L3/2; -L2*1/3; 1.0],
        [L3/2; 0.0; -L2*1/3; 1.0],
        [-L3/2; 0.0; -L2*1/3; 1.0],
    ),  # 3
    (
        [0.0; 0.0; L3/2; 1.0],
        [0.0; 0.0; -L3/2; 1.0],
    ),  # 4
    (
        [0.0; L3/2; -L4/2; 1.0],
        [0.0; -L3/2; -L4/2; 1.0],
        [L3/2; 0.0; -L4/2; 1.0],
        [-L3/2; 0.0; -L4/2; 1.0],
    ),  # 5
    (
        [0.0; 0.0; L3/2; 1.0],
        [0.0; 0.0; -L3/2; 1.0]
    ),  # 6
    (
        [0.0; L3/2; L6*1/3; 1.0],
        [0.0; -L3/2; L6*1/3; 1.0],
        [L3/2; 0.0; L6*1/3; 1.0],
        [-L3/2; 0.0; L6*1/3; 1.0],
        [0.0; L3/2; L6*2/3; 1.0],
        [0.0; -L3/2; L6*2/3; 1.0],
        [L3/2; 0.0; L6*2/3; 1.0],
        [-L3/2; 0.0; L6*2/3; 1.0],
    ),  # 7
    (
        [0.0; 0.0; L6/3; 1.0],
    ),  # 8
)

# 右手
const HTM_BR_Wo = [
    -sqrt(2)/2 sqrt(2)/2 0.0 -L
    -sqrt(2)/2 -sqrt(2)/2 0.0 -h
    0.0 0.0 1.0 H
    0.0 0.0 0.0 1.0
]

const HTM_0_BR = [
    1.0 0.0 0.0 0.0
    0.0 1.0 0.0 0.0
    0.0 0.0 1.0 L0
    0.0 0.0 0.0 1.0
]

# 左手
const HTM_BL_Wo = [
    sqrt(2)/2 sqrt(2)/2 0.0 L
    -sqrt(2)/2 sqrt(2)/2 0.0 -h
    0.0 0.0 1.0 H
    0.0 0.0 0.0 1.0
]

const HTM_0_BL = HTM_0_BR

const HTM_GR_7 = [
    1.0 0.0 0.0 0.0
    0.0 1.0 0.0 0.0
    0.0 0.0 1.0 L6
    0.0 0.0 0.0 1.0
]

const HTM_A = [
    0.0 -1.0 0.0 0.0
    1.0 0.0 0.0 0.0
    0.0 0.0 0.0 0.0
    0.0 0.0 0.0 0.0
]  # 偏微分演算行列

"""DHparamを更新"""
function update_DHparams(q, DHparams=DHparams_neutral)
    #q[2] = q[2] + pi/2
    for i in 1:7
        if i == 2
            DHparams[i].θ = q[i] + pi/2
        else
            DHparams[i].θ = q[i]
        end
    end
    DHparams
end

"""
    HTM(p::DHparam)
同時変換行列(i-1)T(i)
"""
function HTM(p::DHparam{T}) where T
    [
        cos(p.θ) -sin(p.θ) 0.0 p.a
        sin(p.θ)*cos(p.α) cos(p.θ)*cos(p.α) -sin(p.α) -p.d*sin(p.α)
        sin(p.θ)*sin(p.α) cos(p.θ)*sin(p.α) cos(p.α) p.d*cos(p.α)
        0.0 0.0 0.0 1.0
    ]
end

"""同時変換行列を全部計算"""
function calc_HTMs_local_and_global(DHparams::Vector{DHparam{T}}) where T
    HTMs_local = Vector{Matrix{T}}(undef, 10)
    HTMs_global = Vector{Matrix{T}}(undef, 10)
    #_H = Matrix{T}(undef, 4, 4)
    for i in 1:10
        if i == 1
            HTMs_local[i] = HTM_BL_Wo
            HTMs_global[i] = HTM_BL_Wo
        elseif i == 2
            HTMs_local[i] = HTM_0_BL
            HTMs_global[i] =  HTMs_global[i-1] * HTM_0_BL
        elseif i == 10
            HTMs_local[i] = HTM_GR_7
            HTMs_global[i] = HTMs_global[i-1] * HTM_GR_7
        else
            HTMs_local[i] = HTM(DHparams[i-2])
            HTMs_global[i] = HTMs_global[i-1] * HTMs_local[i]
        end
    end
    HTMs_local, HTMs_global
end

"""微分？同時変換行列を全部計算"""
function calc_dHTMs(
    HTMs_local::Vector{Matrix{T}},
    HTMs_global::Vector{Matrix{T}},
) where T
    _dTdq_all = Vector{Vector{Matrix{T}}}(undef, 10)
    for i in 1:7
        _dTdq  = Vector{Matrix{T}}(undef, 10)
        for j in 1:8
            if j < i
                _dTdq[j] = zeros(T, 4, 4)
            elseif j == i
                _dTdq[j] = HTMs_global[j+2] * HTM_A
            else
                _dTdq[j] = _dTdq[j-1] * HTMs_local[j+2]
            end
        end
        _dTdq_all[i] = _dTdq
    end

    Jax_all = Vector{Matrix{T}}(undef, 8)
    Jay_all = Vector{Matrix{T}}(undef, 8)
    Jaz_all = Vector{Matrix{T}}(undef, 8)
    Jo_all = Vector{Matrix{T}}(undef, 8)

    for i in 1:8
        _Jax = Matrix{T}(undef, 4, 7)
        _Jay = Matrix{T}(undef, 4, 7)
        _Jaz = Matrix{T}(undef, 4, 7)
        _Jo = Matrix{T}(undef, 4, 7)
        for j in 1:7
            _Jax[:, j] = _dTdq_all[j][i][:, 1]
            _Jay[:, j] = _dTdq_all[j][i][:, 2]
            _Jaz[:, j] = _dTdq_all[j][i][:, 3]
            _Jo[:, j] = _dTdq_all[j][i][:, 4]
        end
        Jax_all[i] = _Jax
        Jay_all[i] = _Jay
        Jaz_all[i] = _Jaz
        Jo_all[i] = _Jo
    end
    Jax_all, Jay_all, Jaz_all, Jo_all
end

"""Joのヤコビ行列"""
function _calc_Jo_global(Jax, Jay, Jaz, Jo, r_bar)
    z_bar = Jax .* r_bar[1,1] .+ Jay .* r_bar[2,1] .+ Jaz .* r_bar[3,1] .+ Jo
    return z_bar[1:3, :]
end

"""各制御点のヤコビ行列を計算"""
function calc_jacobians(
    Jax_all::Vector{Matrix{T}},
    Jay_all::Vector{Matrix{T}},
    Jaz_all::Vector{Matrix{T}},
    Jo_all::Vector{Matrix{T}},
) where T

    # ジョイント基底位置ベクトルのヤコビ行列を計算
    Jos_joint_all = Vector{Matrix{T}}(undef, 8)
    for i in 1:8
        Jos_joint_all[i] = _calc_Jo_global(
            Jax_all[i], Jay_all[i], Jaz_all[i], Jo_all[i], zeros(T, 4),  # <-正確には[0,0,0,1]
        )
    end

    # 制御点の位置に関するヤコビ行列を計算
    Jos_cpoints_all = Vector{Vector{Matrix{T}}}(undef, 8)
    for i in 1:8
        n = length(cpoints_local[i])
        _J = Vector{Matrix{T}}(undef, n)
        for j in 1:n
            _J[j] = _calc_Jo_global(
                Jax_all[i], Jay_all[i], Jaz_all[i], Jo_all[i], cpoints_local[i][j],
            )
        end
        Jos_cpoints_all[i] = _J
    end

    return Jos_joint_all, Jos_cpoints_all
end

"""制御点の位置と速度を計算

"""
function calc_cpoint_x_and_dx_global(
    HTMs_global::Vector{Matrix{T}},
    Jos_cpoints_all::Vector{Vector{Matrix{T}}},
    dq::Vector{T},
) where T
    cpoints_x_global = Vector{Vector{Vector{T}}}(undef, 8)
    for i in 1:8
        n = length(cpoints_local[i])
        _c = Vector{Vector{T}}(undef, n)
        for j in 1:n
            _c[j] = (HTMs_global[i+2] * cpoints_local[i][j])[1:3]
        end
        cpoints_x_global[i] = _c
    end

    cpoints_dx_global = Vector{Vector{Vector{T}}}(undef, 8)
    for i in 1:8
        n = length(cpoints_local[i])
        _dx = Vector{Vector{T}}(undef, n)
        for j in 1:n
            _dx[j] = Jos_cpoints_all[i][j] * dq
        end
        cpoints_dx_global[i] = _dx
    end

    return cpoints_x_global, cpoints_dx_global
end

"""原点+ジョイント位置を取得"""
function calc_joints_x_and_dx_global(
    HTMs_global::Vector{Matrix{T}},
    Jos_joint_all::Vector{Matrix{T}},
    dq::Vector{T},
) where T
    joints_x_global = Vector{Vector{T}}(undef, 11)
    joints_x_global[1] = zeros(T, 3)
    for i in 2:length(joints_x_global)
        joints_x_global[i] = HTMs_global[i-1][1:3, 4]
    end

    joints_dx_global = Vector{Vector{T}}(undef, 11)
    joints_dx_global[1] = zeros(T, 3)
    joints_dx_global[2] = zeros(T, 3)
    joints_dx_global[3] = zeros(T, 3)
    for i in 4:length(joints_dx_global)
        joints_dx_global[i] = Jos_joint_all[i-3] * dq
    end

    return joints_x_global, joints_dx_global
end

"""全部計算"""
function calc_all(q, dq,)
    DHparams = update_DHparams(q)
    HTMs_local, HTMs_global = calc_HTMs_local_and_global(DHparams)
    Jax_all, Jay_all, Jaz_all, Jo_all = calc_dHTMs(HTMs_local, HTMs_global)
    Jos_joint_all, Jos_cpoint_all = calc_jacobians(Jax_all, Jay_all, Jaz_all, Jo_all)
    cpoints_x_global, cpoints_dx_global = calc_cpoint_x_and_dx_global(
        HTMs_global, Jos_cpoint_all, dq
    )
    joints_x_global, joints_dx_global = calc_joints_x_and_dx_global(
        HTMs_global, Jos_joint_all, dq
    )

    return (
        HTMs_local, HTMs_global,
        Jax_all, Jay_all, Jaz_all, Jo_all,
        Jos_joint_all, Jos_cpoint_all,
        cpoints_x_global, cpoints_dx_global,
        joints_x_global, joints_dx_global,
    )
end

end


text:
module FEMMShellT3FFCompModule

using LinearAlgebra: norm, Transpose, mul!, diag, eigen, I, dot, rank
using Statistics: mean
using FinEtools
import FinEtools.FESetModule: gradN!, nodesperelem, manifdim
using FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.FEMMBaseModule: associategeometry!
import FinEtools.FEMMBaseModule: inspectintegpoints
using FinEtoolsDeforLinear.MatDeforLinearElasticModule: tangentmoduli!, update!, thermalstrain!
using FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, locjac!, add_nnt_ut_only!, add_btsigma!, add_b1tdb2!
using ..FESetShellT3Module: FESetShellT3
using ..TransformerModule: QTEQTransformer, QEQTTransformer, Layup2ElementAngle
using ..CompositeLayupModule: CompositeLayup, thickness, laminate_stiffnesses!, laminate_transverse_stiffness!, plane_stress_T_matrix!, transverse_shear_T_matrix!, laminate_inertia!

const __nn = 3 # number of nodes
const __ndof = 6 # number of degrees of freedom per node

# Formulation for the transverse shear stiffness which averages the
# strain-displacement matrix.
const __TRANSV_SHEAR_FORMULATION_AVERAGE_B = 0
# Formulation for the transverse shear stiffness which averages the
# stiffness matrix.
const __TRANSV_SHEAR_FORMULATION_AVERAGE_K = 1

"""
    FEMMShellT3FFComp{S<:FESetT3, F<:Function} <: AbstractFEMM

Class for the finite element modeling machine of the T3 triangular Flat-Facet
shell with the Discrete Shear Gap technology and a consistent handling of the
normals. This formulation is suitable for modelling of COMPOSITE (layered) materials.

For details for the homogeneous-shell refer to [`FEMMShellT3FF`](@ref).
"""
mutable struct FEMMShellT3FFComp{S<:FESetT3, F<:Function} <: AbstractFEMM
    integdomain::IntegDomain{S, F} # integration domain data
    # Definitions of layups
    layup_groups::Vector{Tuple{CompositeLayup, Vector{FInt}}} # layups: vector of pairs of the composite layup and the group of elements using that layup.
    # Configuration parameters
    transv_shear_formulation::FInt
    drilling_stiffness_scale::Float64
    threshold_angle::Float64
    mult_el_size::Float64
    # Private data
    _associatedgeometry::Bool
    _normals::FFltMat
    _normal_valid::Vector{Bool}
    _layup_group_lookup::FIntVec # look up the layup group for each element
    # The attributes below are buffers used in various operations.
    _loc::FFltMat
    _J0::FFltMat
    _ecoords::FFltMat
    _edisp::FFltVec
    _ecoords_e::FFltMat
    _edisp_e::FFltVec
    _dofnums::FIntMat
    _E_G::FFltMat
    _A_Es::Vector{FFltMat} # transformation nodal-element matrices
    _nvalid::Vector{Bool}
    _T::FFltMat  # element transformation matrix
    _elmat::FFltMat
    _gradN_e::FFltMat
    _Bm::FFltMat
    _Bb::FFltMat
    _Bs::FFltMat
    _DpsBmb::FFltMat
    _DtBs::FFltMat
end

"""
    FEMMShellT3FFComp(integdomain::IntegDomain{S, F}, layup::CompositeLayup) where {S<:FESetT3, F<:Function}

Constructor of the T3FFComp shell FEMM. All elements use a single layup.
"""
function FEMMShellT3FFComp(integdomain::IntegDomain{S,F}, layup::CompositeLayup) where {S<:FESetT3,F<:Function}
    _nnmax = 0
    for j in 1:count(integdomain.fes)
        for k in eachindex(integdomain.fes.conn[j])
            _nnmax = max(_nnmax, integdomain.fes.conn[j][k])
        end
    end
    # When there is only a single layup, all elements belong to a single group
    layup_groups = [(layup, collect(1:count(integdomain.fes)))]
    # Establish the mapping from the elements to the layup group
    _layup_group_lookup = fill(zero(FInt), count(integdomain.fes))
    for j in 1:length(layup_groups)
        eset = layup_groups[j][2]
        for i in eset # Loop over elements in the layup group
            _layup_group_lookup[i] = j
        end
    end
    # Allocate space for the normals
    _normals = fill(0.0, _nnmax, 3)
    _normal_valid = fill(true, _nnmax)
    # Alocate the buffers
    _loc = fill(0.0, 1, 3)
    _J0 = fill(0.0, 3, 2)
    _ecoords = fill(0.0, __nn, 3)
    _edisp = fill(0.0, __nn * __ndof)
    _ecoords_e = fill(0.0, __nn, 2)
    _edisp_e = fill(0.0, __nn * __ndof)
    _dofnums = zeros(FInt, 1, __nn * __ndof)
    _E_G = fill(0.0, 3, 3)
    _A_Es = [fill(0.0, 3, 3), fill(0.0, 3, 3), fill(0.0, 3, 3)]
    _nvalid = fill(false, 3)
    _T = fill(0.0, __nn * __ndof, __nn * __ndof)
    _elmat = fill(0.0, __nn * __ndof, __nn * __ndof)
    _gradN_e = fill(0.0, __nn, 2)
    _Bm = fill(0.0, 3, __nn * __ndof)
    _Bb = fill(0.0, 3, __nn * __ndof)
    _Bs = fill(0.0, 2, __nn * __ndof)
    _DpsBmb = similar(_Bm)
    _DtBs = similar(_Bs)

    @assert delegateof(integdomain.fes) === FESetShellT3()

    return FEMMShellT3FFComp(integdomain, layup_groups,
        __TRANSV_SHEAR_FORMULATION_AVERAGE_B, 1.0, 30.0, 5 / 12 / 1.5,
        false, _normals, _normal_valid,
        _layup_group_lookup,
        _loc, _J0,
        _ecoords, _edisp, _ecoords_e, _edisp_e,
        _dofnums,
        _E_G, _A_Es, _nvalid, _T,
        _elmat,
        _gradN_e,
        _Bm, _Bb, _Bs, _DpsBmb, _DtBs)
end

"""
    make(integdomain, layup)

Make a T3FFComp FEMM from the integration domain and a composite layup.
"""
function make(integdomain, layup)
    return FEMMShellT3FFComp(integdomain, layup)
end

function _compute_nodal_normal!(n, mcsys::CSys, XYZ, J0::FFltMat, labl::FInt)
    updatecsmat!(mcsys, reshape(XYZ, 1, 3), J0, labl);
    n[:] .= mcsys.csmat[:, 3]
    return n
end

function _compute_J!(J0, ecoords)
    # Compute the Jacobian matrix: the Jacobian matrix is constant within the
    # triangle
    # x, y, z = ecoords[2, :].-ecoords[1, :]
    # J0[:, 1] .= (x, y, z)
    # x, y, z = ecoords[3, :].-ecoords[1, :]
    # J0[:, 2] .= (x, y, z)
    for j in 1:3
        J0[j, 1] = ecoords[2, j] - ecoords[1, j]
        J0[j, 2] = ecoords[3, j] - ecoords[1, j]
    end
    return J0
end

function _e_g!(E_G, J0)
    # J0[:, 1] is the tangent to the coordinate curve 1
    for j in 1:3
        E_G[j,1]  = J0[j, 1]
    end
    n = sqrt(E_G[1, 1]^2 + E_G[2, 1]^2 + E_G[3, 1]^2)
    for j in 1:3
        E_G[j,1]  /= n
    end
    # J0[:, 2] is the tangent to the coordinate curve 2
    # Now compute the normal
    E_G[1, 3] = -E_G[3, 1]*J0[2, 2]+E_G[2, 1]*J0[3, 2]
    E_G[2, 3] =  E_G[3, 1]*J0[1, 2]-E_G[1, 1]*J0[3, 2]
    E_G[3, 3] = -E_G[2, 1]*J0[1, 2]+E_G[1, 1]*J0[2, 2]
    n = sqrt(E_G[1, 3]^2 + E_G[2, 3]^2 + E_G[3, 3]^2);
    for j in 1:3
        E_G[j, 3]  /= n
    end
    E_G[1, 2] = -E_G[3, 3]*E_G[2, 1]+E_G[2, 3]*E_G[3, 1]
    E_G[2, 2] =  E_G[3, 3]*E_G[1, 1]-E_G[1, 3]*E_G[3, 1]
    E_G[3, 2] = -E_G[2, 3]*E_G[1, 1]+E_G[1, 3]*E_G[2, 1]
    return E_G
end

function _ecoords_e!(ecoords_e, J0, E_G)
    ecoords_e[1, 1] = ecoords_e[1, 2] = 0.0
    ecoords_e[2, 1] = dot(view(J0, :, 1), view(E_G, :, 1))
    ecoords_e[2, 2] = dot(view(J0, :, 1), view(E_G, :, 2))
    ecoords_e[3, 1] = dot(view(J0, :, 2), view(E_G, :, 1))
    ecoords_e[3, 2] = dot(view(J0, :, 2), view(E_G, :, 2))
    return ecoords_e
end

function _centroid!(centroid, ecoords)
    centroid[1]  = (ecoords[1, 1]+ecoords[2, 1]+ecoords[3, 1])/3
    centroid[2]  = (ecoords[1, 2]+ecoords[2, 2]+ecoords[3, 2])/3
    centroid[3]  = (ecoords[1, 3]+ecoords[2, 3]+ecoords[3, 3])/3
    return centroid
end

struct NodalTriadsE
    r::FFltVec
    nk_e::FFltVec
    nk::FFltVec
    f3_e::FFltVec
end

function NodalTriadsE()
    NodalTriadsE(fill(0.0, 3), fill(0.0, 3), fill(0.0, 3), [0.0, 0.0, 1.0])
end

(o::NodalTriadsE)(A_Es, nvalid, E_G, normals, normal_valid, c) = begin
    # Components of nodal cartesian ordinate systems such that the third
    # direction is the direction of the nodal normal, and the angle to rotate
    # the element normal into the nodal normal is as short as possible; these
    # components are given on the local element coordinate system basis
    # vectors (basis E). If the angle of rotation is excessively small, the
    # nodal matrix is the identity.
    #
    # The array `A_Es` holds the triads that are the nodal-to-element
    # transformation matrices. The array `nvalid` indicates for the three nodes
    # which of the normals is valid (`true`).

    for k in 1:length(A_Es)
        o.nk .= vec(view(normals, c[k], :))
        nvalid[k] = normal_valid[c[k]]
        if nvalid[k]
            # If the nodal normal is valid, pull it back into the element frame.
            mul!(o.nk_e, transpose(E_G), o.nk)
        else
            # Otherwise the element normal replaces the nodal normal.
            o.nk_e .= 0.0; o.nk_e[3] = 1.0
        end
        cross3!(o.r, o.f3_e, o.nk_e)
        if norm(o.r) > 1.0e-12
            rotmat3!(A_Es[k], o.r)
        else
            # A_Es[k] .= I(3)
            A_Es[k] .= 0.0
            A_Es[k][1, 1] = A_Es[k][2, 2] = A_Es[k][3, 3] = 1.0
        end
    end
    return A_Es, nvalid
end

struct TransfmatGToA
    Tblock::FFltMat
end

function TransfmatGToA()
    TransfmatGToA(fill(0.0, 3, 3))
end

(o::TransfmatGToA)(T, A_Es, E_G) = begin
    # Global-to-nodal transformation matrix.

    # The 3x3 blocks consist of the nodal triad expressed on the global basis.
    # The nodal basis vectors in the array `A_Es[i]` are expressed on the
    # element basis, and are then rotated with `E_G` into the global coordinate
    # system.

    # Output
    # - `T` = transformation matrix, input in the global basis, output in the
    #   nodal basis
    T .= 0.0
    for i in 1:__nn
        mul!(o.Tblock,  transpose(A_Es[i]), transpose(E_G))
        offset = (i-1)*__ndof
        r = offset+1:offset+3
        @. T[r, r] = o.Tblock
        r = offset+4:offset+6
        @. T[r, r] = o.Tblock
    end
    return T
end

function _transfmat_a_to_e!(T, A_Es, gradN_e)
    # Nodal-to-element transformation matrix.
    # - `A_Es` = matrix with the nodal triad vectors in columns, components on
    #   the element basis. Its transpose is the element triad on the nodal
    #   basis vectors.
    # - `gradN_e` = basis function gradients on the element basis.
    # Output
    # - `T` = transformation matrix, input in the nodal basis, output in the
    #   element basis: t_e = T * t_n
    # Rotation degrees of freedom: The drilling rotation of the mid surface
    # produced by the 1/2*(v,x - u,y) effect is linked to the out of plane
    # rotations.

    T .= 0.0
    for i in 1:__nn
        roffst = (i-1)*__ndof
        iA_E = A_Es[i]
        iA_E_33 = iA_E[3, 3]
        # T[r, r] .= iA_E
        for cl in 1:3
            for rw in 1:3
                T[roffst+rw, roffst+cl] = iA_E[rw, cl]
            end
        end
        for cl in 1:2
            for rw in 1:2
                T[roffst+3+rw, roffst+3+cl] = iA_E[rw, cl] - (1/iA_E_33)*iA_E[rw, 3]*iA_E[cl, 3]
            end
        end
        m1 = (1/iA_E_33)*iA_E[1, 3]
        m2 = (1/iA_E_33)*iA_E[2, 3]
        for j in 1:__nn
            coffst = (j-1)*__ndof
            for k in 1:3
                a3 = 1/2 * (iA_E[2, k] * gradN_e[j, 1] - iA_E[1, k] * gradN_e[j, 2])
                T[roffst+4, coffst+k] += m1 * a3
                T[roffst+5, coffst+k] += m2 * a3
            end
        end
    end
    return T
end

function _gradN_e_Ae!(gradN_e, ecoords_e)
    # Compute the gradient with respect to the element coordinates, three
    # gradients in two coordinates, and the area of the triangle
    # a, b = ecoords_e[2, :] .- ecoords_e[1, :]
    # c, d = ecoords_e[3, :] .- ecoords_e[1, :]
    # J = (a*d - b*c)
    # gradN_e[:] .= ((b-d)/J, d/J, -b/J, (c-a)/J, -c/J, a/J)
     a = ecoords_e[2, 1] - ecoords_e[1, 1]
     b = ecoords_e[2, 2] - ecoords_e[1, 2]
     c = ecoords_e[3, 1] - ecoords_e[1, 1]
     d = ecoords_e[3, 2] - ecoords_e[1, 2]
     J = (a*d - b*c)
     gradN_e[1, 1] = (b-d)/J
     gradN_e[2, 1] = d/J
     gradN_e[3, 1] = -b/J
     gradN_e[1, 2] = (c-a)/J
     gradN_e[2, 2] = -c/J
     gradN_e[3, 2] = a/J
    return gradN_e, J/2
end

function _add_Bsmat_o!(Bs, ecoords_e, Ae, ordering)
    # Compute the linear transverse shear strain-displacement matrix for one
    # particular ordering of the nodes. The computed entries are ADDED. The
    # matrix is NOT zeroed out initially. That is the responsibility of the
    # caller.

    # Orientation
    s, p, q = ordering
    a = ecoords_e[p, 1] - ecoords_e[s, 1]
    b = ecoords_e[p, 2] - ecoords_e[s, 2]
    c = ecoords_e[q, 1] - ecoords_e[s, 1]
    d = ecoords_e[q, 2] - ecoords_e[s, 2]
    m = (1/2/Ae) # multiplier

    # The first node in the triangle
    # Node s
    co = (s - 1) * 6 # column offset
    Bs[1, co+3] += m*(b-d);                             Bs[1, co+5] += m*(Ae)
    Bs[2, co+3] += m*(c-a); Bs[2, co+4] += m*(-Ae);
    # The other two nodes
    # Node p
    co = (p - 1) * 6 # column offset
    Bs[1, co+3] += m*(d);   Bs[1, co+4] += m*(-b*d/2);  Bs[1, co+5] += m*(a*d/2)
    Bs[2, co+3] += m*(-c);  Bs[2, co+4] += m*(b*c/2);   Bs[2, co+5] += m*(-a*c/2)
    # Node q
    co = (q - 1) * 6 # column offset
    Bs[1, co+3] += m*(-b);  Bs[1, co+4] += m*(b*d/2);   Bs[1, co+5] += m*(-b*c/2)
    Bs[2, co+3] += m*(a);   Bs[2, co+4] += m*(-a*d/2);  Bs[2, co+5] += m*(a*c/2)

    return Bs
end

function _Bsmat!(Bs, ecoords_e, Ae)
    # Compute the linear transverse shear strain-displacement matrix.
    Bs .= 0.0 # Zero out initially, then add the three contributions
    _add_Bsmat_o!(Bs, ecoords_e, Ae, (1, 2, 3))
    _add_Bsmat_o!(Bs, ecoords_e, Ae, (2, 3, 1))
    _add_Bsmat_o!(Bs, ecoords_e, Ae, (3, 1, 2))
    # Now averaging of the three contributions
    Bs .*= (1/3)
    return Bs
end

function _Bmmat!(Bm, gradN)
    # Compute the linear membrane strain-displacement matrix.
    fill!(Bm, 0.0)
    for i in 1:__nn
        Bm[1,6*(i-1)+1] = gradN[i,1];
                                             Bm[2,6*(i-1)+2] = gradN[i,2];
        Bm[3,6*(i-1)+1] = gradN[i,2];        Bm[3,6*(i-1)+2] = gradN[i,1];
    end
end

function _Bbmat!(Bb, gradN)
    # Compute the linear, displacement independent, curvature-displacement/rotation
    # matrix for a shell quadrilateral element with nfens=3 nodes. Displacements and
    # rotations are in a local coordinate system.
    fill!(Bb, 0.0)
    for i in 1:__nn
                                              Bb[1,6*(i-1)+5] = gradN[i,1];
        Bb[2,6*(i-1)+4] = -gradN[i,2];
        Bb[3,6*(i-1)+4] = -gradN[i,1];        Bb[3,6*(i-1)+5] = gradN[i,2];
    end
end

"""
    associategeometry!(self::FEMMShellT3FFComp,  geom::NodalField{FFlt})

Associate geometry with the FEMM.

In this case it means evaluate the nodal normals.
"""
function associategeometry!(self::FEMMShellT3FFComp,  geom::NodalField{FFlt})
    threshold_angle = self.threshold_angle
    J0 = self._J0
    E_G = self._E_G
    normals, normal_valid = self._normals, self._normal_valid
    nnormal = fill(0.0, 3)

    # Compute the normals at the nodes
    for lg in self.layup_groups
        layup = lg[1]
        eset = lg[2]
        for el in eset
            i, j, k = self.integdomain.fes.conn[el]
            J0[:, 1] = geom.values[j, :] - geom.values[i, :]
            J0[:, 2] = geom.values[k, :] - geom.values[i, :]
            for n in [i, j, k]
                _compute_nodal_normal!(nnormal, layup.csys, geom.values[n, :], J0, self.integdomain.fes.label[el])
                normals[n, :] .+= nnormal
            end
        end
    end

    # Normalize to unit length
    for j in 1:size(normals, 1)
        nn = norm(normals[j, :])
        if nn > 0.0
            normals[j, :] ./= nn
        end
    end

    # Now perform a second pass. If the nodal normal differs by a substantial
    # amount from the element normals, zero out the nodal normal to indicate
    # that the nodal normal should not be used at that vertex.
    ntolerance = 1 - sqrt(1 - sin(threshold_angle/180*pi)^2)
    for el in 1:count(self.integdomain.fes)
        i, j, k = self.integdomain.fes.conn[el]
        _compute_J!(J0, geom.values[[i, j, k], :])
        for n in [i, j, k]
            _e_g!(E_G, J0)
            nd = dot(normals[n, :], E_G[:, 3])
            if nd < 1 - ntolerance
                normal_valid[n] = false
            end
        end
    end

    # Mark the engine as being associated with geometry
    self._associatedgeometry = true
    return self
end

"""
    num_normals(self::FEMMShellT3FFComp)

Compute the summary of the nodal normals.
"""
function num_normals(self::FEMMShellT3FFComp)
    normals, normal_valid = self._normals, self._normal_valid
    total_normals = length(normal_valid)
    total_invalid_normals = length([v for v in normal_valid if v != true])
    return total_normals, total_invalid_normals
end

"""
    stiffness(self::FEMMShellT3FFComp, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{TI}, Rfield1::NodalField{TI}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, TI<:Number}

Compute the material stiffness matrix.
"""
function stiffness(self::FEMMShellT3FFComp, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{TI}, Rfield1::NodalField{TI}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, TI<:Number}
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    label = self.integdomain.fes.label
    normals, normal_valid = self._normals, self._normal_valid
    centroid, J0 = self._loc, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e, gradN_e = self._ecoords_e, self._gradN_e
    E_G, A_Es, nvalid, T = self._E_G, self._A_Es, self._nvalid, self._T
    elmat = self._elmat
    transformwith = QTEQTransformer(elmat)
    lla = Layup2ElementAngle()
    _nodal_triads_e! = NodalTriadsE()
    _transfmat_g_to_a! = TransfmatGToA()
    Bm, Bb, Bs, DpsBmb, DtBs = self._Bm, self._Bb, self._Bs, self._DpsBmb, self._DtBs
    A, B, C, BT = zeros(3, 3), zeros(3, 3), zeros(3, 3), zeros(3, 3)
    H = zeros(2, 2)
    sA, sB, sC = zeros(3, 3), zeros(3, 3), zeros(3, 3)
    sH = zeros(2, 2)
    Tps, Tts = zeros(3, 3), zeros(2, 2)
    tps! = QEQTTransformer(Tps)
    tts! = QEQTTransformer(Tts)
    drilling_stiffness_scale = self.drilling_stiffness_scale
    transv_shear_formulation = self.transv_shear_formulation
    mult_el_size = self.mult_el_size
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for lg in self.layup_groups
        layup = lg[1]
        eset = lg[2]
        t = thickness(layup)
        laminate_stiffnesses!(layup, A, B, C)
        laminate_transverse_stiffness!(layup, H)
        for i in eset # Loop over elements in the layup group
            gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
            _centroid!(centroid, ecoords)
            _compute_J!(J0, ecoords)
            _e_g!(E_G, J0)
            _ecoords_e!(ecoords_e, J0, E_G)
            gradN_e, Ae = _gradN_e_Ae!(gradN_e, ecoords_e)
            # Working copies to be transformed
            sA[:] .= A[:]; sB[:] .= B[:]; sC[:] .= C[:];     sH[:] .= H[:]
            # Transform the laminate stiffnesses
            updatecsmat!(layup.csys, reshape(centroid, 1, 3), J0, -1);
            m, n = lla(E_G, layup.csys.csmat)
            plane_stress_T_matrix!(Tps, m, -n)
            tps!(sA, Tps); tps!(sB, Tps); tps!(sC, Tps);
            transverse_shear_T_matrix!(Tts, m, n)
            tts!(sH, Tts)
            # Construct the Stiffness Matrix
            fill!(elmat,  0.0); # Initialize element matrix
            _Bmmat!(Bm, gradN_e)
            add_btdb_ut_only!(elmat, Bm, Ae, sA, DpsBmb)
            _Bbmat!(Bb, gradN_e)
            add_btdb_ut_only!(elmat, Bb, Ae, sC, DpsBmb)
            add_b1tdb2!(elmat, Bm, Bb, Ae, sB, DpsBmb)
            add_b1tdb2!(elmat, Bb, Bm, Ae, sB, DpsBmb)
            # he = sqrt(2*Ae) # we avoid taking the square root here, replacing
            # he^2 with 2*Ae
            if transv_shear_formulation == __TRANSV_SHEAR_FORMULATION_AVERAGE_K
                for  o in [(1, 2, 3), (2, 3, 1), (3, 1, 2)]
                    Bs .= 0.0; _add_Bsmat_o!(Bs, ecoords_e, Ae, o)
                    add_btdb_ut_only!(elmat, Bs, (t^2/(t^2+mult_el_size*2*Ae))*Ae/3, sH, DtBs)
                end
            else
                _Bsmat!(Bs, ecoords_e, Ae)
                add_btdb_ut_only!(elmat, Bs, (t^2/(t^2+mult_el_size*2*Ae))*Ae, sH, DtBs)
            end
            # Complete the elementwise matrix by filling in the lower triangle
            complete_lt!(elmat)
            # Transformation from the nodal (A) to the element (E) basis
            _nodal_triads_e!(A_Es, nvalid, E_G, normals, normal_valid, fes.conn[i])
            _transfmat_a_to_e!(T, A_Es, gradN_e)
            transformwith(elmat, T)
            # Bending diagonal stiffness coefficients
            kavg = mean((elmat[4, 4], elmat[10, 10], elmat[16, 16],
                elmat[5, 5], elmat[11, 11], elmat[17, 17])) * drilling_stiffness_scale
            # Add the artificial drilling stiffness in the nodal basis, but only
            # if the nodal normal is valid.
            nvalid[1] && (elmat[6,6] += kavg)
            nvalid[2] && (elmat[12,12] += kavg)
            nvalid[3] && (elmat[18,18] += kavg)
            # Transform from global (G) into the nodal (A) basis
            _transfmat_g_to_a!(T, A_Es, E_G)
            transformwith(elmat, T)
            # Assembly
            gatherdofnums!(dchi, dofnums, fes.conn[i]);
            assemble!(assembler, elmat, dofnums, dofnums);
        end # Loop over elements
    end # Loop over layup groups
    return makematrix!(assembler);
end

function stiffness(self::FEMMShellT3FFComp, geom0::NodalField{FFlt}, u1::NodalField{TI}, Rfield1::NodalField{TI}, dchi::NodalField{TI}) where {TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom0, u1, Rfield1, dchi);
end

"""
    mass(self::FEMMShellT3FFComp,  assembler::A,  geom::NodalField{FFlt}, dchi::NodalField{TI}) where {A<:AbstractSysmatAssembler, TI<:Number}

Compute the diagonal (lumped) mass matrix

The mass matrix can be expected to be non-singular.
"""
function mass(self::FEMMShellT3FFComp,  assembler::A,  geom0::NodalField{FFlt}, dchi::NodalField{TI}) where {A<:AbstractSysmatAssembler, TI<:Number}
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    label = self.integdomain.fes.label
    normals, normal_valid = self._normals, self._normal_valid
    centroid, J0 = self._loc, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e, gradN_e = self._ecoords_e, self._gradN_e
    E_G, A_Es, nvalid, T = self._E_G, self._A_Es, self._nvalid, self._T
    elmat = self._elmat
    npe = nodesperelem(fes)
    ndn = ndofs(dchi)
    startassembly!(assembler,  size(elmat,1),  size(elmat,2),  count(fes), dchi.nfreedofs,  dchi.nfreedofs);
    for lg in self.layup_groups
            layup = lg[1]
            eset = lg[2]
            mass_density, moment_of_inertia_density = laminate_inertia!(layup)
            for i in eset # Loop over elements in the layup group
                gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
                _centroid!(centroid, ecoords)
                _compute_J!(J0, ecoords)
                _e_g!(E_G, J0)
                _ecoords_e!(ecoords_e, J0, E_G)
                gradN_e, Ae = _gradN_e_Ae!(gradN_e, ecoords_e)
                # Compute the translational and rotational masses corresponding
                # to nodes
                tmass = mass_density*(Ae/3);
                rmass = moment_of_inertia_density*(Ae/3);
                # Fill the elementwise matrix in the global basis: the matrix is
                # lumped, hence diagonal, and no transformation is needed.
                fill!(elmat,  0.0); # Initialize element matrix
                for k in 1:npe
                    # Translation degrees of freedom
                    for d in 1:3
                        c = (k - 1) * __ndof + d
                        elmat[c, c] += tmass
                    end
                    # Bending degrees of freedom
                    for d in 4:6
                        c = (k - 1) * __ndof + d
                        elmat[c, c] += rmass
                    end
                end
                # Assemble
                gatherdofnums!(dchi,  dofnums,  fes.conn[i]);# retrieve degrees of freedom
                assemble!(assembler,  elmat,  dofnums,  dofnums);# assemble symmetric matrix
            end # Loop over elements
        end # Loop over layup groups
    return makematrix!(assembler);
end

function mass(self::FEMMShellT3FFComp,  geom::NodalField{FFlt},  u::NodalField{TI}) where {TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return mass(self, assembler, geom, u);
end

"""
    inspectintegpoints(self::AbstractFEMMDeforLinear,
      geom::NodalField{FFlt},  u::NodalField{TI},
      dT::NodalField{FFlt},
      felist::FIntVec,
      inspector::F,  idat, quantity=:Cauchy;
      context...) where {TI<:Number, F<:Function}

Inspect integration point quantities.

- `geom` - reference geometry field
- `u` - displacement+rotation field
- `dT` - temperature difference field
- `felist` - indexes of the finite elements that are to be inspected:
     The fes to be included are: `fes[felist]`.
- `context`    - structure: see the update!() method of the material.
- `inspector` - function with the signature
        `idat = inspector(idat, j, conn, x, out, loc);`
   where
    `idat` - a structure or an array that the inspector may
           use to maintain some state,  for instance minimum or maximum of
           stress, `j` is the element number, `conn` is the element connectivity,
           `out` is the output of the update!() method,  `loc` is the location
           of the integration point in the *reference* configuration.
### Return
The updated inspector data is returned.
"""
function inspectintegpoints(self::FEMMShellT3FFComp, geom0::NodalField{FFlt},  u::NodalField{TI}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {TI<:Number, F<:Function}
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    label = self.integdomain.fes.label
    normals, normal_valid = self._normals, self._normal_valid
    centroid, J0 = self._loc, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e, gradN_e = self._ecoords_e, self._gradN_e
    E_G, A_Es, nvalid, T = self._E_G, self._A_Es, self._nvalid, self._T
    elmat = self._elmat
    transformwith = QTEQTransformer(elmat)
    _nodal_triads_e! = NodalTriadsE()
    _transfmat_g_to_a! = TransfmatGToA()
    Bm, Bb, Bs, DpsBmb, DtBs = self._Bm, self._Bb, self._Bs, self._DpsBmb, self._DtBs
    A, B, C, BT = zeros(3, 3), zeros(3, 3), zeros(3, 3), zeros(3, 3)
    H = zeros(2, 2)
    sA, sB, sC = zeros(3, 3), zeros(3, 3), zeros(3, 3)
    sH = zeros(2, 2)
    Tps, Tts = zeros(3, 3), zeros(2, 2)
    tps! = QEQTTransformer(Tps)
    tts! = QEQTTransformer(Tts)
    lla = Layup2ElementAngle()
    mult_el_size = self.mult_el_size
    edisp_e = deepcopy(edisp)
    edisp_n = deepcopy(edisp_e)
    out = fill(0.0, 3)
    o2_e = fill(0.0, 2, 2)
    # Sort out  the output requirements
    outputcsys = self.layup_groups[1][1].csys; # default: report the stresses in the material coord system
    for apair in pairs(context)
        sy, val = apair
        if sy == :outputcsys
            outputcsys = val
        end
    end
    BENDING_MOMENT, TRANSVERSE_SHEAR, MEMBRANE_FORCE = 1, 2, 3
    quant = BENDING_MOMENT
    if quantity == :bending || quantity == :moment || quantity == :bending_moment
        quant = BENDING_MOMENT
    end
    if quantity == :transverse_shear || quantity == :transverse || quantity == :shear
        quant = TRANSVERSE_SHEAR
    end
    if quantity == :membrane_force || quantity == :membrane
        quant = MEMBRANE_FORCE
    end
    # Loop over  all the elements and all the quadrature points within them
    for ilist = 1:length(felist) # Loop over elements
        i = felist[ilist];
        # Look up the layup
        layup = self.layup_groups[self._layup_group_lookup[i]][1]
        t = thickness(layup)
        laminate_stiffnesses!(layup, A, B, C)
        laminate_transverse_stiffness!(layup, H)
        gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
        gathervalues_asvec!(u, edisp, fes.conn[i]);
        _centroid!(centroid, ecoords)
        _compute_J!(J0, ecoords)
        _e_g!(E_G, J0)
        _ecoords_e!(ecoords_e, J0, E_G)
        gradN_e, Ae = _gradN_e_Ae!(gradN_e, ecoords_e)
        # Establish nodal triads
        _nodal_triads_e!(A_Es, nvalid, E_G, normals, normal_valid, fes.conn[i])
        # Transform from global into nodal coordinates
        _transfmat_g_to_a!(T, A_Es, E_G)
        mul!(edisp_n, T, edisp)
        # Now treat the transformation from the nodal to the element triad
        _transfmat_a_to_e!(T, A_Es, gradN_e)
         # Transform the nodal vector into the elementwise coordinates
        mul!(edisp_e, T, edisp_n)
        updatecsmat!(outputcsys, centroid, J0, fes.label[i]);
        if dot(view(outputcsys.csmat, :, 3), view(E_G, :, 3)) < 0.95
            @warn "Coordinate systems mismatched?"
        end
        # Established the stiffness matrices
        sA[:] .= A[:]; sB[:] .= B[:]; sC[:] .= C[:];     sH[:] .= H[:]
        # Transform the laminate stiffnesses
        updatecsmat!(layup.csys, centroid, J0, fes.label[i]);
        m, n = lla(E_G, layup.csys.csmat)
        plane_stress_T_matrix!(Tps, m, -n)
        tps!(sA, Tps); tps!(sB, Tps); tps!(sC, Tps);
        transverse_shear_T_matrix!(Tts, m, n)
        tts!(sH, Tts)
        # The output coordinate system
        ocsm, ocsn = lla(E_G, outputcsys.csmat)
        o2_e[1, 1] = o2_e[2, 2]  = ocsm
        o2_e[1, 2] = ocsn; o2_e[2, 1] = -ocsn
        # Compute the Requested Quantity
        _Bbmat!(Bb, gradN_e)
        _Bmmat!(Bm, gradN_e)
        kurv = Bb * edisp_e
        memstr = Bm * edisp_e
        if quant == BENDING_MOMENT
            mom = sB * memstr + sC * kurv
            m = [mom[1] mom[3]; mom[3] mom[2]]
            mo = o2_e' * m * o2_e
            out[:] .= mo[1, 1], mo[2, 2], mo[1, 2]
        end
        if quant == MEMBRANE_FORCE
            frc = sA * memstr + sB * kurv
            f = [frc[1] frc[3]; frc[3] frc[2]]
            fo = o2_e' * f * o2_e
            out[:] .= fo[1, 1], fo[2, 2], fo[1, 2]
            # @infiltrate
        end
        if quant == TRANSVERSE_SHEAR
            _Bsmat!(Bs, ecoords_e, Ae)
            he = sqrt(2*Ae)
            shrstr = Bs * edisp_e
            frc = (t^2/(t^2+mult_el_size*2*Ae))*sH * shrstr
            fo = o2_e' * frc
            out[1:2] .= fo[1], fo[2]
        end
        # Call the inspector
        idat = inspector(idat, i, fes.conn[i], ecoords, out, centroid);
        # end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

function inspectintegpoints(self::FEMMShellT3FFComp, geom::NodalField{FFlt},  u::NodalField{TI}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {TI<:Number, F<:Function}
    dT = NodalField(fill(zero(FFlt), nnodes(geom), 1)) # zero difference in temperature
    return inspectintegpoints(self, geom, u, dT, felist, inspector, idat, quantity; context...);
end

end # module


text:
> z = 4

>1z + 4
 8

>x = 3
 3

>1 < x < 5
 true

> [a 2 3;
       4 5 6;
       7 8 9]
3×3 Array{Int64,2}:
 2  2  3
 4  5  6
 7  8  9

max_stars_repo_name: kyoujuro/Julia-practice

max_issues_repo_name: kyoujuro/Julia-practice

max_forks_repo_name: kyoujuro/Julia-practice


text:
"""
    convolve(d1::T, d2::T) where T<:Distribution -> Distribution

Convolve two distributions of the same type to yield the distribution corresponding to the
sum of independent random variables drawn from the underlying distributions.

The function is only defined in the cases where the convolution has a closed form as
defined here https://en.wikipedia.org/wiki/List_of_convolutions_of_probability_distributions

* `Bernoulli`
* `Binomial`
* `NegativeBinomial`
* `Geometric`
* `Poisson`
* `Normal`
* `Cauchy`
* `Chisq`
* `Exponential`
* `Gamma`
* `MultivariateNormal`
"""
function convolve end

# discrete univariate
function convolve(d1::Bernoulli, d2::Bernoulli)
    _check_convolution_args(d1.p, d2.p)
    return Binomial(2, d1.p)
end

function convolve(d1::Binomial, d2::Binomial)
    _check_convolution_args(d1.p, d2.p)
    return Binomial(d1.n + d2.n, d1.p)
end

function convolve(d1::NegativeBinomial, d2::NegativeBinomial)
    _check_convolution_args(d1.p, d2.p)
    return NegativeBinomial(d1.r + d2.r, d1.p)
end

function convolve(d1::Geometric, d2::Geometric)
    _check_convolution_args(d1.p, d2.p)
    return NegativeBinomial(2, d1.p)
end

convolve(d1::Poisson, d2::Poisson) =  Poisson(d1.λ + d2.λ)

# continuous univariate
convolve(d1::Normal, d2::Normal) = Normal(d1.μ + d2.μ, hypot(d1.σ, d2.σ))
convolve(d1::Cauchy, d2::Cauchy) = Cauchy(d1.μ + d2.μ, d1.σ + d2.σ)
convolve(d1::Chisq, d2::Chisq) = Chisq(d1.ν + d2.ν)

function convolve(d1::Exponential, d2::Exponential)
    _check_convolution_args(d1.θ, d2.θ)
    return Gamma(2, d1.θ)
end

function convolve(d1::Gamma, d2::Gamma)
    _check_convolution_args(d1.θ, d2.θ)
    return Gamma(d1.α + d2.α, d1.θ)
end

# continuous multivariate
# The first two methods exist for performance reasons to avoid unnecessarily converting
# PDMats to a Matrix
function convolve(
    d1::Union{IsoNormal, ZeroMeanIsoNormal, DiagNormal, ZeroMeanDiagNormal},
    d2::Union{IsoNormal, ZeroMeanIsoNormal, DiagNormal, ZeroMeanDiagNormal},
    )
    _check_convolution_shape(d1, d2)
    return MvNormal(d1.μ .+ d2.μ, d1.Σ + d2.Σ)
end

function convolve(
    d1::Union{FullNormal, ZeroMeanFullNormal},
    d2::Union{FullNormal, ZeroMeanFullNormal},
    )
    _check_convolution_shape(d1, d2)
    return MvNormal(d1.μ .+ d2.μ, d1.Σ.mat + d2.Σ.mat)
end

function convolve(d1::MvNormal, d2::MvNormal)
    _check_convolution_shape(d1, d2)
    return MvNormal(d1.μ .+ d2.μ, Matrix(d1.Σ) + Matrix(d2.Σ))
end

function _check_convolution_args(p1, p2)
    p1 ≈ p2 || throw(ArgumentError(
    "$(p1) !≈ $(p2): distribution parameters must be approximately equal",
    ))
end

function _check_convolution_shape(d1, d2)
    length(d1) == length(d2) || throw(ArgumentError("$d1 and $d2 are not the same size"))
end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ 674905ed-75a1-457c-a866-279cdb76d243
#hide
begin
#for sharing environments between notebooks
import Pkg
Pkg.activate(Base.current_project())
end

# ╔═╡ 2c396254-932f-4431-8f94-08dd2766aa8e
#hide
include("../src/Export.jl")

# ╔═╡ 17a13720-29ba-11ec-0ba4-83c9ee9ac582
#export
include("../src/Documenter.jl")

# ╔═╡ 1c5d5f39-0349-454b-aa22-5feda1287b9d
#export
function build(nbsdir)
		Documenter.export2md(nbsdir)
		Documenter.export2readme()
end

# ╔═╡ c96a7fe1-ccba-4499-aa15-95ced972492b
#hide
Export.notebook2script(joinpath("..", "nbs"), joinpath("..", "src"))

# ╔═╡ Cell order:
# ╠═674905ed-75a1-457c-a866-279cdb76d243
# ╠═2c396254-932f-4431-8f94-08dd2766aa8e
# ╠═17a13720-29ba-11ec-0ba4-83c9ee9ac582
# ╠═1c5d5f39-0349-454b-aa22-5feda1287b9d
# ╠═c96a7fe1-ccba-4499-aa15-95ced972492b