text:
# PLANAR MANIPULATIONS:
# =========================================================================================#
"""
    split(f::FusionTree{I, N}, M::Int)
    -> (::FusionTree{I, M}, ::FusionTree{I, N-M+1})

Split a fusion tree `f` into `f1` and `f2` between the `M`th and `M+1`th uncoupled sector,
and return `f1, f2`.

The first tree `f1` has uncoupled sectors as the first `M` uncoupled sectors of the input
tree `f`, and a coupled sector as the `M-1`th internal sector of `f`.

The second tree `f2` has the first uncoupled sector as the coupled sector of `f1`, and the
remaining uncoupled sectors as the remaining `N-M` uncoupled sectors of `f`. Its coupled
sector is the same with that of `f`.
"""
@inline function split(f::FusionTree{I, N}, M::Int) where {I, N}
    if M > N || M < 0
        throw(ArgumentError("M should be between 0 and N = $N"))
    elseif M === N
        (f, FusionTree{I}((f.coupled,), f.coupled, (false,), (), ()))
    elseif M === 1
        isdual1 = (f.isdual[1],)
        isdual2 = Base.setindex(f.isdual, false, 1)
        f1 = FusionTree{I}((f.uncoupled[1],), f.uncoupled[1], isdual1, (), ())
        f2 = FusionTree{I}(f.uncoupled, f.coupled, isdual2, f.innerlines, f.vertices)
        return f1, f2
    elseif M === 0
        f1 = FusionTree{I}((), one(I), (), ())
        uncoupled2 = (one(I), f.uncoupled...)
        coupled2 = f.coupled
        isdual2 = (false, f.isdual...)
        innerlines2 = N >= 2 ? (f.uncoupled[1], f.innerlines...) : ()
        if FusionStyle(I) isa GenericFusion
            vertices2 = (1, f.vertices...)
            return f1, FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2, vertices2)
        else
            return f1, FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2)
        end
    else
        uncoupled1 = ntuple(n->f.uncoupled[n], M)
        isdual1 = ntuple(n->f.isdual[n], M)
        innerlines1 = ntuple(n->f.innerlines[n], max(0, M-2))
        coupled1 = f.innerlines[M-1]
        vertices1 = ntuple(n->f.vertices[n], M-1)

        uncoupled2 = ntuple(N - M + 1) do n
            n == 1 ? f.innerlines[M - 1] : f.uncoupled[M + n - 1]
        end
        isdual2 = ntuple(N - M + 1) do n
            n == 1 ? false : f.isdual[M + n - 1]
        end
        innerlines2 = ntuple(n->f.innerlines[M-1+n], N-M-1)
        coupled2 = f.coupled
        vertices2 = ntuple(n->f.vertices[M-1+n], N-M)

        f1 = FusionTree{I}(uncoupled1, coupled1, isdual1, innerlines1, vertices1)
        f2 = FusionTree{I}(uncoupled2, coupled2, isdual2, innerlines2, vertices2)
        return f1, f2
    end
end

"""
    insertat(f::FusionTree{I, N₁}, i::Int, f2::FusionTree{I, N₂})
    -> <:AbstractDict{<:FusionTree{I, N₁+N₂-1}, <:Number}

Attach a fusion tree `f2` to the uncoupled leg `i` of the fusion tree `f1` and bring it
into a linear combination of fusion trees in canonical form. This requires that
`f2.coupled == f1.uncoupled[i]` and `f1.isdual[i] == false`.
"""
function insertat(f1::FusionTree{I}, i::Int, f2::FusionTree{I, 0}) where {I}
    # this actually removes uncoupled line i, which should be trivial
    (f1.uncoupled[i] == f2.coupled && !f1.isdual[i]) ||
        throw(SectorMismatch("cannot connect $(f2.uncoupled) to $(f1.uncoupled[i])"))
    coeff = Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1,1,1,1]
    uncoupled = TupleLabXD.deleteat(f1.uncoupled, i)
    coupled = f1.coupled
    isdual = TupleLabXD.deleteat(f1.isdual, i)
    if length(uncoupled) <= 2
        inner = ()
    else
        inner = TupleLabXD.deleteat(f1.innerlines, max(1, i-2))
    end
    if length(uncoupled) <= 1
        vertices = ()
    else
        vertices = TupleLabXD.deleteat(f1.vertices, max(1, i-1))
    end
    f = FusionTree(uncoupled, coupled, isdual, inner, vertices)
    return fusiontreedict(I)(f => coeff)
end
function insertat(f1::FusionTree{I}, i::Int, f2::FusionTree{I, 1}) where {I}
    # identity operation
    (f1.uncoupled[i] == f2.coupled && !f1.isdual[i]) ||
        throw(SectorMismatch("cannot connect $(f2.uncoupled) to $(f1.uncoupled[i])"))
    coeff = Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1,1,1,1]
    isdual′ = TupleLabXD.setindex(f1.isdual, f2.isdual[1], i)
    f = FusionTree{I}(f1.uncoupled, f1.coupled, isdual′, f1.innerlines, f1.vertices)
    return fusiontreedict(I)(f => coeff)
end
function insertat(f1::FusionTree{I}, i::Int, f2::FusionTree{I, 2}) where {I}
    # elementary building block
    (f1.uncoupled[i] == f2.coupled && !f1.isdual[i]) ||
        throw(SectorMismatch("cannot connect $(f2.uncoupled) to $(f1.uncoupled[i])"))
    uncoupled = f1.uncoupled
    coupled = f1.coupled
    inner = f1.innerlines
    b, c = f2.uncoupled
    isdual = f1.isdual
    isdualb, isdualc = f2.isdual
    if i == 1
        uncoupled′ = (b, c, tail(uncoupled)...)
        isdual′ = (isdualb, isdualc, tail(isdual)...)
        inner′ = (uncoupled[1], inner...)
        vertices′ = (f2.vertices..., f1.vertices...)
        coeff = Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1,1,1,1]
        f′ = FusionTree(uncoupled′, coupled, isdual′, inner′, vertices′)
        return fusiontreedict(I)(f′ => coeff)
    end
    uncoupled′ = TupleLabXD.insertafter(Base.setindex(uncoupled, b, i), i, (c,))
    isdual′ = TupleLabXD.insertafter(Base.setindex(isdual, isdualb, i), i, (isdualc,))
    a = i == 2 ? uncoupled[1] : inner[i-2]
    d = i == length(f1) ? coupled : inner[i-1]
    e′ = uncoupled[i]
    if FusionStyle(I) isa MultiplicityFreeFusion
        local newtrees
        for e in a ⊗ b
            coeff = conj(Fsymbol(a, b, c, d, e, e′))
            iszero(coeff) && continue
            inner′ = TupleLabXD.insertafter(inner, i-2, (e,))
            f′ = FusionTree(uncoupled′, coupled, isdual′, inner′)
            if @isdefined newtrees
                push!(newtrees, f′=> coeff)
            else
                newtrees = fusiontreedict(I)(f′ => coeff)
            end
        end
        return newtrees
    else
        local newtrees
        κ = f2.vertices[1]
        λ = f1.vertices[i-1]
        for e in a ⊗ b
            inner′ = TupleLabXD.insertafter(inner, i-2, (e,))
            Fmat = Fsymbol(a, b, c, d, e, e′)
            for μ = 1:size(Fmat, 1), ν = 1:size(Fmat, 2)
                coeff = conj(Fmat[μ,ν,κ,λ])
                iszero(coeff) && continue
                vertices′ = TupleLabXD.setindex(f1.vertices, ν, i-1)
                vertices′ = TupleLabXD.insertafter(vertices′, i-2, (μ,))
                f′ = FusionTree(uncoupled′, coupled, isdual′, inner′, vertices′)
                if @isdefined newtrees
                    push!(newtrees, f′=> coeff)
                else
                    newtrees = fusiontreedict(I)(f′ => coeff)
                end
            end
        end
        return newtrees
    end
end
function insertat(f1::FusionTree{I,N₁}, i::Int, f2::FusionTree{I,N₂}) where {I,N₁,N₂}
    F = fusiontreetype(I, N₁ + N₂ - 1)
    (f1.uncoupled[i] == f2.coupled && !f1.isdual[i]) ||
        throw(SectorMismatch("cannot connect $(f2.uncoupled) to $(f1.uncoupled[i])"))
    coeff = Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1,1,1,1]
    T = typeof(coeff)
    if length(f1) == 1
        return fusiontreedict(I){F,T}(f2 => coeff)
    end
    if i == 1
        uncoupled = (f2.uncoupled..., tail(f1.uncoupled)...)
        isdual = (f2.isdual..., tail(f1.isdual)...)
        inner = (f2.innerlines..., f2.coupled, f1.innerlines...)
        vertices = (f2.vertices..., f1.vertices...)
        coupled = f1.coupled
        f′ = FusionTree(uncoupled, coupled, isdual, inner, vertices)
        return fusiontreedict(I){F,T}(f′ => coeff)
    else # recursive definition
        N2 = length(f2)
        f2′, f2′′ = split(f2, N2 - 1)
        local newtrees::fusiontreedict(I){F,T}
        for (f, coeff) in insertat(f1, i, f2′′)
            for (f′, coeff′) in insertat(f, i, f2′)
                if @isdefined newtrees
                    coeff′′ = coeff*coeff′
                    newtrees[f′] = get(newtrees, f′, zero(coeff′′)) + coeff′′
                else
                    newtrees = fusiontreedict(I){F,T}(f′ => coeff*coeff′)
                end
            end
        end
        return newtrees
    end
end

"""
    merge(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂}, c::I, μ = nothing)
    -> <:AbstractDict{<:FusionTree{I, N₁+N₂}, <:Number}

Merge two fusion trees together to a linear combination of fusion trees whose uncoupled
sectors are those of `f1` followed by those of `f2`, and where the two coupled sectors of
`f1` and `f2` are further fused to `c`. In case of `FusionStyle(I) == GenericFusion()`,
also a degeneracy label `μ` for the fusion of the coupled sectors of `f1` and `f2` to `c`
needs to be specified.
"""
function merge(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂},
                    c::I, μ = nothing) where {I, N₁, N₂}
    if FusionStyle(I) isa GenericFusion && μ === nothing
        throw(ArgumentError("vertex label for merging required"))
    end
    if !(c in f1.coupled ⊗ f2.coupled)
        throw(SectorMismatch("cannot fuse sectors $(f1.coupled) and $(f2.coupled) to $c"))
    end
    f0 = FusionTree((f1.coupled, f2.coupled), c, (false, false), (), (μ,))
    f, coeff = first(insertat(f0, 1, f1)) # takes fast path, single output
    @assert coeff == one(coeff)
    return insertat(f, N₁+1, f2)
end
function merge(f1::FusionTree{I, 0}, f2::FusionTree{I, 0}, c::I, μ =nothing) where {I}
    c == one(I) ||
        throw(SectorMismatch("cannot fuse sectors $(f1.coupled) and $(f2.coupled) to $c"))
    return fusiontreedict(I)(f1=>Fsymbol(c, c, c, c, c, c)[1,1,1,1])
end

"""
    bendright(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂}) where {I<:Sector, N₁, N₂}

Bend the last uncoupled space of `f1` upward from right hand side to be the last uncoupled
space of new `f2`. That is, map final splitting vertex (a, b)<-c of `f1` to fusion vertex
a<-(c, dual(b)).

Return a dictionary with new tree pairs as keys and coefficients in the linear combination
as values.

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
function bendright(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂}) where {I<:Sector, N₁, N₂}
    @assert N₁ > 0
    c = f1.coupled
    a = N₁ == 1 ? one(I) : (N₁ == 2 ? f1.uncoupled[1] : f1.innerlines[end])
    b = f1.uncoupled[N₁]

    uncoupled1 = Base.front(f1.uncoupled)
    isdual1 = Base.front(f1.isdual)
    inner1 = N₁ > 2 ? Base.front(f1.innerlines) : ()
    vertices1 = N₁ > 1 ? Base.front(f1.vertices) : ()
    f1′ = FusionTree(uncoupled1, a, isdual1, inner1, vertices1)

    uncoupled2 = (f2.uncoupled..., dual(b))
    isdual2 = (f2.isdual..., !(f1.isdual[N₁]))
    inner2 = N₂ > 1 ? (f2.innerlines..., c) : ()

    if FusionStyle(I) isa MultiplicityFreeFusion
        coeff = sqrtdim(c) * isqrtdim(a) * Bsymbol(a, b, c)
        if f1.isdual[N₁]
            coeff *= conj(frobeniusschur(dual(b)))
        end
        vertices2 = N₂ > 0 ? (f2.vertices..., nothing) : ()
        f2′ = FusionTree(uncoupled2, a, isdual2, inner2, vertices2)
        return SingletonDict( (f1′, f2′) => coeff )
    else
        local newtrees
        Bmat = Bsymbol(a, b, c)
        μ = N₁ > 1 ? f1.vertices[end] : 1
        for ν = 1:size(Bmat, 2)
            coeff = sqrtdim(c) * isqrtdim(a) * Bmat[μ,ν]
            iszero(coeff) && continue
            if f1.isdual[N₁]
                coeff *= conj(frobeniusschur(dual(b)))
            end
            vertices2 = N₂ > 0 ? (f2.vertices..., ν) : ()
            f2′ = FusionTree(uncoupled2, a, isdual2, inner2, vertices2)
            if @isdefined newtrees
                push!(newtrees, (f1′, f2′) => coeff)
            else
                newtrees = FusionTreeDict( (f1′, f2′) => coeff )
            end
        end
        return newtrees
    end
end

"""
    bendleft(f1::FusionTree{I}, f2::FusionTree{I}) where I

Bend the last uncoupled space of fusion tree `f2` downward from right hand side to be the
last uncoupled space of the splitting tree constructed from `f1`. That is, map final fusion
vertex c<-(a, b) of `f2` to splitting vertex (c, dual(b))<-a.

Return a dictionary with new tree pairs as keys and coefficients in the linear combination
as values.

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
function bendleft(f1::FusionTree{I}, f2::FusionTree{I}) where I
    return fusiontreedict(I)((f1′, f2′) => conj(coeff) for
                                ((f2′, f1′), coeff) in bendright(f2, f1))
end

"""
    foldright(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂}) where {I<:Sector, N₁, N₂}

Bend the first uncoupled space of the splitting tree `f1` upward from left hand side
to be the first uncoupled space of new fusion tree constructed from `f2`. That is, map first
splitting vertex ``(a, b)<-c`` of `f1` to fusion vertex ``b<-(dual(a), c)``.

Return a dictionary with new tree pairs as keys and coefficients in the linear combination
as values.

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
function foldright(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂}) where {I<:Sector, N₁, N₂}
    @assert N₁ > 0
    if FusionStyle(I) isa UniqueFusion
        a = f1.uncoupled[1]
        isduala = f1.isdual[1]
        factor = sqrtdim(a)
        if !isduala
            factor *= frobeniusschur(a)
        end
        c1 = dual(a)
        c2 = f1.coupled
        c = first(c1 ⊗ c2)
        fl = FusionTree{I}(Base.tail(f1.uncoupled), c, Base.tail(f1.isdual))
        fr = FusionTree{I}((c1, f2.uncoupled...), c, (!isduala, f2.isdual...))
        return fusiontreedict(I)((fl,fr)=>1)
    else
        a = f1.uncoupled[1]
        isduala = f1.isdual[1]
        factor = sqrtdim(a)
        if !isduala
            factor *= frobeniusschur(a)
        end
        c1 = dual(a)
        c2 = f1.coupled
        hasmultiplicities = FusionStyle(a) isa GenericFusion
        local newtrees
        for c in c1 ⊗ c2
            N₁ == 1 && c != one(c) && continue
            for μ in (hasmultiplicities ? (1:Nsymbol(c1, c2, c)) : (nothing,))
                fc = FusionTree((c1, c2), c, (!isduala, false), (), (μ,))
                for (fl′, coeff1) in insertat(fc, 2, f1)
                    N₁ > 1 && fl′.innerlines[1] != one(I) && continue
                    coupled = fl′.coupled
                    uncoupled = Base.tail(Base.tail(fl′.uncoupled))
                    isdual = Base.tail(Base.tail(fl′.isdual))
                    inner = N₁ <= 3 ? () : Base.tail(Base.tail(fl′.innerlines))
                    vertices = N₁ <= 2 ? () : Base.tail(Base.tail(fl′.vertices))
                    fl = FusionTree{I}(uncoupled, coupled, isdual, inner, vertices)
                    for (fr, coeff2) in insertat(fc, 2, f2)
                        coeff = factor * coeff1 * coeff2
                        if (@isdefined newtrees)
                            newtrees[(fl,fr)] = get(newtrees, (fl, fr), zero(coeff)) + coeff
                        else
                            newtrees = fusiontreedict(I)((fl,fr)=>coeff)
                        end
                    end
                end
            end
        end
        return newtrees
    end
end

"""
    foldleft(f1::FusionTree{I}, f2::FusionTree{I}) where I

Bend the first uncoupled space of the fusion tree `f2` downward from left hand side
to be the first uncoupled space of new fusion tree constructed from `f1`. That is, map first
fusion vertex ``c<-(a, b)`` of `f2` to splitting vertex ``(dual(a), c)<-b``.

Return a dictionary with new tree pairs as keys and coefficients in the linear combination
as values.

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
function foldleft(f1::FusionTree{I}, f2::FusionTree{I}) where I
    return fusiontreedict(I)((f1′, f2′) => conj(coeff) for
                                    ((f2′, f1′), coeff) in foldright(f2, f1))
end

"""
    cycleclockwise(f1::FusionTree{I}, f2::FusionTree{I}) where {I<:Sector}

Clockwise cyclic permutation while preserving (N₁, N₂): foldright & bendleft

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
function cycleclockwise(f1::FusionTree{I}, f2::FusionTree{I}) where {I<:Sector}
    local newtrees
    if length(f1) > 0
        for ((f1a, f2a), coeffa) in foldright(f1, f2)
            for ((f1b, f2b), coeffb) in bendleft(f1a, f2a)
                coeff = coeffa * coeffb
                if (@isdefined newtrees)
                    newtrees[(f1b,f2b)] = get(newtrees, (f1b, f2b), zero(coeff)) + coeff
                else
                    newtrees = fusiontreedict(I)((f1b,f2b)=>coeff)
                end
            end
        end
    else
        for ((f1a, f2a), coeffa) in bendleft(f1, f2)
            for ((f1b, f2b), coeffb) in foldright(f1a, f2a)
                coeff = coeffa * coeffb
                if (@isdefined newtrees)
                    newtrees[(f1b,f2b)] = get(newtrees, (f1b, f2b), zero(coeff)) + coeff
                else
                    newtrees = fusiontreedict(I)((f1b,f2b)=>coeff)
                end
            end
        end
    end
    return newtrees
end

"""
    cycleanticlockwise(f1::FusionTree{I}, f2::FusionTree{I}) where {I<:Sector}

Anticlockwise cyclic permutation while preserving (N₁, N₂): foldleft & bendright

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
function cycleanticlockwise(f1::FusionTree{I}, f2::FusionTree{I}) where {I<:Sector}
    local newtrees
    if length(f2) > 0
        for ((f1a, f2a), coeffa) in foldleft(f1, f2)
            for ((f1b, f2b), coeffb) in bendright(f1a, f2a)
                coeff = coeffa * coeffb
                if (@isdefined newtrees)
                    newtrees[(f1b,f2b)] = get(newtrees, (f1b, f2b), zero(coeff)) + coeff
                else
                    newtrees = fusiontreedict(I)((f1b,f2b)=>coeff)
                end
            end
        end
    else
        for ((f1a, f2a), coeffa) in bendright(f1, f2)
            for ((f1b, f2b), coeffb) in foldleft(f1a, f2a)
                coeff = coeffa * coeffb
                if (@isdefined newtrees)
                    newtrees[(f1b,f2b)] = get(newtrees, (f1b, f2b), zero(coeff)) + coeff
                else
                    newtrees = fusiontreedict(I)((f1b,f2b)=>coeff)
                end
            end
        end
    end
    return newtrees
end

"""
    repartition(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂}, N::Int) where {I, N₁, N₂}
    -> <:AbstractDict{Tuple{FusionTree{I, N}, FusionTree{I, N₁+N₂-N}}, <:Number}

Computes new trees and corresponding coefficients obtained from repartitioning the tree by
bending incoming to outgoing sectors (or vice versa) in order to have `N` outgoing sectors,
corresponding to the first `N` sectors out of the list
``(a_1, a_2, …, a_{N_1}, b_{N_2}^*, …, b_{1}^*)`` and `N₁+N₂-N` incoming sectors,
corresponding to the dual of the last `N₁+N₂-N` sectors from the previous list in reverse
order.

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
@inline function repartition(f1::FusionTree{I, N₁},
                        f2::FusionTree{I, N₂},
                        N::Int) where {I<:Sector, N₁, N₂}
    f1.coupled == f2.coupled || throw(SectorMismatch())
    @assert 0 <= N <= N₁+N₂
    return _recursive_repartition(f1, f2, Val(N))
end

function _recursive_repartition(f1::FusionTree{I, N₁},
                                f2::FusionTree{I, N₂},
                                ::Val{N}) where {I<:Sector, N₁, N₂, N}
    # recursive definition is only way to get correct number of loops for
    # GenericFusion, but is too complex for type inference to handle, so we
    # precompute the parameters of the return type
    F1 = fusiontreetype(I, N)
    F2 = fusiontreetype(I, N₁ + N₂ - N)
    coeff = @inbounds Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1,1,1,1]
    T = typeof(coeff)
    if N == N₁
        return fusiontreedict(I){Tuple{F1, F2}, T}( (f1, f2) => coeff)
    else
        local newtrees::fusiontreedict(I){Tuple{F1, F2}, T}
        for ((f1′, f2′), coeff1) in (N < N₁ ? bendright(f1, f2) : bendleft(f1, f2))
            for ((f1′′, f2′′), coeff2) in _recursive_repartition(f1′, f2′, Val(N))
                if (@isdefined newtrees)
                    push!(newtrees, (f1′′, f2′′) => coeff1*coeff2)
                else
                    newtrees =
                        fusiontreedict(I){Tuple{F1, F2}, T}((f1′′, f2′′) => coeff1*coeff2)
                end
            end
        end
        return newtrees
    end
end

const transposecache = LRU{Any, Any}(; maxsize = 10^5)
const usetransposecache = Ref{Bool}(true)

"""
    linearizepermutation(p1::NTuple{N₁, Int}, p2::NTuple{N₂,Int},
                                    n₁::Int, n₂::Int) where {N₁, N₂}

The input fusion tree pair has `n₁` outgoing and `n₂` incoming sectors. The output fusion
tree has `N₁` outgoing corresponding to `p1` and `N₂` incoming sectors corresponding to `p2`.

First use repartition to map the input fusion-splitting tree to a splitting tree. Then, map
the permutations `p1` and `p2` to an equivalent permutation on the splitting tree, which
is returned.

Note that the sequence of the incoming sectors of the fusion tree is reversed after
repartition into the splitting tree.
"""
function linearizepermutation(p1::NTuple{N₁, Int}, p2::NTuple{N₂,Int},
                                n₁::Int, n₂::Int) where {N₁, N₂}
    p1′ = ntuple(Val(N₁)) do n
        p1[n] > n₁ ? n₁+(n₂+n₁+1-p1[n]) : p1[n]
    end
    p2′ = ntuple(Val(N₂)) do n
        p2[N₂+1-n] > n₁ ? n₁+(n₂+n₁+1-p2[N₂+1-n]) : p2[N₂+1-n]
    end
    return (p1′..., p2′...)
end

"""
    iscyclicpermutation(p)

Check whether `p` is a cyclicpermutation, that is `p[i+1] == p[i]+1` mod `N`.
"""
function iscyclicpermutation(p)
    N = length(p)
    @inbounds for i = 1:N
        p[mod1(i+1, N)] == mod1(p[i] + 1, N) || return false
    end
    return true
end

"""
    iscyclicpermutation(v1, v2)

Check whether `v1` is a cyclic permutation of `v2`.
"""
function iscyclicpermutation(v1, v2)
    length(v1) == length(v2) || return false
    return iscyclicpermutation(indexin(v1, v2))
end

"""
    transpose(f1::FusionTree{I}, f2::FusionTree{I},
            p1::NTuple{N₁, Int}, p2::NTuple{N₂, Int}) where {I, N₁, N₂}
    -> <:AbstractDict{Tuple{FusionTree{I, N₁}, FusionTree{I, N₂}}, <:Number}

Computes new trees and corresponding coefficients obtained from repartitioning and cyclic
permutating the fusion-splitting tree such that sectors `p1` become outgoing and sectors
`p2` become incoming. It is required that the linearized permutation is cyclic to avoid
braiding.

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
function Base.transpose(f1::FusionTree{I}, f2::FusionTree{I},
                    p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {I<:Sector, N₁, N₂}
    N = N₁ + N₂
    @assert length(f1) + length(f2) == N
    p = linearizepermutation(p1, p2, length(f1), length(f2))
    @assert iscyclicpermutation(p)
    if usetransposecache[]
        u = one(I)
        T = eltype(Fsymbol(u, u, u, u, u, u))
        F₁ = fusiontreetype(I, N₁)
        F₂ = fusiontreetype(I, N₂)
        D = fusiontreedict(I){Tuple{F₁, F₂}, T}
        return _get_transpose(D, (f1, f2, p1, p2))
    else
        return _transpose((f1, f2, p1, p2))
    end
end

@noinline function _get_transpose(::Type{D}, @nospecialize(key)) where D
    d::D = get!(transposecache, key) do
        _transpose(key)
    end
    return d
end

const TransposeKey{I<:Sector, N₁, N₂} = Tuple{<:FusionTree{I}, <:FusionTree{I},
                                                IndexTuple{N₁}, IndexTuple{N₂}}

function _transpose((f1, f2, p1, p2)::TransposeKey{I,N₁,N₂}) where {I<:Sector, N₁, N₂}
    N = N₁ + N₂
    p = linearizepermutation(p1, p2, length(f1), length(f2))
    newtrees = repartition(f1, f2, N₁)
    length(p) == 0 && return newtrees
    i1 = findfirst(==(1), p)
    @assert i1 !== nothing
    i1 == 1 && return newtrees
    Nhalf = N >> 1
    while 1 < i1 <= Nhalf
        local newtrees′
        for ((f1a, f2a), coeffa) in newtrees
            for ((f1b, f2b), coeffb) in cycleanticlockwise(f1a, f2a)
                coeff = coeffa * coeffb
                if (@isdefined newtrees′)
                    newtrees′[(f1b, f2b)] = get(newtrees′, (f1b, f2b), zero(coeff)) + coeff
                else
                    newtrees′ = fusiontreedict(I)((f1b, f2b) => coeff)
                end
            end
        end
        newtrees = newtrees′
        i1 -= 1
    end
    while Nhalf < i1
        local newtrees′
        for ((f1a, f2a), coeffa) in newtrees
            for ((f1b, f2b), coeffb) in cycleclockwise(f1a, f2a)
                coeff = coeffa * coeffb
                if (@isdefined newtrees′)
                    newtrees′[(f1b, f2b)] = get(newtrees′, (f1b, f2b), zero(coeff)) + coeff
                else
                    newtrees′ = fusiontreedict(I)((f1b, f2b) => coeff)
                end
            end
        end
        newtrees = newtrees′
        i1 = mod1(i1 + 1, N)
    end
    return newtrees
end

"""
    elementary_trace(f::FusionTree{I, N}, i) where {I<:Sector, N}

Trace the ``i``th and ``i+1``th uncoupled sectors of splitting tree `f`. The trace is taken
as the evaluation map.
"""
function elementary_trace(f::FusionTree{I, N}, i) where {I<:Sector, N}
    (N > 1 && 1 <= i <= N) ||
        throw(ArgumentError("Cannot trace outputs i=$i and i+1 out of only $N outputs"))
    i < N || f.coupled == one(I) ||
        throw(ArgumentError("Cannot trace outputs i=$N and 1 of fusion tree that
                                couples to non-trivial sector"))

    u = one(I)
    T = typeof(Fsymbol(u,u,u,u,u,u)[1,1,1,1])
    F = fusiontreetype(I, N-2)
    newtrees = FusionTreeDict{F,T}()

    j = mod1(i+1, N)
    b = f.uncoupled[i]
    b′ = f.uncoupled[j]
    # if trace is zero, return empty dict
    (b == dual(b′) && f.isdual[i] != f.isdual[j]) || return newtrees
    if i < N
        a = i == 1 ? one(I) : (i == 2 ? f.uncoupled[1] : f.innerlines[i-2])
        d = i == N-1 ? f.coupled : f.innerlines[i]
        a == d || return newtrees
        uncoupled′ = TupleLabXD.deleteat(TupleLabXD.deleteat(f.uncoupled, i+1), i)
        isdual′ = TupleLabXD.deleteat(TupleLabXD.deleteat(f.isdual, i+1), i)
        coupled′ = f.coupled
        if N <= 4
            inner′ = ()
        else
            inner′ = i <= 2 ? Base.tail(Base.tail(f.innerlines)) :
                        TupleLabXD.deleteat(TupleLabXD.deleteat(f.innerlines, i-1), i-2)
        end
        if N <= 3
            vertices′ = ()
        else
            vertices′ = i <= 2 ? Base.tail(Base.tail(f.vertices)) :
                        TupleLabXD.deleteat(TupleLabXD.deleteat(f.vertices, i), i-1)
        end
        f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′)
        coeff = sqrtdim(b)
        if i > 1
            c = f.innerlines[i-1]
            if FusionStyle(I) isa MultiplicityFreeFusion
                coeff *= Fsymbol(a, b, dual(b), a, c, one(I))
            else
                μ = f.vertices[i-1]
                ν = f.vertices[i]
                coeff *= Fsymbol(a, b, dual(b), a, c, one(I))[μ, ν, 1, 1]
            end
        end
        if f.isdual[i]
            coeff *= frobeniusschur(b)
        end
        push!(newtrees, f′ => coeff)
        return newtrees
    else # i == N
        if N == 2
            f′ = FusionTree{I}((), one(I), (), (), ())
            coeff = sqrtdim(b)
            if !(f.isdual[N])
                coeff *= conj(frobeniusschur(b))
            end
            push!(newtrees, f′ => coeff)
            return newtrees
        end
        uncoupled_ = Base.front(f.uncoupled)
        inner_ = Base.front(f.innerlines)
        coupled_ = f.innerlines[end]
        @assert coupled_ == dual(b)
        isdual_ = Base.front(f.isdual)
        vertices_ = Base.front(f.vertices)
        f_ = FusionTree(uncoupled_, coupled_, isdual_, inner_, vertices_)
        fs = FusionTree((b,), b, (!f.isdual[1],), (), ())
        for (f_′, coeff) = merge(fs, f_, one(I), 1)
            f_′.innerlines[1] == one(I) || continue
            uncoupled′ = Base.tail(Base.tail(f_′.uncoupled))
            isdual′ = Base.tail(Base.tail(f_′.isdual))
            inner′ = N <= 4 ? () : Base.tail(Base.tail(f_′.innerlines))
            vertices′ = N <= 3 ? () : Base.tail(Base.tail(f_′.vertices))
            f′ = FusionTree(uncoupled′, one(I), isdual′, inner′, vertices′)
            coeff *= sqrtdim(b)
            if !(f.isdual[N])
                coeff *= conj(frobeniusschur(b))
            end
            newtrees[f′] = get(newtrees, f′, zero(coeff)) + coeff
        end
        return newtrees
    end
end

"""
    planar_trace(f::FusionTree{I,N},
                    q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector, N, N₃}

Take the traces between `q1[k]`th and `q2[k]`th sector of the splitting tree `f`, where
`1<=k<=N₃`. All the traces must be planar, i.e., no tracing lines cross with the outgoing
lines nor cross with each other.
"""
function planar_trace(f::FusionTree{I,N},
                        q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector, N, N₃}

    u = one(I)
    T = typeof(Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1])
    F = fusiontreetype(I, N - 2*N₃)
    newtrees = FusionTreeDict{F,T}()
    N₃ === 0 && return push!(newtrees, f=>one(T))

    for (i,j) in zip(q1, q2)
        (f.uncoupled[i] == dual(f.uncoupled[j]) && f.isdual[i] != f.isdual[j]) ||
            return newtrees
    end
    k = 1
    local i, j
    while k <= N₃
        if mod1(q1[k] + 1, N) == q2[k]
            i = q1[k]
            j = q2[k]
            break
        elseif mod1(q2[k] + 1, N) == q1[k]
            i = q2[k]
            j = q1[k]
            break
        else
            k += 1
        end
    end
    k > N₃ && throw(ArgumentError("Not a planar trace"))

    q1′ = let i = i, j = j
        map(l->(l - (l>i) - (l>j)), TupleLabXD.deleteat(q1, k))
    end
    q2′ = let i = i, j = j
        map(l->(l - (l>i) - (l>j)), TupleLabXD.deleteat(q2, k))
    end
    for (f′, coeff′) in elementary_trace(f, i)
        for (f′′, coeff′′) in planar_trace(f′, q1′, q2′)
            coeff = coeff′ * coeff′′
            if !iszero(coeff)
                newtrees[f′′] = get(newtrees, f′′, zero(coeff)) + coeff
            end
        end
    end
    return newtrees
end

"""
    planar_trace(f1::FusionTree{I}, f2::FusionTree{I},
                    p1::IndexTuple{N₁}, p2::IndexTuple{N₂},
                    q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector, N₁, N₂, N₃}

Take the traces between `q1[k]`th and `q2[k]`th sector of the fusion-splitting tree
and the result is transposed according to `p1` and `p2`. We need to make sure that the
traces are planar without any crossing between lines.

# Arguments:
- `f1`: splitting tree
- `f2`: fusion tree
"""
function planar_trace(f1::FusionTree{I}, f2::FusionTree{I},
                    p1::IndexTuple{N₁}, p2::IndexTuple{N₂},
                    q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {I<:Sector, N₁, N₂, N₃}

    N = N₁ + N₂ + 2N₃
    @assert length(f1) + length(f2) == N
    if N₃ == 0
        return transpose(f1, f2, p1, p2)
    end

    linearindex = (ntuple(identity, Val(length(f1)))...,
                    reverse(length(f1) .+ ntuple(identity, Val(length(f2))))...)

    q1′ = TupleLabXD.getindices(linearindex, q1)
    q2′ = TupleLabXD.getindices(linearindex, q2)
    p1′, p2′ = let q′ = (q1′..., q2′...)
        (map(l-> l - count(l .> q′), TupleLabXD.getindices(linearindex, p1)),
            map(l-> l - count(l .> q′), TupleLabXD.getindices(linearindex, p2)))
    end

    u = one(I)
    T = typeof(Fsymbol(u, u, u, u, u, u)[1, 1, 1, 1])
    F₁ = fusiontreetype(I, N₁)
    F₂ = fusiontreetype(I, N₂)
    newtrees = FusionTreeDict{Tuple{F₁,F₂}, T}()
    for ((f1′, f2′), coeff′) in repartition(f1, f2, N)
        for (f1′′, coeff′′) in planar_trace(f1′, q1′, q2′)
            for (f12′′′, coeff′′′) in transpose(f1′′, f2′, p1′, p2′)
                coeff = coeff′ * coeff′′ * coeff′′′
                if !iszero(coeff)
                    newtrees[f12′′′] = get(newtrees, f12′′′, zero(coeff)) + coeff
                end
            end
        end
    end
    return newtrees
end

# BRAIDING MANIPULATIONS:
# =========================================================================================#
"""
    artin_braid(f::FusionTree, i; inv::Bool = false) -> <:AbstractDict{typeof(f), <:Number}

Perform an elementary braid (Artin generator) of neighbouring uncoupled indices `i` and
`i+1` on a splitting tree `f`, and returns the result as a dictionary of output trees and
corresponding coefficients.

The keyword `inv` determines whether index `i` will braid above or below index `i+1`, i.e.
applying `artin_braid(f′, i; inv = true)` to all the outputs `f′` of
`artin_braid(f, i; inv = false)` and collecting the results should yield a single splitting
tree `f` with coefficient `1`. This keyword has no effect if
`BraidingStyle(sectortype(f)) isa SymmetricBraiding`.
"""
function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I<:Sector, N}
    1 <= i < N ||
        throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
    uncoupled = f.uncoupled
    coupled′ = f.coupled
    isdual′ = TupleLabXD.setindex(f.isdual, f.isdual[i], i+1)
    isdual′ = TupleLabXD.setindex(isdual′, f.isdual[i+1], i)
    inner = f.innerlines
    vertices = f.vertices
    u = one(I)
    oneT = one(eltype(Rsymbol(u,u,u))) * one(eltype(Fsymbol(u,u,u,u,u,u)))
    if i == 1
        a, b = uncoupled[1], uncoupled[2]
        c = N > 2 ? inner[1] : coupled′
        uncoupled′ = TupleLabXD.setindex(uncoupled, b, 1)
        uncoupled′ = TupleLabXD.setindex(uncoupled′, a, 2)
        if FusionStyle(I) isa MultiplicityFreeFusion
            R = oftype(oneT, (inv ? conj(Rsymbol(b, a, c)) : Rsymbol(a, b, c)))
            f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices)
            return fusiontreedict(I)(f′ => R)
        else # GenericFusion
            μ = vertices[1]
            Rmat = inv ? Rsymbol(b, a, c)' : Rsymbol(a, b, c)
            local newtrees
            for ν = 1:size(Rmat, 2)
                R = oftype(oneT, Rmat[μ,ν])
                iszero(R) && continue
                vertices′ = TupleLabXD.setindex(vertices, ν, 1)
                f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices′)
                if (@isdefined newtrees)
                    push!(newtrees, f′ => R)
                else
                    newtrees = fusiontreedict(I)(f′ => R)
                end
            end
            return newtrees
        end
    end
    # case i > 1:
    b = uncoupled[i]
    d = uncoupled[i+1]
    a = i == 2 ? uncoupled[1] : inner[i-2]
    c = inner[i-1]
    e = i == N-1 ? coupled′ : inner[i]
    uncoupled′ = TupleLabXD.setindex(uncoupled, d, i)
    uncoupled′ = TupleLabXD.setindex(uncoupled′, b, i+1)
    if FusionStyle(I) isa UniqueFusion
        inner′ = TupleLabXD.setindex(inner, first(a ⊗ d), i-1)
        bd = first(b ⊗ d)
        R = oftype(oneT, inv ? conj(Rsymbol(d, b, bd)) : Rsymbol(b, d, bd))
        f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′)
        return fusiontreedict(I)(f′ => R)
    elseif FusionStyle(I) isa SimpleFusion
        local newtrees
        for c′ in intersect(a ⊗ d, e ⊗ conj(b)) # c′ is f in the figure
            coeff = oftype(oneT, if inv
                    conj(Rsymbol(d, c, e))*conj(Fsymbol(d, a, b, e, c′, c))*Rsymbol(d, a, c′)
                else
                    Rsymbol(c, d, e)*conj(Fsymbol(d, a, b, e, c′, c))*conj(Rsymbol(a, d, c′))
                end)
            iszero(coeff) && continue
            inner′ = TupleLabXD.setindex(inner, c′, i-1)
            f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′)
            if (@isdefined newtrees)
                push!(newtrees, f′ => coeff)
            else
                newtrees = fusiontreedict(I)(f′ => coeff)
            end
        end
        return newtrees
    else # GenericFusion
        local newtrees
        for c′ in intersect(a ⊗ d, e ⊗ conj(b))
            Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e)
            Rmat2 = inv ? Rsymbol(d, a, c′)' : Rsymbol(a, d, c′) # There's still problem in Jutho's codes
            Fmat = Fsymbol(d, a, b, e, c′, c)
            μ = vertices[i-1]
            ν = vertices[i]
            for σ = 1:Nsymbol(a, d, c′)
                for λ = 1:Nsymbol(c′, b, e)
                    coeff = zero(oneT)
                    for ρ = 1:Nsymbol(d, c, e), κ = 1:Nsymbol(d, a, c′)
                        coeff += Rmat1[ν,ρ]*conj(Fmat[κ,λ,μ,ρ])*conj(Rmat2[σ,κ])
                    end
                    iszero(coeff) && continue
                    vertices′ = TupleLabXD.setindex(vertices, σ, i-1)
                    vertices′ = TupleLabXD.setindex(vertices′, λ, i)
                    inner′ = TupleLabXD.setindex(inner, c′, i-1)
                    f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′)
                    if (@isdefined newtrees)
                        push!(newtrees, f′ => coeff)
                    else
                        newtrees = fusiontreedict(I)(f′ => coeff)
                    end
                end
            end
        end
        return newtrees
    end
end

"""
    permutation2swaps(perm) -> Vector{Int}

Change the permutation `perm` into a vector of Int. Each element `V[i]` of the returned
vector `V` represent a swap operation between `V[i]`th and `V[i]+1`th objects. Applying the
swaps in sequence on `1:length(perm)` gives the `perm`.
"""
function permutation2swaps(perm)
    p = collect(perm)
    @assert isperm(p)
    swaps = Vector{Int}()
    N = length(p)
    for k = 1:N-1
        append!(swaps, p[k]-1:-1:k)
        for l = k+1:N
            if p[l] < p[k]
                p[l] += 1
            end
        end
        p[k] = k
    end
    return swaps
end

"""
    braid(f::FusionTree{<:Sector, N}, levels::NTuple{N, Int}, p::NTuple{N, Int})
    -> <:AbstractDict{typeof(t), <:Number}

Perform a braiding of the uncoupled indices of the fusion tree `f` and return the result as
a `<:AbstractDict` of output trees and corresponding coefficients. The braiding is
determined by specifying that the new sector at position `k` corresponds to the sector that
was originally at the position `i = p[k]`, and assigning to every index `i` of the original
fusion tree a distinct level or depth `levels[i]`. This permutation is then decomposed into
elementary swaps between neighbouring indices, where the swaps are applied as braids such
that if `i` and `j` cross, ``τ_{i,j}`` is applied if `levels[i] < levels[j]` and
``τ_{j,i}^{-1}`` if `levels[i] > levels[j]`. This does not allow to encode the most general
braid, but a general braid can be obtained by combining such operations.
"""
function braid(f::FusionTree{I, N},
                levels::NTuple{N, Int},
                p::NTuple{N, Int}) where {I<:Sector, N}
    TupleLabXD.isperm(p) || throw(ArgumentError("not a valid permutation: $p"))
    if FusionStyle(I) isa UniqueFusion && BraidingStyle(I) isa SymmetricBraiding
        coeff = Rsymbol(one(I), one(I), one(I))
        for i = 2:N
            for j = 1:i-1
                if p[j] > p[i]
                    a, b = f.uncoupled[p[j]], f.uncoupled[p[i]]
                    coeff *= Rsymbol(a, b, first(a ⊗ b))
                end
            end
        end
        uncoupled′ = TupleLabXD._permute(f.uncoupled, p)
        coupled′ = f.coupled
        isdual′ = TupleLabXD._permute(f.isdual, p)
        f′ = FusionTree{I}(uncoupled′, coupled′, isdual′)
        return fusiontreedict(I)(f′ => coeff)
    else
        coeff = Rsymbol(one(I), one(I), one(I))[1,1]
        trees = FusionTreeDict(f => coeff)
        newtrees = empty(trees)
        for s in permutation2swaps(p)
            inv = levels[s] > levels[s+1]
            for (f, c) in trees
                for (f′, c′) in artin_braid(f, s; inv = inv)
                    newtrees[f′] = get(newtrees, f′, zero(coeff)) + c*c′
                end
            end
            l = levels[s]
            levels = TupleLabXD.setindex(levels, levels[s+1], s)
            levels = TupleLabXD.setindex(levels, l, s+1)
            trees, newtrees = newtrees, trees
            empty!(newtrees)
        end
        return trees
    end
end

"""
    permute(f::FusionTree, p::NTuple{N, Int}) -> <:AbstractDict{typeof(t), <:Number}

Perform a permutation of the uncoupled indices of the fusion tree `f` and returns the result
as a `<:AbstractDict` of output trees and corresponding coefficients; this requires that
`BraidingStyle(sectortype(f)) isa SymmetricBraiding`.
"""
function permute(f::FusionTree{I, N}, p::NTuple{N, Int}) where {I<:Sector, N}
    @assert BraidingStyle(I) isa SymmetricBraiding
    return braid(f, ntuple(identity, Val(N)), p)
end

# braid double fusion tree
const braidcache = LRU{Any, Any}(; maxsize = 10^5)
const usebraidcache_abelian = Ref{Bool}(false)
const usebraidcache_nonabelian = Ref{Bool}(true)

"""
    braid(f1::FusionTree{I}, f2::FusionTree{I},
            levels1::IndexTuple, levels2::IndexTuple,
            p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {I<:Sector, N₁, N₂}
    -> <:AbstractDict{Tuple{FusionTree{I, N₁}, FusionTree{I, N₂}}, <:Number}

Input is a fusion-splitting tree that describes the fusion of a set of incoming
uncoupled sectors to a set of outgoing uncoupled sectors, represented using the splitting
tree `f1` and fusion tree `f2`, such that the incoming sectors `f2.uncoupled` are fused to
`f1.coupled == f2.coupled` and then to the outgoing sectors `f1.uncoupled`. Compute new
trees and corresponding coefficients obtained from repartitioning and braiding the tree such
that sectors `p1` become outgoing and sectors `p2` become incoming. The uncoupled indices in
splitting tree `f1` and fusion tree `f2` have levels (or depths) `levels1` and `levels2`
respectively, which determines how indices braid. In particular, if `i` and `j` cross,
``τ_{i,j}`` is applied if `levels[i] < levels[j]` and ``τ_{j,i}^{-1}`` if `levels[i] >
levels[j]`. This does not allow to encode the most general braid, but a general braid can
be obtained by combining such operations.
"""
function braid(f1::FusionTree{I}, f2::FusionTree{I},
                levels1::IndexTuple, levels2::IndexTuple,
                p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {I<:Sector, N₁, N₂}
    @assert length(f1) + length(f2) == N₁ + N₂
    @assert length(f1) == length(levels1) && length(f2) == length(levels2)
    @assert TupleLabXD.isperm((p1..., p2...))
    if FusionStyle(f1) isa UniqueFusion &&
        BraidingStyle(f1) isa SymmetricBraiding
        if usebraidcache_abelian[]
            u = one(I)
            T = Int
            F₁ = fusiontreetype(I, N₁)
            F₂ = fusiontreetype(I, N₂)
            D = SingletonDict{Tuple{F₁, F₂}, T}
            return _get_braid(D, (f1, f2, levels1, levels2, p1, p2))
        else
            return _braid((f1, f2, levels1, levels2, p1, p2))
        end
    else
        if usebraidcache_nonabelian[]
            u = one(I)
            T = typeof(sqrtdim(u)*Fsymbol(u, u, u, u, u, u)[1,1,1,1]*Rsymbol(u, u, u)[1,1])
            F₁ = fusiontreetype(I, N₁)
            F₂ = fusiontreetype(I, N₂)
            D = FusionTreeDict{Tuple{F₁, F₂}, T}
            return _get_braid(D, (f1, f2, levels1, levels2, p1, p2))
        else
            return _braid((f1, f2, levels1, levels2, p1, p2))
        end
    end
end

@noinline function _get_braid(::Type{D}, @nospecialize(key)) where D
    d::D = get!(braidcache, key) do
        _braid(key)
    end
    return d
end

const BraidKey{I<:Sector, N₁, N₂} = Tuple{<:FusionTree{I}, <:FusionTree{I},
                                        IndexTuple, IndexTuple,
                                        IndexTuple{N₁}, IndexTuple{N₂}}

function _braid((f1, f2, l1, l2, p1, p2)::BraidKey{I, N₁, N₂}) where {I<:Sector, N₁, N₂}
    p = linearizepermutation(p1, p2, length(f1), length(f2))
    levels = (l1..., reverse(l2)...)
    local newtrees
    for ((f, f0), coeff1) in repartition(f1, f2, N₁ + N₂)
        for (f′, coeff2) in braid(f, levels, p)
            for ((f1′, f2′), coeff3) in repartition(f′, f0, N₁)
                if @isdefined newtrees
                    newtrees[(f1′, f2′)] = get(newtrees, (f1′, f2′), zero(coeff3)) +
                        coeff1*coeff2*coeff3
                else
                    newtrees = fusiontreedict(I)( (f1′, f2′) => coeff1*coeff2*coeff3 )
                end
            end
        end
    end
    return newtrees
end

"""
    permute(f1::FusionTree{I}, f2::FusionTree{I},
            p1::NTuple{N₁, Int}, p2::NTuple{N₂, Int}) where {I, N₁, N₂}
    -> <:AbstractDict{Tuple{FusionTree{I, N₁}, FusionTree{I, N₂}}, <:Number}

Input is a fusion-splitting tree that describes the fusion of a set of incoming uncoupled
sectors to a set of outgoing uncoupled sectors, represented using the individual trees of
outgoing (`f1`) and incoming sectors (`f2`) respectively (with identical coupled sector
`f1.coupled == f2.coupled`). Computes new trees and corresponding coefficients obtained from
repartitioning and permuting the tree such that sectors `p1` become outgoing and sectors
`p2` become incoming. Work only in the case if `BraidingStyle(I) isa SymmetricBraiding`.
"""
function permute(f1::FusionTree{I}, f2::FusionTree{I},
                    p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {I<:Sector, N₁, N₂}
    @assert BraidingStyle(I) isa SymmetricBraiding
    levels1 = ntuple(identity, length(f1))
    levels2 = length(f1) .+ ntuple(identity, length(f2))
    return braid(f1, f2, levels1, levels2, p1, p2)
end


text:
"""
    CodeBits

A bit sequence type with a maximum length of 16. Internally this is represented
by a 32-bit integer, with the upper 16 bits representing the length and the
lower 16 bits representing the bit pattern.
"""
struct CodeBits
    data::UInt32
    CodeBits(data::UInt32=zero(UInt32)) = new(data)
    CodeBits(length::Integer, word::Unsigned) =
        new(UInt32(length) << 0x10 | (word & (0xFFFF >> UInt8(0x10 - length))))
end

Base.length(c::CodeBits) = Int((c.data >> 0x10) & 0x1F)
Base.size(c::CodeBits) = (length(c),)

word(c::CodeBits) = unsafe_trunc(UInt16, c.data & 0xFFFF)
inc(c::CodeBits) = CodeBits(UInt32(c.data + 0x1))
shl(c::CodeBits) = CodeBits(length(c) + 1, c.data << 0x1)
isvalid(c::CodeBits) = c.data <= 0x0010FFFF

Base.getindex(c::CodeBits, i::Integer) =
    c.data & (oneunit(UInt32) << UInt8(length(c) - i)) != zero(UInt32)

function Base.iterate(c::CodeBits, state::UInt8=0x0)
    state >= length(c) && return nothing
    state += 0x1
    (c[state], state)
end

function Base.show(io::IO, c::CodeBits)
    bits = isvalid(c) ? string(word(c), base=2, pad=length(c)) : "?"
    print(io, "code\"", bits, "\"")
end

macro code_str(p)
  p == "" ? CodeBits() : CodeBits(length(p), parse(UInt16, p, base=2))
end

ext: jl


text:
verbose = true
 using Dates
 if verbose println(now()) end
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 #using RvSpectMLBase
 #using EchelleInstruments#, EchelleInstruments.NEID
 using RvSpectML
 if verbose println("# Loading NeidSolarScripts")    end
 using NeidSolarScripts
 println("# Loading other packages 1/2")
 using ArgParse
 #using Markdown
 using PkgVersion
 # Just for version info (sigh)
 import EchelleInstruments, EchelleCCFs, RvSpectMLBase, RvSpectML, SunAsAStar

println("# Parsing arguments...")
   #lsf_width = 3.0e3
   function parse_commandline()
     s = ArgParseSettings( description = "Make dailiy report from daily rvs.")
     #import_settings!(s, s_files_only, args_only=false)
     add_arg_group!(s, "Files", :argg_files)
     @add_arg_table! s begin
         "input"
             help = "Filename for input with daily rvs (csv)"
             arg_type = String
             #default = "daily_rvs.csv"
         "output"
             help = "Filename for output dailiy RVs (toml)"
             arg_type = String
             default = "daily_summary.toml"
#=
         "--template_file"
            help = "Filename with CCF template stored as 'order_ccfs' (jld2) [TODO]"
            arg_type = String
=#
         "--overwrite"
            help = "Specify it's ok to overwrite the output file."
            #default = true
            action = :store_true
      end
      add_arg_group!(s, "Filter rvs for ", :argg_filter_param)
      @add_arg_table! s begin
         "--target"
            help = "Target field"
            arg_type = String
            default = "Sun"
         "--datestr"
            help = "Filenames that contain date string."
            arg_type = String
        "--max_solar_hour_angle"
           help = "Maximum absolute value of solar hour angle."
           arg_type = Float64
           default = 3.12
         "--max_airmass"
            help = "Maximum airmass."
            arg_type = Float64
            default = 2.0
         "--min_expmeter"
            help = "Minimum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 6e4
#=
         "--min_expmeter_factor"
            help = "Minimum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 0.9
         "--max_expmeter_factor"
            help = "Maximum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 1.1
=#
         "--max_expmeter_frac_rms"
            help = "Maximum fractional RMS of exposure meter flux."
            arg_type = Float64
            default = 0.01
#=
         "--min_pyrhelio_factor"
            help = "Minimum pyrheliometer flux relative to model flux."
            arg_type = Float64
            default = 0.9
         "--max_pyrhelio_factor"
            help = "Maximum pyrheliometer flux relative to model flux."
            arg_type = Float64
            default = 1.1
=#
         "--max_pyrhelio_frac_rms"
            help = "Maximum fractional RMS of pyrheliometer meter flux."
            arg_type = Float64
            default = 0.0035
         #=
         "--min_drp_snr"
            help = "Minimum extracted SNR reported by NEID DRP."
            arg_type = Float64
            #default = 0.0
         =#
            #=
         "--start_time"
            help = "Specify daily start time for CCFs to be processed (HH MM)"
            nargs = 2
            arg_type = Int64
            default = [0, 0]
            #default = [18, 30]
        "--stop_time"
           help = "Specify daily stop time for CCFs to be processed (HH MM)"
           nargs = 2
           arg_type = Int64
           #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
           #default = [22, 30]
           default = [23, 59]
           =#
         "--max_num_spectra"
            help = "Maximum number of spectra to process."
            arg_type = Int
            default = 300   # TODO: Increase  after finish testing
      end

     return parse_args(s)
 end
 args = parse_commandline()

if verbose   println("# Loading other packages 2/2")    end
 using CSV, DataFrames, Query, Dates
 using FileIO, MD5 #, JLD2, SHA
 using StatsBase, Statistics, NaNMath
 using TOML, OrderedCollections

 # Filename arguments
 @assert isfile(args["input"]) || islink(args["input"])
 @assert match(r"\.csv$",args["input"]) != nothing
 daily_rvs_filename = args["input"]

 args["overwrite"] = true   # TODO: Comment after done testing
 if isfile(args["output"]) && !args["overwrite"]
    error("# Output file " * args["output"] * " already exists (size = " * string(filesize(args["output"])) * " ).")
 end
 @assert !isfile(args["output"]) || args["overwrite"]
 @assert match(r"\.toml$",args["output"]) != nothing
 daily_report_filename = args["output"]
 touch(daily_report_filename)

 file_hashes = Dict{String,String}()
 start_processing_time = now()

println("# Reading daily RVs from ", daily_rvs_filename)
df_in = CSV.read(daily_rvs_filename,DataFrame)
warning_str = ""

println("# Filtering for usable observations...")

df_use = df_in |>
 #@filter( args["target"] == nothing || _.target == args["target"] ) |>
 @filter( args["datestr"] == nothing || occursin(args["datestr"],_.Filename) ) |>
 #@filter( _.driftfun == "dailymodel0" ) |>
 @filter( args["max_airmass"] == nothing || _.airmass <= args["max_airmass"] ) |>
 @filter( args["max_solar_hour_angle"] == nothing || abs(_.hour_angle) <= args["max_solar_hour_angle"] ) |>
 #@filter( args["start_time"] == nothing || Time(julian2datetime(_.bjd)) >= start_time ) |>
 #@filter( args["stop_time"] == nothing || Time(julian2datetime(_.bjd)) <= stop_time ) |> # TODO for other instruments may need to deal wtih cross end of 24 UTC
 @filter( args["min_expmeter"] == nothing || _.expmeter_mean >= args["min_expmeter"] ) |>
 DataFrame

if hasproperty(df_use,:pyrflux_mean) && hasproperty(df_use,:pyrflux_rms)
   df_use = df_use |>
      @filter( args["max_pyrhelio_frac_rms"] == nothing || _.pyrflux_rms <= args["max_pyrhelio_frac_rms"] * _.pyrflux_mean ) |>
      DataFrame
else
   add_warning_str = "" * "# Warning:  mean_pyroflux and rms_pyroflux not found.  Reverting to expmeter_mean and expmeter_rms."
   warning_str = "" * add_warning_str
   println(add_warning_str)
   df_use = df_use |>
      @filter( args["max_expmeter_frac_rms"] == nothing || _.expmeter_rms <= args["max_expmeter_frac_rms"] * _.expmeter_mean ) |>
      DataFrame
end

df_use = df_use |>
   #@take(args["max_num_spectra"] ) |> @orderby(_.bjd) |>
   @take(args["max_num_spectra"] ) |> @orderby(_.jd_drp) |>
   DataFrame

println("# Found ", size(df_use,1), " files of ",  size(df_use,1), " to use for RVs.")
daily_out = OrderedDict{String,Any}()
#@assert size(df_use,1) >= 1
if size(df_use,1) >= 1
   println("# Creating daily report.")
   daily_out["title"] = "NEID Solar Observations Daily Report"
   if length(warning_str) > 0
      daily_out["warning"] = warning_str
   end

   obs_date = OrderedDict{String,Any}()
      mean_bjd = mean(df_use.jd_drp)
      obs_date["string"] = isnothing(args["datestr"]) ? Date(julian2datetime(mean_bjd)) : args["datestr"]
      obs_date["mean_bjd"] = mean_bjd
      daily_out["obs_date"] = obs_date

   num_rvs = OrderedDict{String,Any}()
     num_rvs["good"] = size(df_use,1)
     num_rvs["usable"] = size(df_in,1)
     daily_out["num_rvs"] = num_rvs

   rv = OrderedDict{String,Any}()
   rv_drp = OrderedDict{String,Any}()
      drp_vel_to_mps = 1000
      rv_drp["mean_rv"] = mean(df_use.rv_drp)  * drp_vel_to_mps
      rv_drp["median_rv"] = median(df_use.rv_drp) * drp_vel_to_mps
      rv_drp["median_σ_rv"] = median(df_use.σrv_drp) * drp_vel_to_mps
      rv_drp["rms_rvs"] = sqrt(var(df_use.rv_drp,corrected=false)) * drp_vel_to_mps
      (times_binned, rvs_binned) = bin_times_and_rvs_max_Δt(times=df_use.jd_drp, rvs=df_use.rv_drp, Δt_threshold=6/(60*24))
      rv_drp["rms_binned_rvs"] = sqrt(var(rvs_binned,corrected=false)) * drp_vel_to_mps

      trimed_rvs = trim(df_use.rv_drp,prop=0.025) .* drp_vel_to_mps
      if length(collect(trimed_rvs)) >= 2
      rv_drp["winsor_mean_rv"] = mean(trimed_rvs)
      rv_drp["winsor_rms_rv"] = sqrt(mean((trimed_rvs.- rv_drp["winsor_mean_rv"]).^2))
      end
      trimed_binned_rvs = trim(rvs_binned,prop=0.025) .* drp_vel_to_mps
      if length(collect(trimed_binned_rvs)) >= 2
      mean_binned_rvs = mean(trimed_binned_rvs)
      rv_drp["winsor_rms_binned_rvs"] = sqrt(mean((trimed_binned_rvs.- mean_binned_rvs).^2))
      end
   rv["drp"] = rv_drp

   if hasproperty(df_use,:rv_template)
   rv_template= OrderedDict{String,Any}()
      rv_template["mean_rv"] = mean(df_use.rv_template)
      rv_template["median_rv"] = median(df_use.rv_template)
      rv_template["median_σ_rv"] = median(df_use.σrv_template)
      rv_template["rms_rvs"] = sqrt(var(df_use.rv_template,corrected=false))
      (times_binned, rvs_binned) = bin_times_and_rvs_max_Δt(times=df_use.jd_drp, rvs=df_use.rv_template, Δt_threshold=6/(60*24))
      rv_template["rms_binned_rvs"] = sqrt(var(rvs_binned,corrected=false))

      trimed_rvs = trim(df_use.rv_template,prop=0.025)
      if length(collect(trimed_rvs)) >= 2
      rv_template["winsor_mean_rv"] = mean(trimed_rvs)
      rv_template["winsor_rms_rv"] = sqrt(mean((trimed_rvs.- rv_template["winsor_mean_rv"]).^2))
      end
      trimed_binned_rvs = trim(rvs_binned,prop=0.025)
      if length(collect(trimed_binned_rvs)) >= 2
      mean_binned_rvs = mean(trimed_binned_rvs)
      rv_template["winsor_rms_binned_rvs"] = sqrt(mean((trimed_binned_rvs.- mean_binned_rvs).^2))
      end
      rv["template"] = rv_template
   end

   if hasproperty(df_use,:rv_gaussian)
   rv_gauss= OrderedDict{String,Any}()
      rv_gauss["mean_rv"] = mean(df_use.rv_gaussian)
      rv_gauss["median_rv"] = median(df_use.rv_gaussian)
      rv_gauss["median_σ_rv"] = median(df_use.σrv_gaussian)
      rv_gauss["rms_rvs"] = sqrt(var(df_use.rv_gaussian,corrected=false))
      (times_binned, rvs_binned) = bin_times_and_rvs_max_Δt(times=df_use.jd_drp, rvs=df_use.rv_gaussian, Δt_threshold=6/(60*24))
      rv_gauss["rms_binned_rvs"] = sqrt(var(rvs_binned,corrected=false))

      trimed_rvs = trim(df_use.rv_gaussian,prop=0.025)
      if length(collect(trimed_rvs)) >= 2
      rv_gauss["winsor_mean_rv"] = mean(trimed_rvs)
      rv_gauss["winsor_rms_rv"] = sqrt(mean((trimed_rvs.- rv_gauss["winsor_mean_rv"]).^2))
      end
      trimed_binned_rvs = trim(rvs_binned,prop=0.025)
      if length(collect(trimed_binned_rvs)) >= 2
      mean_binned_rvs = mean(trimed_binned_rvs)
      rv_gauss["winsor_rms_binned_rvs"] = sqrt(mean((trimed_binned_rvs.- mean_binned_rvs).^2))
      end
      rv["gaussian"] = rv_gauss
   end

   #rv = OrderedDict{String,Any}("drp"=>rv_drp, "ccf_template"=>rv_template, "ccf_gaussian"=>rv_gauss )
   daily_out["rv"] = rv

   if hasproperty(df_use, :CaIIHK) && hasproperty(df_use, :σ_CaIIHK)
      indicator = OrderedDict{String,Any}()
      indicator["CaIIHK"] = mean(df_use.CaIIHK)
      indicator["σ_CaIIHK"] = median(df_use.σ_CaIIHK)
      indicator["Ha06_1"] = mean(df_use.Ha06_1)
      indicator["σ_Ha06_1"] = median(df_use.σ_Ha06_1)
      indicator["Ha06_2"] = mean(df_use.Ha06_2)
      indicator["σ_Ha06_2"] = median(df_use.σ_Ha06_2)
      indicator["Ha16_1"] = mean(df_use.Ha16_1)
      indicator["σ_Ha16_1"] = median(df_use.σ_Ha16_1)
      indicator["Ha16_2"] = mean(df_use.Ha16_2)
      indicator["σ_Ha16_2"] = median(df_use.σ_Ha16_2)
      #indicator["NaI"] = df_use.NaI
      #indicator["σ_NaI"] = df_use.σ_NaI
      daily_out["indicators"] = indicator
   end

   quality = OrderedDict{String,Any}()
      quality["extsnr"] = median(df_use.drpextsnr)
      if hasproperty(df_use, :expmeter_mean) && hasproperty(df_use, :pyrflux_mean)
      quality["median_expmeter_over_pyrheliometer"] = median(df_use.expmeter_mean ./ df_use.pyrflux_mean)
      quality["min_expmeter_over_pyrheliometer"] = minimum(df_use.expmeter_mean ./ df_use.pyrflux_mean)
      quality["max_expmeter_over_pyrheliometer"] = maximum(df_use.expmeter_mean ./ df_use.pyrflux_mean)
      end
   daily_out["quality"] = quality
   environment = OrderedDict{String,Any}()
      if hasproperty(df_use, :envohum)   environment["humidity"] = median(df_use.envohum) end
      if hasproperty(df_use, :envotmp)   environment["temperature"] = median(df_use.envotmp) end
      if hasproperty(df_use, :envwinds)  environment["wind_speed"] = median(df_use.envwinds) end
      if hasproperty(df_use, :envwindd)  environment["wind_direction"] = median(df_use.envwindd) end
   daily_out["environment"] = environment

   report = OrderedDict{String,Any}()
      report["hostname"]=gethostname()
      report["date_run"]=now()
      report["input_file"]=daily_rvs_filename
      report["input_ctime"] = ctime(daily_rvs_filename)
      report["input_md5"] = bytes2hex(open(md5,daily_rvs_filename))
      report_versions = OrderedDict{String,Any}()
      report_versions["NeidSolarScripts"] = string(PkgVersion.Version(NeidSolarScripts))
      report_versions["EchelleCCFs"] = string(PkgVersion.Version(EchelleCCFs))
      report_versions["EchelleInstruments"] = string(PkgVersion.Version(EchelleInstruments))
      report_versions["RvSpectML"] = string(PkgVersion.Version(RvSpectML))
      report_versions["RvSpectMLBase"] = string(PkgVersion.Version(RvSpectMLBase))
      report_versions["SunAsAStar"] = string(PkgVersion.Version(SunAsAStar))
      report["versions"] = report_versions
   daily_out["report"] = report

   println("# Writing daily report.")
   open(daily_report_filename, "w") do io
        TOML.print(io, daily_out) end
else
   touch(daily_report_filename)
end

#=
   report_str = """
# NEID Solar Observations Daily Report $(args["datestr"])
- Usable files: $(size(manifest_use,1)) of $(size(manifest,1)) solar observations
- Daily mean JD: $daily_mean_bjd   $(julian2datetime(daily_mean_bjd))
- Daily mean RV: $daily_mean_rv
- Daily median RV: $daily_median_rv
- Median of σ_RVs: $daily_median_σ_rv
- Daily RMS of RVs: $daily_rms_rvs

---
This report generated on $(gethostname()) at $(now()) based on CCFs calculated at $(input_data["start_processing_time"]).
Input file: $(input_data["daily_ccf_filename"])
Input md5sum: $input_md5
"""
else
   report_str = """
# NEID Solar Observations Daily Report $(args["datestr"])
- Usable files: $(size(manifest_use,1)) of $(size(manifest,1)) solar observations

---
This report generated on $(gethostname()) at $(now()) based on CCFs calculated at $(input_data["start_processing_time"]).
Input file: $(input_data["daily_ccf_filename"])
Input md5sum: $input_md5
"""
end

write(daily_report_filename, report_str)
=#


text:
using BenchmarkTools, Gillespie, StatsBase

function ssa_wrapper(x, x0, F, nu, parms, t_final, n_trial)
  for i in 1:n_trial
    copy!(x, x0)
    ssa(x, F, nu, parms, t_final)
  end

  return nothing
end

function summarize(results, benchmark_set)

  for m in benchmark_set
    file = open("./gillespie/$(m).txt", "a")

    result = results[m]

    println(file, "model: $(m)")
    println(file)
    println(file, "n = $(length(result.times))")
    println(file, "mean:               $(mean(result.times) * 1e-6) ms")
    println(file, "std:                $(std(result.times) * 1e-6) ms")
    println(file, "median:             $(median(result.times) * 1e-6) ms")
    println(file, "mad (normal=true):  $(mad(result.times, normalize=true) * 1e-6) ms")
    println(file, "mad (normal=false): $(mad(result.times, normalize=false) * 1e-6) ms")
    println(file)

    close(file)
  end

  return nothing
end

# initialize test suite
suite = BenchmarkGroup()

# load the following variables:
# - parameters
# - MODELS
# - SEED
include("benchmark_parameters.jl")

# build handy collections
const benchmark_set = intersect(ARGS, MODELS)

# build up the benchmark suite
for m in benchmark_set
  # unpack parameters
  t_final, n_saves, n_trial, n_sample, t_limit = parameters[m]

  # load model
  include(joinpath(m, "gillespie.jl"))

  x = similar(x0)

  suite[m] = @benchmarkable ssa_wrapper($x, $x0, $F, $nu, $parms, $t_final, $n_trial) samples=n_sample seconds=t_limit setup=(srand($SEED))
end

results = run(suite, verbose = true)

# summarize(results, benchmark_set)
BenchmarkTools.save("gillespie.json", results)


text:
function vor( u :: Array{Float64,2}, v :: Array{Float64,2}, m :: MITgcmDatas)
  dudy = dimDiff(u, 2, OnePointR()) /m.dyspacing
  dvdx = dimDiff(v, 1, OnePointR()) /m.dxspacing
  return (dvdx - dudy)
end

function vor( u :: Array{Float64,3}, v :: Array{Float64,3}, m :: MITgcmDatas)
  vv = zeros(size(u))
  for iz = 1: size(u,3)
    vv[:,:,iz] = vor(u[:,:,iz],v[:,:,iz],m)
  end
  return vv
end

#= include("intercept.jl") =#
function MITgcmPV( u :: Array{Float64,3}, v :: Array{Float64,3},
             T :: Array{Float64,3}, m :: MITgcmDatas;
             dzspacing :: Float64 = 2.5, f :: Float64 = 7.29e-5, BbyT :: Float64 = 2.0e-4 * 9.8)
  parIndent = OnePointR()
  dudy = begin
    fdy =  dimDiff(u, 2, OnePointLR()) /(m.dyspacing)
    dimAverage(fdy,1, parIndent)
    #= dimMoves(fdy,2,OnePointR()) =#
  end
  dvdx = begin
    fdx = dimDiff(v,1,OnePointLR())/(m.dxspacing)
    dimAverage(fdx, 2,parIndent)
    #= dimMoves(fdx,1,OnePointR()) =#
  end
  dudz = begin
    fdz = -dimDiff(u,3,OnePointLR())/dzspacing
    dimAverage(fdz,1, parIndent)
    #= dimMoves(fdz,2,OnePointR()) =#
  end
  dvdz = begin
    fdz = -dimDiff(v,3,OnePointLR())/dzspacing
    dimAverage(fdz,2, parIndent)
    #= dimMoves(fdz,1,OnePointR()) =#
  end
  avor = ( -dvdz, dudz, f+(dvdx - dudy))

  dB = begin
    TT = map( x-> BbyT * dimDiff(T,x, OnePointLR()), (1,2,3))
    (TT[1]/m.dxspacing, TT[2]/m.dyspacing, -TT[3]/dzspacing)
  end

  return mapreduce(+, zip(avor, dB))do x
    x[1].*x[2]
  end
end

export vor
export MITgcmPV, MITgcmVor

include("intercept.jl")
function MFPV( u :: Array{Float64,3}, v :: Array{Float64,3},
             T :: Array{Float64,3}, m :: MITgcmDatas;
             dzspacing :: Float64 = 2.5, f :: Float64 = 7.29e-5, BbyT :: Float64 = 2.0e-4 * 9.8,
             whether_smooth_T :: Bool = true
             )
  dudy = begin
    fdy =  dimDiff(u, 2, OnePointLR()) /(m.dyspacing)
    shift = GridMoving(uGrid, GridPosition(mGrid.x, mGrid.y - 0.5),2.0,2.0)
    intercept(fdy, shift)/(m.dyspacing)
  end
  dvdx = begin
    fdx = dimDiff(v,1,OnePointLR())/(m.dxspacing)
    shift = GridMoving(vGrid, GridPosition(mGrid.x - 0.5, mGrid.y),2.0,2.0)
    intercept(fdx, shift)
  end
  dudz = begin
    fdz = -dimDiff(u,3,OnePointLR())/dzspacing
    shift = GridMoving(uGrid, mGrid,2.,2.)
    intercept(fdz,shift)
  end
  dvdz = begin
    fdz = -dimDiff(v,3,OnePointLR())/dzspacing
    shift = GridMoving(vGrid, mGrid,2.,2.)
    intercept(fdz,shift)/dzspacing
  end
  avor = ( -dvdz, dudz, f+(dvdx - dudy))

  dB = begin
    T1 = intercept(T,GridMoving(0.0, 0.0, 2.0,2.0);  keep_when_notmoving=!whether_smooth_T)
    TT = map( x-> BbyT * dimDiff(T1,x, OnePointLR()), (1,2,3))
    (TT[1]/m.dxspacing, TT[2]/m.dyspacing, -TT[3]/dzspacing)
  end

  return mapreduce(+, zip(avor, dB))do x
    x[1].*x[2]
  end
end
export MFPV


text:
function sumlengths(points)
    T = eltype(points[1])
    result = zeros(T, length(points))
    for i=1:length(points)
        i0 = max(i-1,1)
        p1, p2 = points[i0], points[i]
        if !(any(map(isnan, p1)) || any(map(isnan, p2)))
            result[i] = result[i0] + norm(p1-p2)
        else
            result[i] = result[i0]
        end
    end
    result
end

function _default{T<:Point}(position::Union{VecTypes{T}, MatTypes{T}}, s::style"lines", data::Dict)
    pv = value(position)
    p_vec = if isa(position, GPUArray)
        position
    else
        const_lift(position) do p
            pv = vec(p)
            if length(pv) < 4 # geometryshader doesn't work with less then 4
                return [pv..., fill(T(NaN), 4-length(pv))...]
            else
                return pv
            end
        end
    end
    _startend = const_lift(p_vec) do vec
        l = length(vec)
        map(1:l) do i
            (i == 1 || isnan(vec[max(i-1, 1)])) && return Float32(0) # start
            (i == l || isnan(vec[min(i+1, l)])) && return Float32(1) # end
            Float32(2) # segment
        end
    end
    @gen_defaults! data begin
        dims::Vec{2, Int32} = ndims(pv) == 1 ? (length(pv), 1) : size(pv)
        vertex              = p_vec  => GLBuffer
        color               = default(RGBA, s, 1) => GLBuffer
        thickness::Float32  = 2f0
        pattern             = nothing
        preferred_camera    = :orthographic_pixel
        boundingbox         = GLBoundingBox(to_cpu_mem(value(p_vec)))
        indices             = const_lift(length, p_vec) => to_indices
        shader              = GLVisualizeShader("fragment_output.frag", "util.vert", "lines.vert", "lines.geom", "lines.frag")
        gl_primitive        = GL_LINE_STRIP_ADJACENCY
        startend            = _startend => GLBuffer
    end
    if pattern != nothing
        @gen_defaults! data begin
            pattern_length = last(pattern)
            lastlen   = const_lift(sumlengths, p_vec) => GLBuffer
            maxlength = const_lift(last, lastlen)
        end
    end
    data
end

to_points{T}(x::Vector{LineSegment{T}}) = reinterpret(T, x, (length(x)*2,))

_default{T <: Point}(positions::VecTypes{LineSegment{T}}, s::Style, data::Dict) =
    _default(const_lift(to_points, positions), style"linesegment"(), data)

function _default{T <: Point}(positions::VecTypes{T}, s::style"linesegment", data::Dict)
    @gen_defaults! data begin
        vertex              = positions           => GLBuffer
        color               = default(RGBA, s, 1) => GLBuffer
        thickness::Float32  = 2f0                 => GLBuffer
        shape               = RECTANGLE
        pattern             = nothing
        indices             = const_lift(length, positions) => to_indices
        preferred_camera    = :orthographic_pixel
        boundingbox         = GLBoundingBox(to_cpu_mem(value(positions)))
        shader              = GLVisualizeShader("fragment_output.frag", "util.vert", "line_segment.vert", "line_segment.geom", "lines.frag")
        gl_primitive        = GL_LINES
    end
end

function _default{T <: AbstractFloat}(positions::Vector{T}, range::Range, s::style"lines", data::Dict)
    length(positions) != length(range) && throw(
        DimensionMismatsch("length of $(typeof(positions)) $(length(positions)) and $(typeof(range)) $(length(range)) must match")
    )
    _default(points2f0(positions, range), s, data)
end

function _default{G<:GeometryPrimitive{3}}(
        geometry::TOrSignal{G}, s::style"lines", data::Dict
    )
    wireframe(geometry, data)
end
function _default(
        geometry::TOrSignal{GLNormalMesh}, s::style"lines", data::Dict
    )
    wireframe(geometry, data)
end
function wireframe(
        geometry, data::Dict
    )
    points = const_lift(geometry) do g
        decompose(Point3f0, g) # get the point representation of the geometry
    end
    # Get line index representation
    indices = decompose(Face{2, GLuint, -1}, value(geometry))
    data[:indices] = reinterpret(GLuint, indices)
    _default(points, style"linesegment"(), data)
end

immutable GridPreRender end

@compat function (::GridPreRender)()
    glEnable(GL_CULL_FACE)
    glCullFace(GL_BACK)
end

function _default{T<:AABB}(c::TOrSignal{T}, ::Style{:grid}, data)
    @gen_defaults! data begin
        primitive::GLPlainMesh = c
        bg_color = RGBA{Float32}(0.99,0.99,0.99,1)
        grid_color = RGBA{Float32}(0.8,0.8,0.8,1)
        grid_thickness = Vec3f0(0.999)
        gridsteps = Vec3f0(5)
        shader = GLVisualizeShader("fragment_output.frag", "grid.vert", "grid.frag")
        boundingbox = c
        prerender = GridPreRender()
        postrender = () -> glDisable(GL_CULL_FACE);
    end
end


text:
using SimPN
using Plots
plotly()

t, u = SimPN.use_integrator(dt=0.001)
plot(t, u)

max_stars_repo_name: praneethnamburi/blender-ScriptViz

max_issues_repo_name: praneethnamburi/blender-ScriptViz

max_forks_repo_name: praneethnamburi/blender-ScriptViz


text:
using Plots,DifferentialEquations
Plots.PlotlyJSBackend()
#=
Objective: This script intends to reproduce the Fast Offshore Wells Model.
Paper: Fast Offshore Wells Model (FOWM): A practical dynamic model for multiphase oil production systems in deepwater and ultra-deepwater scenarios
Authors: Fabio C.Diehl Thiago K.Anzai Cristina S.Almeida Oscar F.Von Meien Saul S.Neto Vinicius R.Rosa Mario C.M.M.Campos Filipe Reolon Giovani Gerevini Cassiano Ranza MarceloFarenzena Jorge O.Trierweiler
=#
function fowm(du,u,p,t)
    #=
    States:
        01 - m_ga:      Gas Annular
        02 - m_gt:      Gas Tubing
        03 - m_lt:      Liq Tubing
        04 - m_gb:      Gas Bubble
        05 - m_gr:      Gas Flowline
        06 - m_lr:      Liq Flowline
    =#
    m_ga,m_gt,m_lt,m_gb,m_gr,m_lr = u
    #=
    Parameters:
    - Operation Parameter:
        01 - z: Choke Valve Opening Fraction
    - Well Parameters
        01 - ρ_l:       Liquid Density
        02 - P_r:       Reservoir's Pressure
        03 - P_s:       After Choke Pressure
        04 - α_gw:      Gas Mass Fraction at Reservoir's Condition
        05 - ρ_mres:    Reservoir's Density
        06 - M:         Gas Molecular Weight
        07 - T:         Average Temperature
        08 - L_r:       Riser Length
        09 - L_fl:      Flowline Length
        10 - L_t:       Tubing Length
        11 - L_a:       Annular Length
        12 - H_t:       Distance Xmas Tree - Bottom Hole
        13 - H_pdg:     Distance Xmas Tree - PDG Point
        14 - H_vgl:     Distance Xmas Tree - Gas Lift Injection Point
        15 - D_ss:      Subsea Pipeline Diameter
        16 - D_t:       Tubing Diameter
        17 - D_a:       Annular Diameter
        18 - θ:         Annular Diameter
    - Estimated Parameters in Diehl et al.
        01 - m_still:   Minimum Liquid Mass at Subsea Pipeline
        02 - C_g:       Virtual Valve Flow Constant
        03 - C_out:     Choke Valve Flow Constant
        04 - V_eb:      Bubble Volume
        05 - E:         Mass Fraction of Gas Bypassing Virtual Valve
        06 - K_w:       Xmas Tree Flow Constant
        07 - K_a:       Annular-Tubing Flow Constant
        08 - K_r:       Reservoir's Flow Constant
        09 - ωᵤ:        Bubble Location (Assistant Parameter)
    =#
    z,ρ_l,P_r,P_s,α_gw,ρ_mres,M,T,L_r,L_fl,L_t,L_a,H_t,
    H_pdg,H_vgl,D_ss,D_t,D_a,θ,m_still,C_g,C_out,V_eb,E,K_w,K_a,K_r,ωᵤ = p

    # Universal Gas Constant
    R      = 8314 # J/kmolK
    g      = 9.81
    #=
    First Rank Variables:
        01 - V_t:       Tubing Volume
        02 - V_ss:      Bubble Downstream Volume
        03 - V_a:       Annular Volume
        04 - A_ss:      Subsea Pipeline Cross Sectional Area
        05 - V_gt:      Tubing Gas Volume
        06 - ρ_mt:      Mixture Density
        07 - ρ_gt:      Gas Density
    =#
    V_t    = π*D_t^2/4*L_t
    V_ss   = π*D_ss^2/4*L_r + π*D_ss^2/4*L_fl
    V_a    = π*D_a^2/4*L_a
    A_ss   = π*D_ss^2/4
    V_gt   = V_t - m_lt/ρ_l
    ρ_mt   = (m_gt+m_lt)/V_t
    ρ_gt   = m_gt/V_gt
    #=
    Second Rank Variables:
        01 - P_ai:      Annular Gas Injection Point Pressure
        02 - P_tt:      Tubing's Topside Pressure
        03 - P_tb:      Tubing Gas Injection Point Pressure
        04 - P_pdg:     Pressure at PDG Point
        05 - P_bh:      Pressure at Bottom Hole
        06 - P_rt:      Pressure at the Riser's Topside
        07 - P_rb:      Pressure before the Bubble
        08 - P_eb:      Bubble Pressure
        09 - ρ_ai:      Gas Density at Annular
        10 - α_gt:      Gas Mass Fraction at Tubing
        11 - α_gr:      Gas Mass Fraction at Subsea Pipeline
        12 - α_lr:      Liquid Mass Fraction at Subsea Pipeline

    =#
    P_ai   = m_ga*(R*T/V_a/M+g*L_a/V_a)
    P_tt   = ρ_gt*R*T/M
    P_tb   = P_tt + ρ_mt*g*H_vgl
    P_pdg  = P_tb +ρ_mres*g*(H_pdg-H_vgl)
    P_bh   = P_pdg + ρ_mres*g*(H_t-H_pdg)
    P_rt   = m_gr*R*T/M/(ωᵤ*V_ss-(m_lr+m_still)/ρ_l)
    P_rb   = P_rt + (m_lr+m_still)*g*sin(θ)/A_ss
    P_eb   = m_gb*R*T/M/V_eb
    ρ_ai   = M*P_ai/R/T
    α_gt   = m_gt/(m_gt+m_lt)
    α_gr   = m_gr/(m_gr+m_lr)
    α_lr   = 1-α_gr
    #=
    Mass Flow Rates:
        01 - W_iv:      Gas Annular to tubing
        02 - W_r:       Reservoir to Bottom Hole
        03 - W_whg:     Gas at Xmas Tree
        04 - W_whl:     Liquid at Xmas Tree
        05 - W_gc:      Gas Annular
        06 - W_g:       Gas at Virtual Valve
        07 - W_gout:    Gas through Choke
        08 - W_lout:    Liquid through Choke
    =#
    W_iv   = K_a*√(ρ_ai*max(0,P_ai-P_tb))
    W_r    = K_r*(1 - 0.2*P_bh/P_r-(0.8*P_bh/P_r)^2)
    W_whg  = K_w*α_gt*√(ρ_l*max(0,P_tt-P_rb))
    W_whl  = K_w*(1-α_gt)*√(ρ_l*max(0,P_tt-P_rb))
    W_gc   = W_iv
    W_g    = max(0,C_g*(P_eb-P_rb))
    W_gout = α_gr*C_out*z*√(ρ_l*max(0,P_rt-P_s))
    W_lout = α_lr*C_out*z*√(ρ_l*max(0,P_rt-P_s))
    # States' Derivative
    du[1]  = δm_ga = W_gc - W_iv
    du[2]  = δm_gt = W_r*α_gw+W_iv-W_whg
    du[3]  = δm_lt = W_r*(1-α_gw) - W_whl
    du[4]  = δm_gb = (1-E)*W_whg - W_g
    du[5]  = δm_gr = E*W_whg + W_g - W_gout
    du[6]  = δm_lr = W_whl - W_lout
end

function simulate_fowm(u₀,tf,z,Δt=[])
    #=
    Args:
        - u₀: Initial Condition
        - tf: Final Simulation Time (hours)
        - z:  Choke Valve Opening Fraction
        - Δt: Time Step to Save Results (hours)
    =#
    tspan = (0,tf*3600)
    p = [z,ρ_l,P_r,P_s,α_gw,ρ_mres,M,T,L_r,L_fl,L_t,L_a,H_t,H_pdg,H_vgl,D_ss,D_t,D_a,θ,m_still,C_g,C_out,V_eb,E,K_w,K_a,K_r,ωᵤ]
    prob = ODEProblem(fowm,u₀,tspan,p)
    sol = solve(prob,Tsit5(),saveat = Δt)
    return sol
end

# Parameters for Simulation (Diehl et al. Table 2 Well A and Table 3 Real Well A)
ρ_l         =   900.0
P_r         =   225e5
P_s         =   10e5
α_gw        =   0.0188
ρ_mres      =   892.0
M           =   18.0
T           =   298.0
L_r         =   1569.0
L_fl        =   2928.0
L_t         =   1639.0
L_a         =   1118.0
H_t         =   1279.0
H_pdg       =   1117.0
H_vgl       =   916.0
D_ss        =   0.15
D_t         =   0.15
D_a         =   0.14
θ           =   π*60/180
m_still     =   6.222e1
C_g         =   1.137e-3
C_out       =   2.039e-3
V_eb        =   6.098e1
E           =   1.545e-1
K_w         =   6.876e-4
K_a         =   2.293e-5
K_r         =   1.269e2
ωᵤ          =   2.780e0

# Initial Conditions
x₀ = repeat([10000.0],6)
# Simulation Time and Time-step in Hours
tf = 30.0
Δt = 1/60*15*3600
# Choke Opening Fraction
z = 0.5
# Simulating
sol = simulate_fowm(x₀,tf,z,Δt)

m_ga = sol[1,:]
m_gt = sol[2,:]
m_lt = sol[3,:]
m_gb = sol[4,:]
m_gr = sol[5,:]
m_lr = sol[6,:]

V_t    = π*D_t^2/4*L_t
V_ss   = π*D_ss^2/4*L_r + π*D_ss^2/4*L_fl
V_a    = π*D_a^2/4*L_a
A_ss   = π*D_ss^2/4
V_gt   = V_t .- m_lt./ρ_l
ρ_mt   = (m_gt+m_lt)/V_t
ρ_gt   = m_gt./V_gt

R      = 8314 # J/kmolK
g      = 9.81

P_ai   =  m_ga.*(R.*T./V_a./M+g.*L_a./V_a)
P_tt   =  ρ_gt.*R.*T./M
P_tb   =  P_tt .+ ρ_mt.*g.*H_vgl
P_pdg  =  P_tb .+ρ_mres.*g.*(H_pdg.-H_vgl)
P_bh   =  P_pdg .+ ρ_mres.*g.*(H_t-H_pdg)
P_rt   =  m_gr.*R.*T./M./(ωᵤ.*V_ss.-(m_lr.+m_still)./ρ_l)
P_rb   =  P_rt .+ (m_lr.+m_still).*g.*sin.(θ)./A_ss
P_eb   =  m_gb.*R.*T./M./V_eb
ρ_ai   =  M.*P_ai./R./T
α_gt   =  m_gt./(m_gt.+m_lt)
α_gr   =  m_gr./(m_gr.+m_lr)
α_lr   =  1 .-α_gr

t = sol.t

plotly()
plot(t,P_pdg)


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

#ToDo: Move functionality of reshape_variate and reshape_variates to package ValueShapes.

reshape_variate(shape::Missing, v::Any) = throw(ArgumentError("Can't reshape variate without shape information"))

function reshape_variate(shape::AbstractValueShape, v::Any)
    v_shape = valshape(v)
    if !(v_shape <= shape)
        throw(ArgumentError("Shape of variate doesn't match target variate shape, with variate of type $(typeof(v)) and expected shape $(shape)"))
    end
    v
end

function reshape_variate(shape::ArrayShape{<:Real,1}, v::Any)
    unshaped_v = unshaped(v)::AbstractVector{<:Real}
    reshape_variate(shape, unshaped_v)
end

function reshape_variate(shape::AbstractValueShape, v::AbstractVector{<:Real})
    _reshape_realvec(shape, v)
end

function reshape_variate(shape::ArrayShape{<:Real,1}, v::AbstractVector{<:Real})
    _reshape_realvec(shape, v)
end

function _reshape_realvec(shape::AbstractValueShape, v::AbstractVector{<:Real})
    ndof = length(eachindex(v))
    ndof_expected = totalndof(shape)
    if ndof != ndof_expected
        throw(ArgumentError("Invalid length ($ndof) of variate, target shape  $(shape) has $ndof_expected degrees of freedom"))
    end
    strip_shapedasnt(shape(v))
end

reshape_variates(shape::Missing, vs::AbstractVector{<:Any}) = v

function reshape_variates(shape::AbstractValueShape, vs::AbstractVector{<:Any})
    v_elshape = elshape(vs)
    if !(v_elshape <= shape)
        throw(ArgumentError("Shape of variates doesn't match target variate shape, with variates of type $(eltype(vs)) and expected element shape $(shape)"))
    end
    vs
end

function reshape_variates(shape::ArrayShape{<:Real,1}, vs::AbstractVector{<:Any})
    unshaped_vs = unshaped.(vs)::AbstractVector{<:AbstractVector{<:Real}}
    reshape_variates(shape, unshaped_vs)
end

function reshape_variates(shape::AbstractValueShape, vs::AbstractVector{<:AbstractVector{<:Real}})
    _reshape_realvecs(shape, vs)
end

function reshape_variates(shape::ArrayShape{<:Real,1}, vs::AbstractVector{<:AbstractVector{<:Real}})
    _reshape_realvecs(shape, vs)
end

function _reshape_realvecs(shape::AbstractValueShape, vs::AbstractVector{<:AbstractVector{<:Real}})
    ndof = first(innersize(vs))
    ndof_expected = totalndof(shape)
    if ndof != ndof_expected
        throw(ArgumentError("Invalid length ($ndof) of variates, target shape $(shape) has $ndof_expected degrees of freedom"))
    end
    shape.(vs)
end

function check_variate end

function ChainRulesCore.rrule(::typeof(check_variate), shape::Any, v::Any)
    result = check_variate(shape, v)
    _check_variate_pullback(ΔΩ) = (NoTangent(), NoTangent(), ZeroTangent())
    return result, _check_variate_pullback
end

check_variate(shape::ScalarShape{Real}, v::Real) = nothing

function check_variate(shape::ArrayShape{<:Real,N}, v::AbstractArray{<:Real,N}) where N
    ndof = length(eachindex(v))
    ndof_expected = totalndof(shape)
    if ndof != ndof_expected
        throw(ArgumentError("Invalid length ($ndof) of variate, target shape  $(shape) has $ndof_expected degrees of freedom"))
    end
    nothing
end

function check_variate(shape::NamedTupleShape, v::ShapedAsNT)
    if !(valshape(v) <= shape)
        throw(ArgumentError("Shape of variate incompatible with target variate shape, with variate of type $(typeof(v)) and expected shape $(shape)"))
    end
    nothing
end

function check_variate(shape::NamedTupleShape{names}, v::NamedTuple{names}) where names
    nothing
end

function check_variate(shape::Any, v::Any)
    throw(ArgumentError("Shape of variate incompatible with target variate shape, with variate of type $(typeof(v)) and expected shape $(shape)"))
end

function check_variate(shape::Missing, v::Any)
    throw(ArgumentError("Cannot evaluate without value shape information"))
end


text:
function constructSimModelMesh(model::SimModel3D{T},grid_size) where T
  origin = model.cell_origins[1]
  a_vec = model.cell_origins[end,1,1]-origin+Point3D{T}(model.prim_vec[1,:])
  b_vec = model.cell_origins[1,end,1]-origin+Point3D{T}(model.prim_vec[2,:])
  c_vec = model.cell_origins[1,1,end]-origin+Point3D{T}(model.prim_vec[3,:])
  points = Array{T,3}(grid_size,grid_size,grid_size)
  for (ia,a) in enumerate(linspace(0,1,grid_size)),(ib,b) in enumerate(linspace(0,1,grid_size)),(ic,c) in enumerate(linspace(0,1,grid_size))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    w = zero(Complex{T})
    for wfc in model.wfcs
      w+=wave_equation(point,wfc)
    end
    points[ia,ib,ic] = real(w)
  end
  return points
end

function calculateDipoleMesh(model::SimModel3D{T},cs,grid_size) where T
  origin = model.cell_origins[1]
  a_vec = model.cell_origins[end,1,1]-origin+Point3D{T}(model.prim_vec[1,:])
  b_vec = model.cell_origins[1,end,1]-origin+Point3D{T}(model.prim_vec[2,:])
  c_vec = model.cell_origins[1,1,end]-origin+Point3D{T}(model.prim_vec[3,:])
  dip_mesh = Array{Tuple{Point3D{T},Point3D{T}},3}(grid_size,grid_size,grid_size)
  for (ia,a) in enumerate(linspace(0,1,grid_size)),(ib,b) in enumerate(linspace(0,1,grid_size)),(ic,c) in enumerate(linspace(0,1,grid_size))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    dip = Point3D{T}(0.0)
    for (wfc1,c1) in zip(model.wfcs,cs)
      for (wfc2,c2) in zip(model.wfcs,cs)
        dip += conj(c1*wave_equation(point,wfc1))*c2*wave_equation(point,wfc2)*point
      end
    end
    dip_mesh[ia,ib,ic] = (point,dip)
  end
  return dip_mesh
end

function calculateTotDipole(model::SimModel3D{T},cs,grid_size) where T
  origin = model.cell_origins[1]
  a_vec = model.cell_origins[end,1,1]-origin+Point3D{T}(model.prim_vec[1,:])
  b_vec = model.cell_origins[1,end,1]-origin+Point3D{T}(model.prim_vec[2,:])
  c_vec = model.cell_origins[1,1,end]-origin+Point3D{T}(model.prim_vec[3,:])
  dip_x = zero(typeof(cs[1]))
  dip_y = zero(typeof(cs[1]))
  dip_z = zero(typeof(cs[1]))
  n1 = zero(T)
  n2 = zero(T)
  for (ia,a) in enumerate(linspace(0,1,grid_size)),(ib,b) in enumerate(linspace(0,1,grid_size)),(ic,c) in enumerate(linspace(0,1,grid_size))
    point = origin+(a*a_vec+b*b_vec+c*c_vec)
    for i1=1:length(model.wfcs)
      wfc1,c1 = model.wfcs[i1],cs[i1]
      w1 = wave_equation(point,wfc1)
      n1 += norm(w1)^2
      for i2=1:length(model.wfcs)
        wfc2,c2 = model.wfcs[i2],cs[i2]
        w2 = wave_equation(point,wfc2)
        n2+=norm(w2)^2
        fac = conj(c1*w1)*c2*w2
        dip_x += fac*point.x
        dip_y += fac*point.y
        dip_z += fac*point.z
      end
    end
  end
  n = sqrt(n1*n2)
  return Point3D{T}(real(dip_x)/n,real(dip_y)/n,real(dip_z)/n)
end


text:
"""
    PDDL.Arrays

Extends PDDL with array-valued fluents. Register by calling `PDDL.Arrays.register!`.
Attach to a specific `domain` by calling `PDDL.Arrays.attach!(domain)`.
"""
module Arrays

using ..PDDL
import ..PDDL: defaultval, valterm

# Array constructors
new_array(val, dims...) = fill!(Array{Any}(undef, dims), val)
new_matrix(val, h, w) = fill!(Matrix{Any}(undef, h, w), val)
new_vector(val, n) = fill!(Vector{Any}(undef, n), val)
new_bit_array(val::Bool, dims...) = val ? trues(dims) : falses(dims)
new_bit_matrix(val::Bool, h, w) = val ? trues(h, w) : falses(h, w)
new_bit_vector(val::Bool, n) = val ? trues(n) : falses(n)
vec(xs...) = collect(Any, xs)
bit_vec(xs...) = BitVector(collect(xs))
mat(vs::Vector{Any}...) = Matrix{Any}(hcat(vs...))
bit_mat(vs::BitVector...) = hcat(vs...)::BitMatrix
# Array accessors
get_index(a::AbstractArray{T,N}, idxs::Vararg{Int,N}) where {T,N} =
    getindex(a, idxs...)
set_index(a::AbstractArray{T,N}, val, idxs::Vararg{Int,N}) where {T,N} =
    setindex!(copy(a), val, idxs...)
# Array dimensions
length(a::AbstractArray) = Base.length(a)
height(m::AbstractMatrix) = size(m, 1)
width(m::AbstractMatrix) = size(m, 2)
ndims(a::AbstractArray) = Base.ndims(a)
# Push and pop
push(v::AbstractVector, x) = push!(copy(v), x)
pop(v::AbstractVector) = (v = copy(v); pop!(v); v)
# Transformations
_transpose(v::AbstractVector) = permutedims(v)
_transpose(m::AbstractMatrix) = permutedims(m)

defaultval(::Val{:array}) = Array{Any}(undef, ())
defaultval(::Val{:vector}) = Vector{Any}(undef, 0)
defaultval(::Val{:matrix}) = Matrix{Any}(undef, 0, 0)
defaultval(::Val{Symbol("bit-array")}) = falses()
defaultval(::Val{Symbol("bit-vector")}) = falses(0)
defaultval(::Val{Symbol("bit-matrix")}) = falses(0, 0)

valterm(v::AbstractVector) = Compound(:vec, Const.(v))
valterm(v::BitVector) = Compound(Symbol("bit-vec"), Const.(Int.(v)))
valterm(m::AbstractMatrix) =
    Compound(:transpose, [Compound(:mat, valterm.(eachrow(m)))])
valterm(m::BitMatrix) =
    Compound(:transpose,
        [Compound(Symbol("bit-mat"), valterm.(BitVector.(eachrow(m))))])

const DATATYPES = Dict(
    "array" => Array{Any},
    "vector" => Vector{Any},
    "matrix" => Matrix{Any},
    "bit-array" => BitArray,
    "bit-vector" => BitVector,
    "bit-matrix" => BitMatrix
)

const FUNCTIONS = Dict(
    # Array constructors
    "new-array" => new_array,
    "new-vector" => new_vector,
    "new-matrix" => new_matrix,
    "new-bit-array" => new_bit_array,
    "new-bit-vector" => new_bit_vector,
    "new-bit-matrix" => new_bit_matrix,
    "vec" => vec,
    "bit-vec" => bit_vec,
    "mat" => mat,
    "bit-mat" => bit_mat,
    # Array accessors
    "get-index" => get_index,
    "set-index" => set_index,
    # Array dimensions
    "length" => length,
    "height" => height,
    "width" => width,
    "ndims" => ndims,
    # Push and pop
    "push" => push,
    "pop" => pop,
    # Transformations
    "transpose" => _transpose
)

function register!()
    for (name, ty) in DATATYPES
        PDDL.register!(:datatype, name, ty)
    end
    for (name, f) in FUNCTIONS
        PDDL.register!(:function, name, f)
    end
    return nothing
end

function deregister!()
    for (name, ty) in DATATYPES
        PDDL.deregister!(:datatype, name)
    end
    for (name, f) in FUNCTIONS
        PDDL.deregister!(:function, name)
    end
    return nothing
end

function attach!(domain::GenericDomain)
    for (name, ty) in DATATYPES
        PDDL.attach!(domain, :datatype, name, ty)
    end
    for (name, f) in FUNCTIONS
        PDDL.attach!(domain, :function, name, f)
    end
    return nothing
end

end


text:
names = [
    "test_phasespace.jl",
    "test_spectralanalysis.jl",

    "test_timeevolution_schroedinger.jl",
    "test_timeevolution_master.jl",
    "test_timeevolution_mcwf.jl",
    "test_timeevolution_bloch_redfield.jl",

    "test_timeevolution_twolevel.jl",
    "test_timeevolution_pumpedcavity.jl",

    "test_steadystate.jl",
    "test_timecorrelations.jl",

    "test_semiclassical.jl",

    "test_stochastic_definitions.jl",
    "test_stochastic_schroedinger.jl",
    "test_stochastic_master.jl",
    "test_stochastic_semiclassical.jl",

    "test_timeevolution_abstractdata.jl"
]

detected_tests = filter(
    name->startswith(name, "test_") && endswith(name, ".jl"),
    readdir("."))

unused_tests = setdiff(detected_tests, names)
if length(unused_tests) != 0
    error("The following tests are not used:\n", join(unused_tests, "\n"))
end

unavailable_tests = setdiff(names, detected_tests)
if length(unavailable_tests) != 0
    error("The following tests could not be found:\n", join(unavailable_tests, "\n"))
end

for name=names
    if startswith(name, "test_") && endswith(name, ".jl")
        include(name)
    end
end


text:
"""
    Periodic <: Kernel

TODO
"""
struct Periodic <: Kernel
    # TODO
end


text:
"""
    SEq1

Solutions of Equations in One Variable
here are functions to find root:
1. Bisection
2. Fixed-Point Iteration
3. Newton Method
4. The Secant Method
5. The False Position Method
6.  Müller’s Method
"""
module SEq1
using ForwardDiff

"""
    Bisection(f::Function, a::Real, b::Real; TOL::Float64=0.0000001, N::Int=40, output_seq::Bool=false)

Use Bisection To find a solution to.

``f(x) = 0``

You should use the function, and a interval ``[a, b]``, tolerance TOL, maximum number of iterations N.

Then you will get approximate solution ``p`` or can't find in maximum N and return -1.

you can set output_seq=true to output sequence.
"""
@inline function Bisection(f::Function, a::Real, b::Real; TOL::Float64=0.0000001, N::Int=40, output_seq::Bool=false)
    i = 1
    fa = f(a)
    seq = []
    if f(a)*f(b) > 0
        return -1
    else
        while i <= N
            p = a + (b - a)/2
            fp = f(p)
            push!(seq, [p, f(p)])
            if fp == 0 || (b - a)/2 < TOL
                if output_seq
                    return seq
                else
                    return p
                end
                break
            end
            i += 1
            if fa*fp > 0
                a = p
                fa = fp
            else
                b = p
            end
        end
        return -1
    end
end

"""
    fixed_point(f::Function, p₀::Real; TOL::Float64=0.0000001, N::Int=30, output_seq::Bool=false)

Use fixed point Iteration to find root.

The number p is a **fixed point** for a given function g if ``g(p) = p``.

input a function, initial p₀, optional Args: tolerancs TOL, maximum number of iterations N, output_seq.

Output a root p or sequence.
"""
@inline function fixed_point(f::Function, p₀::Real; TOL::Float64=0.0000001, N::Int=30, output_seq::Bool=false)
    i = 1
    seq = []
    while i <= N
        p = f(p₀)
        push!(seq, p₀)
        if abs(p - p₀) < TOL
            if output_seq
                return seq
            else
                return p
            end
            break
        end
        i += 1
        p₀ = p
    end
    return -1
end

"""
    Newton(f::Function, p₀::Real; TOL::Float64=0.0000001, N::Int=20, output_seq=false)

Use Newton method to find root.

you should use a function and initial p₀, then you can set Args or use default.

"""
@inline function Newton(f::Function, p₀::Real; TOL::Float64=0.0000001, N::Int=20, output_seq::Bool=false)
    i = 1
    seq = []
    f1(x) = ForwardDiff.derivative(x -> f(x), x)
    while i <= N
        push!(seq, p₀)
        p = p₀ - f(p₀)/f1(p₀)
        if abs(p - p₀) < TOL
            if output_seq
                return seq
            else
                return p
            end
            break
        end
        p₀ = p
        i += 1
    end
    return -1
end

"""
    Secant(f::Function, p₀::Real, )

Use the Secant Method to find root
A function and p₀, p₁(p₀<p₁), then you can set Args or use default.
"""
@inline function Secant(f::Function, p₀::Real, p₁::Real; TOL::Float64=0.0000001, N::Int=30, output_seq::Bool=false)
    i = 2
    seq = []
    a₀ = f(p₀)
    a₁ = f(p₁)
    push!(seq, p₀)
    push!(seq, p₁)
    while i <= N
        p = p₁ - a₁*(p₁ - p₀)/(a₁ - a₀)
        push!(seq, p)
        if abs(p - p₁) < TOL
            if output_seq
                return seq
            else
                return p
            end
            break
        end
        i += 1
        p₀ = p₁
        a₀ = a₁
        p₁ = p
        a₁ = f(p)
    end
end

"""
    FalsePos(f::Function, p₀::Real, p₁::Real; TOL::Float64=0.00000001, N::Int=40, output_seq::Bool=false)

Use the false position method to find root.
a function and interval ``[a, b]``, then you can set Args or use default.
"""
@inline function FalsePos(f::Function, p₀::Real, p₁::Real; TOL::Float64=0.00000001, N::Int=40, output_seq::Bool=false)
    i = 2
    seq = []
    q₀ = f(p₀)
    q₁ = f(p₁)
    push!(seq, p₀)
    push!(seq, p₁)
    while i <= N
        p = p₁ - q₁*(p₁ - p₀)/(q₁ - q₀)
        push!(seq, p)
        if abs(p - p₁) < TOL
            if output_seq
                return seq
            else
                return p
            end
            break
        end
        i += 1
        q = f(p)
        if q*q₁ < 0
            p₀ = p₁
            q₀ = q₁
        end
        p₁ = p
        q₁ = q
    end
    return -1
end

"""
    ModifiedNewton(f::Function, p₀::Real; TOL::Float64=0.00000001, N::Int=20, output_seq::Bool=false)

Newton and ModifiedNewton, Both methods are rapidly convergent to the actual zero.
"""
@inline function ModifiedNewton(f::Function, p₀::Real; TOL::Float64=0.00000001, N::Int=20, output_seq::Bool=false)
    i = 1
    seq = []
    push!(seq, p₀)
    f1(x) = ForwardDiff.derivative(x -> f(x), x)
    f11(x) = ForwardDiff.derivative(x -> f1(x), x)
    while i <= N
        p = p₀ - f(p₀)*f1(p₀)/((f1(p₀))^2 - f(p₀)*f11(p₀))
        push!(seq, p)
        if abs(p - p₀) < TOL
            if output_seq
                return seq
            else
                return p
            end
            break
        end
        i += 1
        p₀ = p
    end
    return -1
end

"""
    Muller(f::Function, p₀::Real, p₁::Real, p₂::Real; TOL::Float64=0.00000001, N::Int=20, output_seq::Bool=false)
"""
@inline function Muller(f::Function, p₀::Real, p₁::Real, p₂::Real; TOL::Float64=0.00000001, N::Int=20, output_seq::Bool=false)
    h₁ = p₁ - p₀
    h₂ = p₂ - p₁
    δ₁ = (f(p₁) - f(p₀))/h₁
    δ₂ = (f(p₂) - f(p₁))/h₂
    d = (δ₂ - δ₁)/(h₁ + h₂)
    i = 3
    seq = []
    while i <= N
        b = δ₂ + h₂*d
        D = sqrt(complex(b^2 - 4*f(p₂)*d))
        if abs(b - D) < abs(b + D)
            E = b + D
        else
            E = b - D
        end
        h = -2 * f(p₂)/E
        p = p₂ + h
        push!(seq, p)
        if abs(h) < TOL
            if output_seq
                return seq
            else
                return p
            end
            break
        end
        p₀ = p₁
        p₁ = p₂
        p₂ = p
        h₁ = p₁ - p₀
        h₂ = p₂ - p₁
        δ₁ = (f(p₁) - f(p₀))/h₁
        δ₂ = (f(p₂) - f(p₁))/h₂
        d = (δ₂ - δ₁)/(h₁ + h₂)
        i += 1
    end
    return -1
end

end


text: println("Building BookBayes in: ", __DIR__)


text:
using Documenter
using NoiseCC

makedocs(
    modules = [NoiseCC],
    sitename = "NoiseCC",
    authors = "Fu Yin",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = Any[
        "Home" => "index.md",
        "Framwork" => "Framwork.md",
        "FFT" => "FFT.md",
        "CC" => "CC.md",
        "FFT_CC" => "FFT_CC.md",
        "STACK" => "STACK.md",
        "Plot" => "plot.md",
        "References" => "references.md",
        ],
)

# deploydocs(
#     repo = "github.com/OUCyf/NoiseCC.jl.git",
#     target = "build",
#     deps   = nothing,
#     make   = nothing,
# )

max_stars_repo_name: OUCyf/NoiseCC

max_issues_repo_name: OUCyf/NoiseCC

max_forks_repo_name: OUCyf/NoiseCC